Actual source code: tsadapt.c
petsc-3.14.6 2021-03-30
2: #include <petsc/private/tsimpl.h>
4: PetscClassId TSADAPT_CLASSID;
6: static PetscFunctionList TSAdaptList;
7: static PetscBool TSAdaptPackageInitialized;
8: static PetscBool TSAdaptRegisterAllCalled;
10: PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt);
11: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt);
12: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt);
13: PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt);
14: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt);
15: PETSC_EXTERN PetscErrorCode TSAdaptCreate_History(TSAdapt);
17: /*@C
18: TSAdaptRegister - adds a TSAdapt implementation
20: Not Collective
22: Input Parameters:
23: + name_scheme - name of user-defined adaptivity scheme
24: - routine_create - routine to create method context
26: Notes:
27: TSAdaptRegister() may be called multiple times to add several user-defined families.
29: Sample usage:
30: .vb
31: TSAdaptRegister("my_scheme",MySchemeCreate);
32: .ve
34: Then, your scheme can be chosen with the procedural interface via
35: $ TSAdaptSetType(ts,"my_scheme")
36: or at runtime via the option
37: $ -ts_adapt_type my_scheme
39: Level: advanced
41: .seealso: TSAdaptRegisterAll()
42: @*/
43: PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
44: {
48: TSAdaptInitializePackage();
49: PetscFunctionListAdd(&TSAdaptList,sname,function);
50: return(0);
51: }
53: /*@C
54: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
56: Not Collective
58: Level: advanced
60: .seealso: TSAdaptRegisterDestroy()
61: @*/
62: PetscErrorCode TSAdaptRegisterAll(void)
63: {
67: if (TSAdaptRegisterAllCalled) return(0);
68: TSAdaptRegisterAllCalled = PETSC_TRUE;
69: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
70: TSAdaptRegister(TSADAPTBASIC, TSAdaptCreate_Basic);
71: TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP);
72: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
73: TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);
74: TSAdaptRegister(TSADAPTHISTORY,TSAdaptCreate_History);
75: return(0);
76: }
78: /*@C
79: TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
80: called from PetscFinalize().
82: Level: developer
84: .seealso: PetscFinalize()
85: @*/
86: PetscErrorCode TSAdaptFinalizePackage(void)
87: {
91: PetscFunctionListDestroy(&TSAdaptList);
92: TSAdaptPackageInitialized = PETSC_FALSE;
93: TSAdaptRegisterAllCalled = PETSC_FALSE;
94: return(0);
95: }
97: /*@C
98: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
99: called from TSInitializePackage().
101: Level: developer
103: .seealso: PetscInitialize()
104: @*/
105: PetscErrorCode TSAdaptInitializePackage(void)
106: {
110: if (TSAdaptPackageInitialized) return(0);
111: TSAdaptPackageInitialized = PETSC_TRUE;
112: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
113: TSAdaptRegisterAll();
114: PetscRegisterFinalize(TSAdaptFinalizePackage);
115: return(0);
116: }
118: /*@C
119: TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE
121: Logicially Collective on TSAdapt
123: Input Parameter:
124: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
125: - type - either TSADAPTBASIC or TSADAPTNONE
127: Options Database:
128: . -ts_adapt_type <basic or dsp or none> - to set the adapter type
130: Level: intermediate
132: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
133: @*/
134: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
135: {
136: PetscBool match;
137: PetscErrorCode ierr,(*r)(TSAdapt);
142: PetscObjectTypeCompare((PetscObject)adapt,type,&match);
143: if (match) return(0);
144: PetscFunctionListFind(TSAdaptList,type,&r);
145: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
146: if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
147: PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
148: PetscObjectChangeTypeName((PetscObject)adapt,type);
149: (*r)(adapt);
150: return(0);
151: }
153: /*@C
154: TSAdaptGetType - gets the TS adapter method type (as a string).
156: Not Collective
158: Input Parameter:
159: . adapt - The TS adapter, most likely obtained with TSGetAdapt()
161: Output Parameter:
162: . type - The name of TS adapter method
164: Level: intermediate
166: .seealso TSAdaptSetType()
167: @*/
168: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
169: {
173: *type = ((PetscObject)adapt)->type_name;
174: return(0);
175: }
177: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
178: {
183: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
184: return(0);
185: }
187: /*@C
188: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
190: Collective on PetscViewer
192: Input Parameters:
193: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
194: some related function before a call to TSAdaptLoad().
195: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
196: HDF5 file viewer, obtained from PetscViewerHDF5Open()
198: Level: intermediate
200: Notes:
201: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
203: Notes for advanced users:
204: Most users should not need to know the details of the binary storage
205: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
206: But for anyone who's interested, the standard binary matrix storage
207: format is
208: .vb
209: has not yet been determined
210: .ve
212: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
213: @*/
214: PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
215: {
217: PetscBool isbinary;
218: char type[256];
223: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
224: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
226: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
227: TSAdaptSetType(adapt,type);
228: if (adapt->ops->load) {
229: (*adapt->ops->load)(adapt,viewer);
230: }
231: return(0);
232: }
234: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
235: {
237: PetscBool iascii,isbinary,isnone,isglee;
241: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
244: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
245: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
246: if (iascii) {
247: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
248: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
249: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTGLEE,&isglee);
250: if (!isnone) {
251: if (adapt->always_accept) {PetscViewerASCIIPrintf(viewer," always accepting steps\n");}
252: PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);
253: PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
254: PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);
255: PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);
256: PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);
257: PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);
258: PetscViewerASCIIPrintf(viewer," maximum solution absolute value to be ignored %g\n",(double)adapt->ignore_max);
259: }
260: if (isglee) {
261: if (adapt->glee_use_local) {
262: PetscViewerASCIIPrintf(viewer," GLEE uses local error control\n");
263: } else {
264: PetscViewerASCIIPrintf(viewer," GLEE uses global error control\n");
265: }
266: }
267: if (adapt->ops->view) {
268: PetscViewerASCIIPushTab(viewer);
269: (*adapt->ops->view)(adapt,viewer);
270: PetscViewerASCIIPopTab(viewer);
271: }
272: } else if (isbinary) {
273: char type[256];
275: /* need to save FILE_CLASS_ID for adapt class */
276: PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
277: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
278: } else if (adapt->ops->view) {
279: (*adapt->ops->view)(adapt,viewer);
280: }
281: return(0);
282: }
284: /*@
285: TSAdaptReset - Resets a TSAdapt context.
287: Collective on TS
289: Input Parameter:
290: . adapt - the TSAdapt context obtained from TSAdaptCreate()
292: Level: developer
294: .seealso: TSAdaptCreate(), TSAdaptDestroy()
295: @*/
296: PetscErrorCode TSAdaptReset(TSAdapt adapt)
297: {
302: if (adapt->ops->reset) {(*adapt->ops->reset)(adapt);}
303: return(0);
304: }
306: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
307: {
311: if (!*adapt) return(0);
313: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return(0);}
315: TSAdaptReset(*adapt);
317: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
318: PetscViewerDestroy(&(*adapt)->monitor);
319: PetscHeaderDestroy(adapt);
320: return(0);
321: }
323: /*@
324: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
326: Collective on TSAdapt
328: Input Arguments:
329: + adapt - adaptive controller context
330: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
332: Options Database Keys:
333: . -ts_adapt_monitor - to turn on monitoring
335: Level: intermediate
337: .seealso: TSAdaptChoose()
338: @*/
339: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
340: {
346: if (flg) {
347: if (!adapt->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);}
348: } else {
349: PetscViewerDestroy(&adapt->monitor);
350: }
351: return(0);
352: }
354: /*@C
355: TSAdaptSetCheckStage - Set a callback to check convergence for a stage
357: Logically collective on TSAdapt
359: Input Arguments:
360: + adapt - adaptive controller context
361: - func - stage check function
363: Arguments of func:
364: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
366: + adapt - adaptive controller context
367: . ts - time stepping context
368: - accept - pending choice of whether to accept, can be modified by this routine
370: Level: advanced
372: .seealso: TSAdaptChoose()
373: @*/
374: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*))
375: {
379: adapt->checkstage = func;
380: return(0);
381: }
383: /*@
384: TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of
385: any error or stability condition not meeting the prescribed goal.
387: Logically collective on TSAdapt
389: Input Arguments:
390: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
391: - flag - whether to always accept steps
393: Options Database Keys:
394: . -ts_adapt_always_accept - to always accept steps
396: Level: intermediate
398: .seealso: TSAdapt, TSAdaptChoose()
399: @*/
400: PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)
401: {
405: adapt->always_accept = flag;
406: return(0);
407: }
409: /*@
410: TSAdaptSetSafety - Set safety factors
412: Logically collective on TSAdapt
414: Input Arguments:
415: + adapt - adaptive controller context
416: . safety - safety factor relative to target error/stability goal
417: - reject_safety - extra safety factor to apply if the last step was rejected
419: Options Database Keys:
420: + -ts_adapt_safety <safety> - to set safety factor
421: - -ts_adapt_reject_safety <reject_safety> - to set reject safety factor
423: Level: intermediate
425: .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose()
426: @*/
427: PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)
428: {
433: if (safety != PETSC_DEFAULT && safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be non negative",(double)safety);
434: if (safety != PETSC_DEFAULT && safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be less than one",(double)safety);
435: if (reject_safety != PETSC_DEFAULT && reject_safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reject safety factor %g must be non negative",(double)reject_safety);
436: if (reject_safety != PETSC_DEFAULT && reject_safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reject safety factor %g must be less than one",(double)reject_safety);
437: if (safety != PETSC_DEFAULT) adapt->safety = safety;
438: if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety;
439: return(0);
440: }
442: /*@
443: TSAdaptGetSafety - Get safety factors
445: Not Collective
447: Input Arguments:
448: . adapt - adaptive controller context
450: Ouput Arguments:
451: . safety - safety factor relative to target error/stability goal
452: + reject_safety - extra safety factor to apply if the last step was rejected
454: Level: intermediate
456: .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose()
457: @*/
458: PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety)
459: {
464: if (safety) *safety = adapt->safety;
465: if (reject_safety) *reject_safety = adapt->reject_safety;
466: return(0);
467: }
469: /*@
470: TSAdaptSetMaxIgnore - Set error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value). A negative value (default) of the threshold leads to considering all solution components.
472: Logically collective on TSAdapt
474: Input Arguments:
475: + adapt - adaptive controller context
476: - max_ignore - threshold for solution components that are ignored during error estimation
478: Options Database Keys:
479: . -ts_adapt_max_ignore <max_ignore> - to set the threshold
481: Level: intermediate
483: .seealso: TSAdapt, TSAdaptGetMaxIgnore(), TSAdaptChoose()
484: @*/
485: PetscErrorCode TSAdaptSetMaxIgnore(TSAdapt adapt,PetscReal max_ignore)
486: {
490: adapt->ignore_max = max_ignore;
491: return(0);
492: }
494: /*@
495: TSAdaptGetMaxIgnore - Get error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value).
497: Not Collective
499: Input Arguments:
500: . adapt - adaptive controller context
502: Ouput Arguments:
503: . max_ignore - threshold for solution components that are ignored during error estimation
505: Level: intermediate
507: .seealso: TSAdapt, TSAdaptSetMaxIgnore(), TSAdaptChoose()
508: @*/
509: PetscErrorCode TSAdaptGetMaxIgnore(TSAdapt adapt,PetscReal *max_ignore)
510: {
514: *max_ignore = adapt->ignore_max;
515: return(0);
516: }
519: /*@
520: TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
522: Logically collective on TSAdapt
524: Input Arguments:
525: + adapt - adaptive controller context
526: . low - admissible decrease factor
527: - high - admissible increase factor
529: Options Database Keys:
530: . -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors
532: Level: intermediate
534: .seealso: TSAdaptChoose(), TSAdaptGetClip(), TSAdaptSetScaleSolveFailed()
535: @*/
536: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
537: {
542: if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
543: if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
544: if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be greater than one",(double)high);
545: if (low != PETSC_DEFAULT) adapt->clip[0] = low;
546: if (high != PETSC_DEFAULT) adapt->clip[1] = high;
547: return(0);
548: }
550: /*@
551: TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
553: Not Collective
555: Input Arguments:
556: . adapt - adaptive controller context
558: Ouput Arguments:
559: + low - optional, admissible decrease factor
560: - high - optional, admissible increase factor
562: Level: intermediate
564: .seealso: TSAdaptChoose(), TSAdaptSetClip(), TSAdaptSetScaleSolveFailed()
565: @*/
566: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
567: {
572: if (low) *low = adapt->clip[0];
573: if (high) *high = adapt->clip[1];
574: return(0);
575: }
577: /*@
578: TSAdaptSetScaleSolveFailed - Scale step by this factor if solve fails
580: Logically collective on TSAdapt
582: Input Arguments:
583: + adapt - adaptive controller context
584: - scale - scale
586: Options Database Keys:
587: . -ts_adapt_scale_solve_failed <scale> - to set scale step by this factor if solve fails
589: Level: intermediate
591: .seealso: TSAdaptChoose(), TSAdaptGetScaleSolveFailed(), TSAdaptGetClip()
592: @*/
593: PetscErrorCode TSAdaptSetScaleSolveFailed(TSAdapt adapt,PetscReal scale)
594: {
598: if (scale != PETSC_DEFAULT && scale <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be positive",(double)scale);
599: if (scale != PETSC_DEFAULT && scale > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be less than one",(double)scale);
600: if (scale != PETSC_DEFAULT) adapt->scale_solve_failed = scale;
601: return(0);
602: }
604: /*@
605: TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size
607: Not Collective
609: Input Arguments:
610: . adapt - adaptive controller context
612: Ouput Arguments:
613: . scale - scale factor
615: Level: intermediate
617: .seealso: TSAdaptChoose(), TSAdaptSetScaleSolveFailed(), TSAdaptSetClip()
618: @*/
619: PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal *scale)
620: {
624: if (scale) *scale = adapt->scale_solve_failed;
625: return(0);
626: }
628: /*@
629: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
631: Logically collective on TSAdapt
633: Input Arguments:
634: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
635: . hmin - minimum time step
636: - hmax - maximum time step
638: Options Database Keys:
639: + -ts_adapt_dt_min <min> - to set minimum time step
640: - -ts_adapt_dt_max <max> - to set maximum time step
642: Level: intermediate
644: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
645: @*/
646: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
647: {
653: if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
654: if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
655: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
656: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
657: hmin = adapt->dt_min;
658: hmax = adapt->dt_max;
659: if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must greater than minimum time step %g",(double)hmax,(double)hmin);
660: return(0);
661: }
663: /*@
664: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
666: Not Collective
668: Input Arguments:
669: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
671: Output Arguments:
672: + hmin - minimum time step
673: - hmax - maximum time step
675: Level: intermediate
677: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
678: @*/
679: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
680: {
686: if (hmin) *hmin = adapt->dt_min;
687: if (hmax) *hmax = adapt->dt_max;
688: return(0);
689: }
691: /*
692: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
694: Collective on TSAdapt
696: Input Parameter:
697: . adapt - the TSAdapt context
699: Options Database Keys:
700: + -ts_adapt_type <type> - algorithm to use for adaptivity
701: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
702: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
703: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
704: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
705: . -ts_adapt_dt_min <min> - minimum timestep to use
706: . -ts_adapt_dt_max <max> - maximum timestep to use
707: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
708: . -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
709: - -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver
711: Level: advanced
713: Notes:
714: This function is automatically called by TSSetFromOptions()
716: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
717: TSAdaptSetClip(), TSAdaptSetScaleSolveFailed(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
718: */
719: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
720: {
722: char type[256] = TSADAPTBASIC;
723: PetscReal safety,reject_safety,clip[2],scale,hmin,hmax;
724: PetscBool set,flg;
725: PetscInt two;
729: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
730: * function can only be called from inside TSSetFromOptions() */
731: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
732: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
733: if (flg || !((PetscObject)adapt)->type_name) {
734: TSAdaptSetType(adapt,type);
735: }
737: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
738: if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}
740: safety = adapt->safety; reject_safety = adapt->reject_safety;
741: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
742: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
743: if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}
745: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
746: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
747: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
748: if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}
750: hmin = adapt->dt_min; hmax = adapt->dt_max;
751: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
752: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
753: if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}
755: PetscOptionsReal("-ts_adapt_max_ignore","Adaptor ignores (absolute) solution values smaller than this value","",adapt->ignore_max,&adapt->ignore_max,&set);
756: PetscOptionsBool("-ts_adapt_glee_use_local","GLEE adaptor uses local error estimation for step control","",adapt->glee_use_local,&adapt->glee_use_local,&set);
758: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","TSAdaptSetScaleSolveFailed",adapt->scale_solve_failed,&scale,&set);
759: if (set) {TSAdaptSetScaleSolveFailed(adapt,scale);}
761: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
762: if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
764: PetscOptionsInt("-ts_adapt_time_step_increase_delay","Number of timesteps to delay increasing the time step after it has been decreased due to failed solver","TSAdaptSetTimeStepIncreaseDelay",adapt->timestepjustdecreased_delay,&adapt->timestepjustdecreased_delay,NULL);
766: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
767: if (set) {TSAdaptSetMonitor(adapt,flg);}
769: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
770: PetscOptionsTail();
771: return(0);
772: }
774: /*@
775: TSAdaptCandidatesClear - clear any previously set candidate schemes
777: Logically collective on TSAdapt
779: Input Argument:
780: . adapt - adaptive controller
782: Level: developer
784: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
785: @*/
786: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
787: {
792: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
793: return(0);
794: }
796: /*@C
797: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
799: Logically collective on TSAdapt
801: Input Arguments:
802: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
803: . name - name of the candidate scheme to add
804: . order - order of the candidate scheme
805: . stageorder - stage order of the candidate scheme
806: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
807: . cost - relative measure of the amount of work required for the candidate scheme
808: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
810: Note:
811: This routine is not available in Fortran.
813: Level: developer
815: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
816: @*/
817: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
818: {
819: PetscInt c;
823: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
824: if (inuse) {
825: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
826: adapt->candidates.inuse_set = PETSC_TRUE;
827: }
828: /* first slot if this is the current scheme, otherwise the next available slot */
829: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
831: adapt->candidates.name[c] = name;
832: adapt->candidates.order[c] = order;
833: adapt->candidates.stageorder[c] = stageorder;
834: adapt->candidates.ccfl[c] = ccfl;
835: adapt->candidates.cost[c] = cost;
836: adapt->candidates.n++;
837: return(0);
838: }
840: /*@C
841: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
843: Not Collective
845: Input Arguments:
846: . adapt - time step adaptivity context
848: Output Arguments:
849: + n - number of candidate schemes, always at least 1
850: . order - the order of each candidate scheme
851: . stageorder - the stage order of each candidate scheme
852: . ccfl - the CFL coefficient of each scheme
853: - cost - the relative cost of each scheme
855: Level: developer
857: Note:
858: The current scheme is always returned in the first slot
860: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
861: @*/
862: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
863: {
866: if (n) *n = adapt->candidates.n;
867: if (order) *order = adapt->candidates.order;
868: if (stageorder) *stageorder = adapt->candidates.stageorder;
869: if (ccfl) *ccfl = adapt->candidates.ccfl;
870: if (cost) *cost = adapt->candidates.cost;
871: return(0);
872: }
874: /*@C
875: TSAdaptChoose - choose which method and step size to use for the next step
877: Collective on TSAdapt
879: Input Arguments:
880: + adapt - adaptive contoller
881: - h - current step size
883: Output Arguments:
884: + next_sc - optional, scheme to use for the next step
885: . next_h - step size to use for the next step
886: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
888: Note:
889: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
890: being retried after an initial rejection.
892: Level: developer
894: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
895: @*/
896: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
897: {
899: PetscInt ncandidates = adapt->candidates.n;
900: PetscInt scheme = 0;
901: PetscReal wlte = -1.0;
902: PetscReal wltea = -1.0;
903: PetscReal wlter = -1.0;
911: if (next_sc) *next_sc = 0;
913: /* Do not mess with adaptivity while handling events*/
914: if (ts->event && ts->event->status != TSEVENT_NONE) {
915: *next_h = h;
916: *accept = PETSC_TRUE;
917: return(0);
918: }
920: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
921: if (scheme < 0 || (ncandidates > 0 && scheme >= ncandidates)) SETERRQ2(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",scheme,ncandidates-1);
922: if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
923: if (next_sc) *next_sc = scheme;
925: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
926: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
927: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
928: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
929: PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
930: PetscReal b = adapt->matchstepfac[1];
931: if (t < tmax && tend > tmax) *next_h = hmax;
932: if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
933: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
934: }
936: if (adapt->monitor) {
937: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
938: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
939: if (wlte < 0) {
940: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h);
941: } else {
942: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e wlte=%5.3g wltea=%5.3g wlter=%5.3g\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h,(double)wlte,(double)wltea,(double)wlter);
943: }
944: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
945: }
946: return(0);
947: }
949: /*@
950: TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
951: before increasing the time step.
953: Logicially Collective on TSAdapt
955: Input Arguments:
956: + adapt - adaptive controller context
957: - cnt - the number of timesteps
959: Options Database Key:
960: . -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase
962: Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
963: The successful use of this option is problem dependent
965: Developer Note: there is no theory to support this option
967: Level: advanced
969: .seealso:
970: @*/
971: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
972: {
974: adapt->timestepjustdecreased_delay = cnt;
975: return(0);
976: }
978: /*@
979: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible)
981: Collective on TSAdapt
983: Input Arguments:
984: + adapt - adaptive controller context
985: . ts - time stepper
986: . t - Current simulation time
987: - Y - Current solution vector
989: Output Arguments:
990: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
992: Level: developer
994: .seealso:
995: @*/
996: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
997: {
998: PetscErrorCode ierr;
999: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
1006: if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
1007: if (snesreason < 0) {
1008: *accept = PETSC_FALSE;
1009: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
1010: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
1011: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
1012: if (adapt->monitor) {
1013: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1014: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected t=%-11g+%10.3e, nonlinear solve failures %D greater than current TS allowed\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)ts->time_step,ts->num_snes_failures);
1015: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1016: }
1017: }
1018: } else {
1019: *accept = PETSC_TRUE;
1020: TSFunctionDomainError(ts,t,Y,accept);
1021: if (*accept && adapt->checkstage) {
1022: (*adapt->checkstage)(adapt,ts,t,Y,accept);
1023: if (!*accept) {
1024: PetscInfo1(ts,"Step=%D, solution rejected by user function provided by TSSetFunctionDomainError()\n",ts->steps);
1025: if (adapt->monitor) {
1026: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1027: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected by user function provided by TSSetFunctionDomainError()\n",((PetscObject)adapt)->type_name,ts->steps);
1028: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1029: }
1030: }
1031: }
1032: }
1034: if (!(*accept) && !ts->reason) {
1035: PetscReal dt,new_dt;
1036: TSGetTimeStep(ts,&dt);
1037: new_dt = dt * adapt->scale_solve_failed;
1038: TSSetTimeStep(ts,new_dt);
1039: adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
1040: if (adapt->monitor) {
1041: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1042: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected (%s) t=%-11g+%10.3e retrying with dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,SNESConvergedReasons[snesreason],(double)ts->ptime,(double)dt,(double)new_dt);
1043: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1044: }
1045: }
1046: return(0);
1047: }
1049: /*@
1050: TSAdaptCreate - create an adaptive controller context for time stepping
1052: Collective
1054: Input Parameter:
1055: . comm - The communicator
1057: Output Parameter:
1058: . adapt - new TSAdapt object
1060: Level: developer
1062: Notes:
1063: TSAdapt creation is handled by TS, so users should not need to call this function.
1065: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1066: @*/
1067: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1068: {
1070: TSAdapt adapt;
1074: *inadapt = NULL;
1075: TSAdaptInitializePackage();
1077: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
1079: adapt->always_accept = PETSC_FALSE;
1080: adapt->safety = 0.9;
1081: adapt->reject_safety = 0.5;
1082: adapt->clip[0] = 0.1;
1083: adapt->clip[1] = 10.;
1084: adapt->dt_min = 1e-20;
1085: adapt->dt_max = 1e+20;
1086: adapt->ignore_max = -1.0;
1087: adapt->glee_use_local = PETSC_TRUE;
1088: adapt->scale_solve_failed = 0.25;
1089: /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1090: to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1091: adapt->matchstepfac[0] = 0.01; /* allow 1% step size increase in the last step */
1092: adapt->matchstepfac[1] = 2.0; /* halve last step if it is greater than what remains divided this factor */
1093: adapt->wnormtype = NORM_2;
1094: adapt->timestepjustdecreased_delay = 0;
1096: *inadapt = adapt;
1097: return(0);
1098: }