Actual source code: tsadapt.c
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 Parameters:
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 Parameters:
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 Parameters:
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: {
378: adapt->checkstage = func;
379: return(0);
380: }
382: /*@
383: TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of
384: any error or stability condition not meeting the prescribed goal.
386: Logically collective on TSAdapt
388: Input Parameters:
389: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
390: - flag - whether to always accept steps
392: Options Database Keys:
393: . -ts_adapt_always_accept - to always accept steps
395: Level: intermediate
397: .seealso: TSAdapt, TSAdaptChoose()
398: @*/
399: PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)
400: {
404: adapt->always_accept = flag;
405: return(0);
406: }
408: /*@
409: TSAdaptSetSafety - Set safety factors
411: Logically collective on TSAdapt
413: Input Parameters:
414: + adapt - adaptive controller context
415: . safety - safety factor relative to target error/stability goal
416: - reject_safety - extra safety factor to apply if the last step was rejected
418: Options Database Keys:
419: + -ts_adapt_safety <safety> - to set safety factor
420: - -ts_adapt_reject_safety <reject_safety> - to set reject safety factor
422: Level: intermediate
424: .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose()
425: @*/
426: PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)
427: {
432: if (safety != PETSC_DEFAULT && safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be non negative",(double)safety);
433: if (safety != PETSC_DEFAULT && safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be less than one",(double)safety);
434: 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);
435: 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);
436: if (safety != PETSC_DEFAULT) adapt->safety = safety;
437: if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety;
438: return(0);
439: }
441: /*@
442: TSAdaptGetSafety - Get safety factors
444: Not Collective
446: Input Parameter:
447: . adapt - adaptive controller context
449: Output Parameters:
450: . safety - safety factor relative to target error/stability goal
451: + reject_safety - extra safety factor to apply if the last step was rejected
453: Level: intermediate
455: .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose()
456: @*/
457: PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety)
458: {
463: if (safety) *safety = adapt->safety;
464: if (reject_safety) *reject_safety = adapt->reject_safety;
465: return(0);
466: }
468: /*@
469: 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.
471: Logically collective on TSAdapt
473: Input Parameters:
474: + adapt - adaptive controller context
475: - max_ignore - threshold for solution components that are ignored during error estimation
477: Options Database Keys:
478: . -ts_adapt_max_ignore <max_ignore> - to set the threshold
480: Level: intermediate
482: .seealso: TSAdapt, TSAdaptGetMaxIgnore(), TSAdaptChoose()
483: @*/
484: PetscErrorCode TSAdaptSetMaxIgnore(TSAdapt adapt,PetscReal max_ignore)
485: {
489: adapt->ignore_max = max_ignore;
490: return(0);
491: }
493: /*@
494: 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).
496: Not Collective
498: Input Parameter:
499: . adapt - adaptive controller context
501: Output Parameter:
502: . max_ignore - threshold for solution components that are ignored during error estimation
504: Level: intermediate
506: .seealso: TSAdapt, TSAdaptSetMaxIgnore(), TSAdaptChoose()
507: @*/
508: PetscErrorCode TSAdaptGetMaxIgnore(TSAdapt adapt,PetscReal *max_ignore)
509: {
513: *max_ignore = adapt->ignore_max;
514: return(0);
515: }
517: /*@
518: TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
520: Logically collective on TSAdapt
522: Input Parameters:
523: + adapt - adaptive controller context
524: . low - admissible decrease factor
525: - high - admissible increase factor
527: Options Database Keys:
528: . -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors
530: Level: intermediate
532: .seealso: TSAdaptChoose(), TSAdaptGetClip(), TSAdaptSetScaleSolveFailed()
533: @*/
534: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
535: {
540: if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
541: if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
542: if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be greater than one",(double)high);
543: if (low != PETSC_DEFAULT) adapt->clip[0] = low;
544: if (high != PETSC_DEFAULT) adapt->clip[1] = high;
545: return(0);
546: }
548: /*@
549: TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
551: Not Collective
553: Input Parameter:
554: . adapt - adaptive controller context
556: Output Parameters:
557: + low - optional, admissible decrease factor
558: - high - optional, admissible increase factor
560: Level: intermediate
562: .seealso: TSAdaptChoose(), TSAdaptSetClip(), TSAdaptSetScaleSolveFailed()
563: @*/
564: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
565: {
570: if (low) *low = adapt->clip[0];
571: if (high) *high = adapt->clip[1];
572: return(0);
573: }
575: /*@
576: TSAdaptSetScaleSolveFailed - Scale step by this factor if solve fails
578: Logically collective on TSAdapt
580: Input Parameters:
581: + adapt - adaptive controller context
582: - scale - scale
584: Options Database Keys:
585: . -ts_adapt_scale_solve_failed <scale> - to set scale step by this factor if solve fails
587: Level: intermediate
589: .seealso: TSAdaptChoose(), TSAdaptGetScaleSolveFailed(), TSAdaptGetClip()
590: @*/
591: PetscErrorCode TSAdaptSetScaleSolveFailed(TSAdapt adapt,PetscReal scale)
592: {
596: if (scale != PETSC_DEFAULT && scale <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be positive",(double)scale);
597: if (scale != PETSC_DEFAULT && scale > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be less than one",(double)scale);
598: if (scale != PETSC_DEFAULT) adapt->scale_solve_failed = scale;
599: return(0);
600: }
602: /*@
603: TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size
605: Not Collective
607: Input Parameter:
608: . adapt - adaptive controller context
610: Output Parameter:
611: . scale - scale factor
613: Level: intermediate
615: .seealso: TSAdaptChoose(), TSAdaptSetScaleSolveFailed(), TSAdaptSetClip()
616: @*/
617: PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal *scale)
618: {
622: if (scale) *scale = adapt->scale_solve_failed;
623: return(0);
624: }
626: /*@
627: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
629: Logically collective on TSAdapt
631: Input Parameters:
632: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
633: . hmin - minimum time step
634: - hmax - maximum time step
636: Options Database Keys:
637: + -ts_adapt_dt_min <min> - to set minimum time step
638: - -ts_adapt_dt_max <max> - to set maximum time step
640: Level: intermediate
642: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
643: @*/
644: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
645: {
651: if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
652: if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
653: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
654: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
655: hmin = adapt->dt_min;
656: hmax = adapt->dt_max;
657: 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);
658: return(0);
659: }
661: /*@
662: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
664: Not Collective
666: Input Parameter:
667: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
669: Output Parameters:
670: + hmin - minimum time step
671: - hmax - maximum time step
673: Level: intermediate
675: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
676: @*/
677: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
678: {
684: if (hmin) *hmin = adapt->dt_min;
685: if (hmax) *hmax = adapt->dt_max;
686: return(0);
687: }
689: /*
690: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
692: Collective on TSAdapt
694: Input Parameter:
695: . adapt - the TSAdapt context
697: Options Database Keys:
698: + -ts_adapt_type <type> - algorithm to use for adaptivity
699: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
700: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
701: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
702: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
703: . -ts_adapt_dt_min <min> - minimum timestep to use
704: . -ts_adapt_dt_max <max> - maximum timestep to use
705: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
706: . -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
707: - -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver
709: Level: advanced
711: Notes:
712: This function is automatically called by TSSetFromOptions()
714: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
715: TSAdaptSetClip(), TSAdaptSetScaleSolveFailed(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
716: */
717: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
718: {
720: char type[256] = TSADAPTBASIC;
721: PetscReal safety,reject_safety,clip[2],scale,hmin,hmax;
722: PetscBool set,flg;
723: PetscInt two;
727: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
728: * function can only be called from inside TSSetFromOptions() */
729: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
730: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
731: if (flg || !((PetscObject)adapt)->type_name) {
732: TSAdaptSetType(adapt,type);
733: }
735: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
736: if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}
738: safety = adapt->safety; reject_safety = adapt->reject_safety;
739: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
740: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
741: if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}
743: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
744: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
745: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
746: if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}
748: hmin = adapt->dt_min; hmax = adapt->dt_max;
749: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
750: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
751: if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}
753: PetscOptionsReal("-ts_adapt_max_ignore","Adaptor ignores (absolute) solution values smaller than this value","",adapt->ignore_max,&adapt->ignore_max,&set);
754: PetscOptionsBool("-ts_adapt_glee_use_local","GLEE adaptor uses local error estimation for step control","",adapt->glee_use_local,&adapt->glee_use_local,&set);
756: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","TSAdaptSetScaleSolveFailed",adapt->scale_solve_failed,&scale,&set);
757: if (set) {TSAdaptSetScaleSolveFailed(adapt,scale);}
759: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
760: if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
762: 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);
764: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
765: if (set) {TSAdaptSetMonitor(adapt,flg);}
767: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
768: PetscOptionsTail();
769: return(0);
770: }
772: /*@
773: TSAdaptCandidatesClear - clear any previously set candidate schemes
775: Logically collective on TSAdapt
777: Input Parameter:
778: . adapt - adaptive controller
780: Level: developer
782: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
783: @*/
784: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
785: {
790: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
791: return(0);
792: }
794: /*@C
795: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
797: Logically collective on TSAdapt
799: Input Parameters:
800: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
801: . name - name of the candidate scheme to add
802: . order - order of the candidate scheme
803: . stageorder - stage order of the candidate scheme
804: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
805: . cost - relative measure of the amount of work required for the candidate scheme
806: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
808: Note:
809: This routine is not available in Fortran.
811: Level: developer
813: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
814: @*/
815: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
816: {
817: PetscInt c;
821: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
822: if (inuse) {
823: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
824: adapt->candidates.inuse_set = PETSC_TRUE;
825: }
826: /* first slot if this is the current scheme, otherwise the next available slot */
827: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
829: adapt->candidates.name[c] = name;
830: adapt->candidates.order[c] = order;
831: adapt->candidates.stageorder[c] = stageorder;
832: adapt->candidates.ccfl[c] = ccfl;
833: adapt->candidates.cost[c] = cost;
834: adapt->candidates.n++;
835: return(0);
836: }
838: /*@C
839: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
841: Not Collective
843: Input Parameter:
844: . adapt - time step adaptivity context
846: Output Parameters:
847: + n - number of candidate schemes, always at least 1
848: . order - the order of each candidate scheme
849: . stageorder - the stage order of each candidate scheme
850: . ccfl - the CFL coefficient of each scheme
851: - cost - the relative cost of each scheme
853: Level: developer
855: Note:
856: The current scheme is always returned in the first slot
858: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
859: @*/
860: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
861: {
864: if (n) *n = adapt->candidates.n;
865: if (order) *order = adapt->candidates.order;
866: if (stageorder) *stageorder = adapt->candidates.stageorder;
867: if (ccfl) *ccfl = adapt->candidates.ccfl;
868: if (cost) *cost = adapt->candidates.cost;
869: return(0);
870: }
872: /*@C
873: TSAdaptChoose - choose which method and step size to use for the next step
875: Collective on TSAdapt
877: Input Parameters:
878: + adapt - adaptive contoller
879: . ts - time stepper
880: - h - current step size
882: Output Parameters:
883: + next_sc - optional, scheme to use for the next step
884: . next_h - step size to use for the next step
885: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
887: Note:
888: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
889: being retried after an initial rejection.
891: Level: developer
893: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
894: @*/
895: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
896: {
898: PetscInt ncandidates = adapt->candidates.n;
899: PetscInt scheme = 0;
900: PetscReal wlte = -1.0;
901: PetscReal wltea = -1.0;
902: PetscReal wlter = -1.0;
910: if (next_sc) *next_sc = 0;
912: /* Do not mess with adaptivity while handling events*/
913: if (ts->event && ts->event->status != TSEVENT_NONE) {
914: *next_h = h;
915: *accept = PETSC_TRUE;
916: return(0);
917: }
919: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
920: 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);
921: if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
922: if (next_sc) *next_sc = scheme;
924: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
925: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
926: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
927: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
928: PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
929: PetscReal b = adapt->matchstepfac[1];
930: if (t < tmax && tend > tmax) *next_h = hmax;
931: if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
932: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
933: }
935: if (adapt->monitor) {
936: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
937: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
938: if (wlte < 0) {
939: 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);
940: } else {
941: 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);
942: }
943: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
944: }
945: return(0);
946: }
948: /*@
949: TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
950: before increasing the time step.
952: Logicially Collective on TSAdapt
954: Input Parameters:
955: + adapt - adaptive controller context
956: - cnt - the number of timesteps
958: Options Database Key:
959: . -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase
961: Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
962: The successful use of this option is problem dependent
964: Developer Note: there is no theory to support this option
966: Level: advanced
968: .seealso:
969: @*/
970: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
971: {
973: adapt->timestepjustdecreased_delay = cnt;
974: return(0);
975: }
977: /*@
978: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible)
980: Collective on TSAdapt
982: Input Parameters:
983: + adapt - adaptive controller context
984: . ts - time stepper
985: . t - Current simulation time
986: - Y - Current solution vector
988: Output Parameter:
989: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
991: Level: developer
993: .seealso:
994: @*/
995: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
996: {
997: PetscErrorCode ierr;
998: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
1005: if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
1006: if (snesreason < 0) {
1007: *accept = PETSC_FALSE;
1008: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
1009: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
1010: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
1011: if (adapt->monitor) {
1012: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1013: 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);
1014: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1015: }
1016: }
1017: } else {
1018: *accept = PETSC_TRUE;
1019: TSFunctionDomainError(ts,t,Y,accept);
1020: if (*accept && adapt->checkstage) {
1021: (*adapt->checkstage)(adapt,ts,t,Y,accept);
1022: if (!*accept) {
1023: PetscInfo1(ts,"Step=%D, solution rejected by user function provided by TSSetFunctionDomainError()\n",ts->steps);
1024: if (adapt->monitor) {
1025: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1026: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected by user function provided by TSSetFunctionDomainError()\n",((PetscObject)adapt)->type_name,ts->steps);
1027: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1028: }
1029: }
1030: }
1031: }
1033: if (!(*accept) && !ts->reason) {
1034: PetscReal dt,new_dt;
1035: TSGetTimeStep(ts,&dt);
1036: new_dt = dt * adapt->scale_solve_failed;
1037: TSSetTimeStep(ts,new_dt);
1038: adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
1039: if (adapt->monitor) {
1040: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1041: 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);
1042: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1043: }
1044: }
1045: return(0);
1046: }
1048: /*@
1049: TSAdaptCreate - create an adaptive controller context for time stepping
1051: Collective
1053: Input Parameter:
1054: . comm - The communicator
1056: Output Parameter:
1057: . adapt - new TSAdapt object
1059: Level: developer
1061: Notes:
1062: TSAdapt creation is handled by TS, so users should not need to call this function.
1064: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1065: @*/
1066: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1067: {
1069: TSAdapt adapt;
1073: *inadapt = NULL;
1074: TSAdaptInitializePackage();
1076: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
1078: adapt->always_accept = PETSC_FALSE;
1079: adapt->safety = 0.9;
1080: adapt->reject_safety = 0.5;
1081: adapt->clip[0] = 0.1;
1082: adapt->clip[1] = 10.;
1083: adapt->dt_min = 1e-20;
1084: adapt->dt_max = 1e+20;
1085: adapt->ignore_max = -1.0;
1086: adapt->glee_use_local = PETSC_TRUE;
1087: adapt->scale_solve_failed = 0.25;
1088: /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1089: to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1090: adapt->matchstepfac[0] = 0.01; /* allow 1% step size increase in the last step */
1091: adapt->matchstepfac[1] = 2.0; /* halve last step if it is greater than what remains divided this factor */
1092: adapt->wnormtype = NORM_2;
1093: adapt->timestepjustdecreased_delay = 0;
1095: *inadapt = adapt;
1096: return(0);
1097: }