1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 typedef struct { 4 PetscBool always_accept; 5 PetscReal clip[2]; /* admissible decrease/increase factors */ 6 PetscReal safety; /* safety factor relative to target error */ 7 PetscReal reject_safety; /* extra safety factor if the last step was rejected */ 8 Vec Y; 9 } TSAdapt_Basic; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "TSAdaptChoose_Basic" 13 static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte) 14 { 15 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 16 PetscInt order = PETSC_DECIDE; 17 PetscReal enorm = -1; 18 PetscReal enorma,enormr; 19 PetscReal safety = basic->safety; 20 PetscReal hfac_lte,h_lte; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 *next_sc = 0; /* Reuse the same order scheme */ 25 26 if (ts->ops->evaluatewlte) { 27 ierr = TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);CHKERRQ(ierr); 28 if (enorm >= 0 && order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order); 29 } else if (ts->ops->evaluatestep) { 30 if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered"); 31 if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n); 32 if (!basic->Y) {ierr = VecDuplicate(ts->vec_sol,&basic->Y);CHKERRQ(ierr);} 33 order = adapt->candidates.order[0]; 34 ierr = TSEvaluateStep(ts,order-1,basic->Y,NULL);CHKERRQ(ierr); 35 ierr = TSErrorWeightedNorm(ts,ts->vec_sol,basic->Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 36 } 37 38 if (enorm < 0) { 39 *accept = PETSC_TRUE; 40 *next_h = h; /* Reuse the old step */ 41 *wlte = -1; /* Weighted local truncation error was not evaluated */ 42 PetscFunctionReturn(0); 43 } 44 45 /* Determine whether the step is accepted of rejected */ 46 if (enorm > 1) { 47 if (!*accept) safety *= basic->reject_safety; /* The last attempt also failed, shorten more aggressively */ 48 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 49 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);CHKERRQ(ierr); 50 *accept = PETSC_TRUE; 51 } else if (basic->always_accept) { 52 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);CHKERRQ(ierr); 53 *accept = PETSC_TRUE; 54 } else { 55 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr); 56 *accept = PETSC_FALSE; 57 } 58 } else { 59 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr); 60 *accept = PETSC_TRUE; 61 } 62 63 /* The optimal new step based purely on local truncation error for this step. */ 64 if (enorm > 0) 65 hfac_lte = safety * PetscPowReal(enorm,((PetscReal)-1)/order); 66 else 67 hfac_lte = safety * PETSC_INFINITY; 68 h_lte = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]); 69 70 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 71 *wlte = enorm; 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "TSAdaptReset_Basic" 77 static PetscErrorCode TSAdaptReset_Basic(TSAdapt adapt) 78 { 79 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 ierr = VecDestroy(&basic->Y);CHKERRQ(ierr); 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "TSAdaptDestroy_Basic" 89 static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt) 90 { 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 ierr = TSAdaptReset_Basic(adapt);CHKERRQ(ierr); 95 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 #undef __FUNCT__ 100 #define __FUNCT__ "TSAdaptSetFromOptions_Basic" 101 static PetscErrorCode TSAdaptSetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSAdapt adapt) 102 { 103 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 104 PetscErrorCode ierr; 105 PetscInt two; 106 PetscBool set; 107 108 PetscFunctionBegin; 109 ierr = PetscOptionsHead(PetscOptionsObject,"Basic adaptive controller options");CHKERRQ(ierr); 110 two = 2; 111 ierr = PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase factor in step size","TSAdaptBasicSetClip",basic->clip,&two,&set);CHKERRQ(ierr); 112 if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip"); 113 if (set) {ierr = TSAdaptBasicSetClip(adapt,basic->clip[0],basic->clip[1]);CHKERRQ(ierr);} 114 ierr = PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,NULL);CHKERRQ(ierr); 115 ierr = PetscOptionsReal("-ts_adapt_basic_reject_safety","Extra safety factor to apply if the last step was rejected","",basic->reject_safety,&basic->reject_safety,NULL);CHKERRQ(ierr); 116 ierr = PetscOptionsBool("-ts_adapt_basic_always_accept","Always accept the step regardless of whether local truncation error meets goal","",basic->always_accept,&basic->always_accept,NULL);CHKERRQ(ierr); 117 ierr = PetscOptionsTail();CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "TSAdaptView_Basic" 123 static PetscErrorCode TSAdaptView_Basic(TSAdapt adapt,PetscViewer viewer) 124 { 125 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 126 PetscErrorCode ierr; 127 PetscBool iascii; 128 129 PetscFunctionBegin; 130 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 131 if (iascii) { 132 if (basic->always_accept) {ierr = PetscViewerASCIIPrintf(viewer," Basic: always accepting steps\n");CHKERRQ(ierr);} 133 ierr = PetscViewerASCIIPrintf(viewer," Basic: clip fastest decrease %g, fastest increase %g\n",(double)basic->clip[0],(double)basic->clip[1]);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(viewer," Basic: safety factor %g, extra factor after step rejection %g\n",(double)basic->safety,(double)basic->reject_safety);CHKERRQ(ierr); 135 } 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNCT__ 140 #define __FUNCT__ "TSAdaptCreate_Basic" 141 /*MC 142 TSADAPTBASIC - Basic adaptive controller for time stepping 143 144 Level: intermediate 145 146 .seealso: TS, TSAdapt, TSSetAdapt() 147 M*/ 148 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 149 { 150 PetscErrorCode ierr; 151 TSAdapt_Basic *a; 152 153 PetscFunctionBegin; 154 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr); 155 adapt->data = (void*)a; 156 adapt->ops->choose = TSAdaptChoose_Basic; 157 adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic; 158 adapt->ops->destroy = TSAdaptDestroy_Basic; 159 adapt->ops->view = TSAdaptView_Basic; 160 161 a->clip[0] = 0.1; 162 a->clip[1] = 10.; 163 a->safety = 0.9; 164 a->reject_safety = 0.5; 165 a->always_accept = PETSC_FALSE; 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "TSAdaptBasicSetClip" 171 /*@ 172 TSAdaptBasicSetClip - Sets the admissible decrease/increase factor in step size 173 174 Collective on TSAdapt 175 176 Input Arguments: 177 + adapt - adaptive controller context 178 . low - admissible decrease factor 179 + high - admissible increase factor 180 181 Level: intermediate 182 183 .seealso: TSAdaptChoose() 184 @*/ 185 PetscErrorCode TSAdaptBasicSetClip(TSAdapt adapt,PetscReal low,PetscReal high) 186 { 187 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 188 PetscBool isbasic; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 193 if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low); 194 if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low); 195 if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high); 196 ierr = PetscObjectTypeCompare((PetscObject)adapt,TSADAPTBASIC,&isbasic);CHKERRQ(ierr); 197 if (isbasic) { 198 if (low != PETSC_DEFAULT) basic->clip[0] = low; 199 if (high != PETSC_DEFAULT) basic->clip[1] = high; 200 } 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "TSAdaptBasicGetClip" 206 /*@ 207 TSAdaptBasicGetClip - Gets the admissible decrease/increase factor in step size 208 209 Collective on TSAdapt 210 211 Input Arguments: 212 . adapt - adaptive controller context 213 214 Ouput Arguments: 215 + low - optional, admissible decrease factor 216 - high - optional, admissible increase factor 217 218 Level: intermediate 219 220 .seealso: TSAdaptChoose() 221 @*/ 222 PetscErrorCode TSAdaptBasicGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high) 223 { 224 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 225 PetscBool isbasic; 226 PetscErrorCode ierr; 227 228 PetscFunctionBegin; 229 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 230 if (low) PetscValidRealPointer(low,2); 231 if (high) PetscValidRealPointer(high,3); 232 ierr = PetscObjectTypeCompare((PetscObject)adapt,TSADAPTBASIC,&isbasic);CHKERRQ(ierr); 233 if (isbasic) { 234 if (low) *low = basic->clip[0]; 235 if (high) *high = basic->clip[1]; 236 } else { 237 if (low) *low = 0; 238 if (high) *high = PETSC_MAX_REAL; 239 } 240 PetscFunctionReturn(0); 241 } 242