Actual source code: tsadapt.c
petsc-3.12.5 2020-03-29
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,PETSC_FALSE);
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()
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 geather 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()
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: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
580: Logically collective on TSAdapt
582: Input Arguments:
583: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
584: . hmin - minimum time step
585: - hmax - maximum time step
587: Options Database Keys:
588: + -ts_adapt_dt_min <min> - to set minimum time step
589: - -ts_adapt_dt_max <max> - to set maximum time step
591: Level: intermediate
593: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
594: @*/
595: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
596: {
602: if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
603: if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
604: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
605: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
606: hmin = adapt->dt_min;
607: hmax = adapt->dt_max;
608: if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must geather than minimum time step %g",(double)hmax,(double)hmin);
609: return(0);
610: }
612: /*@
613: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
615: Not Collective
617: Input Arguments:
618: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
620: Output Arguments:
621: + hmin - minimum time step
622: - hmax - maximum time step
624: Level: intermediate
626: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
627: @*/
628: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
629: {
635: if (hmin) *hmin = adapt->dt_min;
636: if (hmax) *hmax = adapt->dt_max;
637: return(0);
638: }
640: /*
641: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
643: Collective on TSAdapt
645: Input Parameter:
646: . adapt - the TSAdapt context
648: Options Database Keys:
649: + -ts_adapt_type <type> - algorithm to use for adaptivity
650: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
651: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
652: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
653: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
654: . -ts_adapt_dt_min <min> - minimum timestep to use
655: . -ts_adapt_dt_max <max> - maximum timestep to use
656: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
657: . -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
658: - -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver
660: Level: advanced
662: Notes:
663: This function is automatically called by TSSetFromOptions()
665: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
666: TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
667: */
668: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
669: {
671: char type[256] = TSADAPTBASIC;
672: PetscReal safety,reject_safety,clip[2],hmin,hmax;
673: PetscBool set,flg;
674: PetscInt two;
678: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
679: * function can only be called from inside TSSetFromOptions() */
680: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
681: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
682: if (flg || !((PetscObject)adapt)->type_name) {
683: TSAdaptSetType(adapt,type);
684: }
686: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
687: if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}
689: safety = adapt->safety; reject_safety = adapt->reject_safety;
690: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
691: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
692: if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}
694: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
695: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
696: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
697: if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}
699: hmin = adapt->dt_min; hmax = adapt->dt_max;
700: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
701: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
702: if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}
704: PetscOptionsReal("-ts_adapt_max_ignore","Adaptor ignores (absolute) solution values smaller than this value","",adapt->ignore_max,&adapt->ignore_max,&set);
705: PetscOptionsBool("-ts_adapt_glee_use_local","GLEE adaptor uses local error estimation for step control","",adapt->glee_use_local,&adapt->glee_use_local,&set);
707: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);
709: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
710: if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
712: 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);
714: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
715: if (set) {TSAdaptSetMonitor(adapt,flg);}
717: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
718: PetscOptionsTail();
719: return(0);
720: }
722: /*@
723: TSAdaptCandidatesClear - clear any previously set candidate schemes
725: Logically collective on TSAdapt
727: Input Argument:
728: . adapt - adaptive controller
730: Level: developer
732: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
733: @*/
734: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
735: {
740: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
741: return(0);
742: }
744: /*@C
745: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
747: Logically collective on TSAdapt
749: Input Arguments:
750: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
751: . name - name of the candidate scheme to add
752: . order - order of the candidate scheme
753: . stageorder - stage order of the candidate scheme
754: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
755: . cost - relative measure of the amount of work required for the candidate scheme
756: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
758: Note:
759: This routine is not available in Fortran.
761: Level: developer
763: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
764: @*/
765: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
766: {
767: PetscInt c;
771: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
772: if (inuse) {
773: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
774: adapt->candidates.inuse_set = PETSC_TRUE;
775: }
776: /* first slot if this is the current scheme, otherwise the next available slot */
777: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
779: adapt->candidates.name[c] = name;
780: adapt->candidates.order[c] = order;
781: adapt->candidates.stageorder[c] = stageorder;
782: adapt->candidates.ccfl[c] = ccfl;
783: adapt->candidates.cost[c] = cost;
784: adapt->candidates.n++;
785: return(0);
786: }
788: /*@C
789: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
791: Not Collective
793: Input Arguments:
794: . adapt - time step adaptivity context
796: Output Arguments:
797: + n - number of candidate schemes, always at least 1
798: . order - the order of each candidate scheme
799: . stageorder - the stage order of each candidate scheme
800: . ccfl - the CFL coefficient of each scheme
801: - cost - the relative cost of each scheme
803: Level: developer
805: Note:
806: The current scheme is always returned in the first slot
808: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
809: @*/
810: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
811: {
814: if (n) *n = adapt->candidates.n;
815: if (order) *order = adapt->candidates.order;
816: if (stageorder) *stageorder = adapt->candidates.stageorder;
817: if (ccfl) *ccfl = adapt->candidates.ccfl;
818: if (cost) *cost = adapt->candidates.cost;
819: return(0);
820: }
822: /*@C
823: TSAdaptChoose - choose which method and step size to use for the next step
825: Collective on TSAdapt
827: Input Arguments:
828: + adapt - adaptive contoller
829: - h - current step size
831: Output Arguments:
832: + next_sc - optional, scheme to use for the next step
833: . next_h - step size to use for the next step
834: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
836: Note:
837: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
838: being retried after an initial rejection.
840: Level: developer
842: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
843: @*/
844: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
845: {
847: PetscInt ncandidates = adapt->candidates.n;
848: PetscInt scheme = 0;
849: PetscReal wlte = -1.0;
850: PetscReal wltea = -1.0;
851: PetscReal wlter = -1.0;
859: if (next_sc) *next_sc = 0;
861: /* Do not mess with adaptivity while handling events*/
862: if (ts->event && ts->event->status != TSEVENT_NONE) {
863: *next_h = h;
864: *accept = PETSC_TRUE;
865: return(0);
866: }
868: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
869: 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);
870: if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
871: if (next_sc) *next_sc = scheme;
873: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
874: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
875: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
876: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
877: PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
878: PetscReal b = adapt->matchstepfac[1];
879: if (t < tmax && tend > tmax) *next_h = hmax;
880: if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
881: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
882: }
884: if (adapt->monitor) {
885: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
886: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
887: if (wlte < 0) {
888: 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);
889: } else {
890: 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);
891: }
892: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
893: }
894: return(0);
895: }
897: /*@
898: TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
899: before increasing the time step.
901: Logicially Collective on TSAdapt
903: Input Arguments:
904: + adapt - adaptive controller context
905: - cnt - the number of timesteps
907: Options Database Key:
908: . -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase
910: Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
911: The successful use of this option is problem dependent
913: Developer Note: there is no theory to support this option
915: Level: advanced
917: .seealso:
918: @*/
919: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
920: {
922: adapt->timestepjustdecreased_delay = cnt;
923: return(0);
924: }
926: /*@
927: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible)
929: Collective on TSAdapt
931: Input Arguments:
932: + adapt - adaptive controller context
933: . ts - time stepper
934: . t - Current simulation time
935: - Y - Current solution vector
937: Output Arguments:
938: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
940: Level: developer
942: .seealso:
943: @*/
944: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
945: {
946: PetscErrorCode ierr;
947: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
954: if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
955: if (snesreason < 0) {
956: *accept = PETSC_FALSE;
957: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
958: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
959: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
960: if (adapt->monitor) {
961: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
962: 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);
963: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
964: }
965: }
966: } else {
967: *accept = PETSC_TRUE;
968: TSFunctionDomainError(ts,t,Y,accept);
969: if (*accept && adapt->checkstage) {
970: (*adapt->checkstage)(adapt,ts,t,Y,accept);
971: if (!*accept) {
972: PetscInfo1(ts,"Step=%D, solution rejected by user function provided by TSSetFunctionDomainError()\n",ts->steps);
973: if (adapt->monitor) {
974: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
975: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected by user function provided by TSSetFunctionDomainError()\n",((PetscObject)adapt)->type_name,ts->steps);
976: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
977: }
978: }
979: }
980: }
982: if (!(*accept) && !ts->reason) {
983: PetscReal dt,new_dt;
984: TSGetTimeStep(ts,&dt);
985: new_dt = dt * adapt->scale_solve_failed;
986: TSSetTimeStep(ts,new_dt);
987: adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
988: if (adapt->monitor) {
989: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
990: 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);
991: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
992: }
993: }
994: return(0);
995: }
997: /*@
998: TSAdaptCreate - create an adaptive controller context for time stepping
1000: Collective
1002: Input Parameter:
1003: . comm - The communicator
1005: Output Parameter:
1006: . adapt - new TSAdapt object
1008: Level: developer
1010: Notes:
1011: TSAdapt creation is handled by TS, so users should not need to call this function.
1013: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1014: @*/
1015: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1016: {
1018: TSAdapt adapt;
1022: *inadapt = NULL;
1023: TSAdaptInitializePackage();
1025: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
1027: adapt->always_accept = PETSC_FALSE;
1028: adapt->safety = 0.9;
1029: adapt->reject_safety = 0.5;
1030: adapt->clip[0] = 0.1;
1031: adapt->clip[1] = 10.;
1032: adapt->dt_min = 1e-20;
1033: adapt->dt_max = 1e+20;
1034: adapt->ignore_max = -1.0;
1035: adapt->glee_use_local = PETSC_TRUE;
1036: adapt->scale_solve_failed = 0.25;
1037: /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1038: to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1039: adapt->matchstepfac[0] = 0.01; /* allow 1% step size increase in the last step */
1040: adapt->matchstepfac[1] = 2.0; /* halve last step if it is greater than what remains divided this factor */
1041: adapt->wnormtype = NORM_2;
1042: adapt->timestepjustdecreased_delay = 0;
1044: *inadapt = adapt;
1045: return(0);
1046: }