Actual source code: adaptbasic.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
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;
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 safety = basic->safety;
19: PetscReal hfac_lte,h_lte;
23: *next_sc = 0; /* Reuse the same order scheme */
25: if (ts->ops->evaluatewlte) {
26: TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);
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) {VecDuplicate(ts->vec_sol,&basic->Y);}
32: order = adapt->candidates.order[0];
33: TSEvaluateStep(ts,order-1,basic->Y,NULL);
34: TSErrorWeightedNorm(ts,ts->vec_sol,basic->Y,adapt->wnormtype,&enorm);
35: }
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: return(0);
42: }
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: PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);
49: *accept = PETSC_TRUE;
50: } else if (basic->always_accept) {
51: PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);
52: *accept = PETSC_TRUE;
53: } else {
54: PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);
55: *accept = PETSC_FALSE;
56: }
57: } else {
58: PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);
59: *accept = PETSC_TRUE;
60: }
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]);
69: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
70: *wlte = enorm;
71: return(0);
72: }
76: static PetscErrorCode TSAdaptReset_Basic(TSAdapt adapt)
77: {
78: TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
82: VecDestroy(&basic->Y);
83: return(0);
84: }
88: static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt)
89: {
93: TSAdaptReset_Basic(adapt);
94: PetscFree(adapt->data);
95: return(0);
96: }
100: static PetscErrorCode TSAdaptSetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
101: {
102: TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
104: PetscInt two;
105: PetscBool set;
108: PetscOptionsHead(PetscOptionsObject,"Basic adaptive controller options");
109: two = 2;
110: PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase factor in step size","TSAdaptBasicSetClip",basic->clip,&two,&set);
111: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip");
112: if (set) {TSAdaptBasicSetClip(adapt,basic->clip[0],basic->clip[1]);}
113: PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,NULL);
114: PetscOptionsReal("-ts_adapt_basic_reject_safety","Extra safety factor to apply if the last step was rejected","",basic->reject_safety,&basic->reject_safety,NULL);
115: 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);
116: PetscOptionsTail();
117: return(0);
118: }
122: static PetscErrorCode TSAdaptView_Basic(TSAdapt adapt,PetscViewer viewer)
123: {
124: TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
126: PetscBool iascii;
129: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
130: if (iascii) {
131: if (basic->always_accept) {PetscViewerASCIIPrintf(viewer," Basic: always accepting steps\n");}
132: PetscViewerASCIIPrintf(viewer," Basic: clip fastest decrease %g, fastest increase %g\n",(double)basic->clip[0],(double)basic->clip[1]);
133: PetscViewerASCIIPrintf(viewer," Basic: safety factor %g, extra factor after step rejection %g\n",(double)basic->safety,(double)basic->reject_safety);
134: }
135: return(0);
136: }
140: /*MC
141: TSADAPTBASIC - Basic adaptive controller for time stepping
143: Level: intermediate
145: .seealso: TS, TSAdapt, TSSetAdapt()
146: M*/
147: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
148: {
150: TSAdapt_Basic *a;
153: PetscNewLog(adapt,&a);
154: adapt->data = (void*)a;
155: adapt->ops->choose = TSAdaptChoose_Basic;
156: adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic;
157: adapt->ops->destroy = TSAdaptDestroy_Basic;
158: adapt->ops->view = TSAdaptView_Basic;
160: a->clip[0] = 0.1;
161: a->clip[1] = 10.;
162: a->safety = 0.9;
163: a->reject_safety = 0.5;
164: a->always_accept = PETSC_FALSE;
165: return(0);
166: }
170: /*@
171: TSAdaptBasicSetClip - Sets the admissible decrease/increase factor in step size
173: Collective on TSAdapt
175: Input Arguments:
176: + adapt - adaptive controller context
177: . low - admissible decrease factor
178: + high - admissible increase factor
180: Level: intermediate
182: .seealso: TSAdaptChoose()
183: @*/
184: PetscErrorCode TSAdaptBasicSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
185: {
186: TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
187: PetscBool isbasic;
192: if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
193: if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
194: if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high);
195: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTBASIC,&isbasic);
196: if (isbasic) {
197: if (low != PETSC_DEFAULT) basic->clip[0] = low;
198: if (high != PETSC_DEFAULT) basic->clip[1] = high;
199: }
200: return(0);
201: }
205: /*@
206: TSAdaptBasicGetClip - Gets the admissible decrease/increase factor in step size
208: Collective on TSAdapt
210: Input Arguments:
211: . adapt - adaptive controller context
213: Ouput Arguments:
214: + low - optional, admissible decrease factor
215: - high - optional, admissible increase factor
217: Level: intermediate
219: .seealso: TSAdaptChoose()
220: @*/
221: PetscErrorCode TSAdaptBasicGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
222: {
223: TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
224: PetscBool isbasic;
231: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTBASIC,&isbasic);
232: if (isbasic) {
233: if (low) *low = basic->clip[0];
234: if (high) *high = basic->clip[1];
235: } else {
236: if (low) *low = 0;
237: if (high) *high = PETSC_MAX_REAL;
238: }
239: return(0);
240: }