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