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