Actual source code: tsadapt.c
petsc-3.10.5 2019-03-28
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);
16: /*@C
17: TSAdaptRegister - adds a TSAdapt implementation
19: Not Collective
21: Input Parameters:
22: + name_scheme - name of user-defined adaptivity scheme
23: - routine_create - routine to create method context
25: Notes:
26: TSAdaptRegister() may be called multiple times to add several user-defined families.
28: Sample usage:
29: .vb
30: TSAdaptRegister("my_scheme",MySchemeCreate);
31: .ve
33: Then, your scheme can be chosen with the procedural interface via
34: $ TSAdaptSetType(ts,"my_scheme")
35: or at runtime via the option
36: $ -ts_adapt_type my_scheme
38: Level: advanced
40: .keywords: TSAdapt, register
42: .seealso: TSAdaptRegisterAll()
43: @*/
44: PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
45: {
49: TSAdaptInitializePackage();
50: PetscFunctionListAdd(&TSAdaptList,sname,function);
51: return(0);
52: }
54: /*@C
55: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
57: Not Collective
59: Level: advanced
61: .keywords: TSAdapt, register, all
63: .seealso: TSAdaptRegisterDestroy()
64: @*/
65: PetscErrorCode TSAdaptRegisterAll(void)
66: {
70: if (TSAdaptRegisterAllCalled) return(0);
71: TSAdaptRegisterAllCalled = PETSC_TRUE;
72: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
73: TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);
74: TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP);
75: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
76: TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);
77: return(0);
78: }
80: /*@C
81: TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
82: called from PetscFinalize().
84: Level: developer
86: .keywords: Petsc, destroy, package
87: .seealso: PetscFinalize()
88: @*/
89: PetscErrorCode TSAdaptFinalizePackage(void)
90: {
94: PetscFunctionListDestroy(&TSAdaptList);
95: TSAdaptPackageInitialized = PETSC_FALSE;
96: TSAdaptRegisterAllCalled = PETSC_FALSE;
97: return(0);
98: }
100: /*@C
101: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
102: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
103: TSAdaptCreate() when using static libraries.
105: Level: developer
107: .keywords: TSAdapt, initialize, package
108: .seealso: PetscInitialize()
109: @*/
110: PetscErrorCode TSAdaptInitializePackage(void)
111: {
115: if (TSAdaptPackageInitialized) return(0);
116: TSAdaptPackageInitialized = PETSC_TRUE;
117: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
118: TSAdaptRegisterAll();
119: PetscRegisterFinalize(TSAdaptFinalizePackage);
120: return(0);
121: }
123: /*@C
124: TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE
126: Logicially Collective on TSAdapt
128: Input Parameter:
129: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
130: - type - either TSADAPTBASIC or TSADAPTNONE
132: Options Database:
133: . -ts_adapt_type <basic or dsp or none> - to set the adapter type
135: Level: intermediate
137: .keywords: TSAdapt, create
139: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
140: @*/
141: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
142: {
143: PetscBool match;
144: PetscErrorCode ierr,(*r)(TSAdapt);
149: PetscObjectTypeCompare((PetscObject)adapt,type,&match);
150: if (match) return(0);
151: PetscFunctionListFind(TSAdaptList,type,&r);
152: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
153: if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
154: PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
155: PetscObjectChangeTypeName((PetscObject)adapt,type);
156: (*r)(adapt);
157: return(0);
158: }
160: /*@C
161: TSAdaptGetType - gets the TS adapter method type (as a string).
163: Not Collective
165: Input Parameter:
166: . adapt - The TS adapter, most likely obtained with TSGetAdapt()
168: Output Parameter:
169: . type - The name of TS adapter method
171: Level: intermediate
173: .keywords: TSAdapt, get, type
174: .seealso TSAdaptSetType()
175: @*/
176: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
177: {
181: *type = ((PetscObject)adapt)->type_name;
182: return(0);
183: }
185: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
186: {
191: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
192: return(0);
193: }
195: /*@C
196: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
198: Collective on PetscViewer
200: Input Parameters:
201: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
202: some related function before a call to TSAdaptLoad().
203: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
204: HDF5 file viewer, obtained from PetscViewerHDF5Open()
206: Level: intermediate
208: Notes:
209: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
211: Notes for advanced users:
212: Most users should not need to know the details of the binary storage
213: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
214: But for anyone who's interested, the standard binary matrix storage
215: format is
216: .vb
217: has not yet been determined
218: .ve
220: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
221: @*/
222: PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
223: {
225: PetscBool isbinary;
226: char type[256];
231: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
232: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
234: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
235: TSAdaptSetType(adapt,type);
236: if (adapt->ops->load) {
237: (*adapt->ops->load)(adapt,viewer);
238: }
239: return(0);
240: }
242: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
243: {
245: PetscBool iascii,isbinary,isnone;
249: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
252: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
253: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
254: if (iascii) {
255: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
256: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
257: if (!isnone) {
258: if (adapt->always_accept) {PetscViewerASCIIPrintf(viewer," always accepting steps\n");}
259: PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);
260: PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
261: PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);
262: PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);
263: PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);
264: PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);
265: }
266: if (adapt->ops->view) {
267: PetscViewerASCIIPushTab(viewer);
268: (*adapt->ops->view)(adapt,viewer);
269: PetscViewerASCIIPopTab(viewer);
270: }
271: } else if (isbinary) {
272: char type[256];
274: /* need to save FILE_CLASS_ID for adapt class */
275: PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
276: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
277: } else if (adapt->ops->view) {
278: (*adapt->ops->view)(adapt,viewer);
279: }
280: return(0);
281: }
283: /*@
284: TSAdaptReset - Resets a TSAdapt context.
286: Collective on TS
288: Input Parameter:
289: . adapt - the TSAdapt context obtained from TSAdaptCreate()
291: Level: developer
293: .seealso: TSAdaptCreate(), TSAdaptDestroy()
294: @*/
295: PetscErrorCode TSAdaptReset(TSAdapt adapt)
296: {
301: if (adapt->ops->reset) {(*adapt->ops->reset)(adapt);}
302: return(0);
303: }
305: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
306: {
310: if (!*adapt) return(0);
312: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return(0);}
314: TSAdaptReset(*adapt);
316: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
317: PetscViewerDestroy(&(*adapt)->monitor);
318: PetscHeaderDestroy(adapt);
319: return(0);
320: }
322: /*@
323: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
325: Collective on TSAdapt
327: Input Arguments:
328: + adapt - adaptive controller context
329: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
331: Options Database Keys:
332: . -ts_adapt_monitor - to turn on monitoring
334: Level: intermediate
336: .seealso: TSAdaptChoose()
337: @*/
338: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
339: {
345: if (flg) {
346: if (!adapt->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);}
347: } else {
348: PetscViewerDestroy(&adapt->monitor);
349: }
350: return(0);
351: }
353: /*@C
354: TSAdaptSetCheckStage - Set a callback to check convergence for a stage
356: Logically collective on TSAdapt
358: Input Arguments:
359: + adapt - adaptive controller context
360: - func - stage check function
362: Arguments of func:
363: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
365: + adapt - adaptive controller context
366: . ts - time stepping context
367: - accept - pending choice of whether to accept, can be modified by this routine
369: Level: advanced
371: .seealso: TSAdaptChoose()
372: @*/
373: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*))
374: {
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 Arguments:
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 Arguments:
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 Arguments:
447: . adapt - adaptive controller context
449: Ouput Arguments:
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: TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
471: Logically collective on TSAdapt
473: Input Arguments:
474: + adapt - adaptive controller context
475: . low - admissible decrease factor
476: - high - admissible increase factor
478: Options Database Keys:
479: . -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors
481: Level: intermediate
483: .seealso: TSAdaptChoose(), TSAdaptGetClip()
484: @*/
485: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
486: {
491: if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
492: if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
493: if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high);
494: if (low != PETSC_DEFAULT) adapt->clip[0] = low;
495: if (high != PETSC_DEFAULT) adapt->clip[1] = high;
496: return(0);
497: }
499: /*@
500: TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
502: Not Collective
504: Input Arguments:
505: . adapt - adaptive controller context
507: Ouput Arguments:
508: + low - optional, admissible decrease factor
509: - high - optional, admissible increase factor
511: Level: intermediate
513: .seealso: TSAdaptChoose(), TSAdaptSetClip()
514: @*/
515: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
516: {
521: if (low) *low = adapt->clip[0];
522: if (high) *high = adapt->clip[1];
523: return(0);
524: }
526: /*@
527: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
529: Logically collective on TSAdapt
531: Input Arguments:
532: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
533: . hmin - minimum time step
534: - hmax - maximum time step
536: Options Database Keys:
537: + -ts_adapt_dt_min <min> - to set minimum time step
538: - -ts_adapt_dt_max <max> - to set maximum time step
540: Level: intermediate
542: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
543: @*/
544: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
545: {
551: if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
552: if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
553: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
554: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
555: hmin = adapt->dt_min;
556: hmax = adapt->dt_max;
557: 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);
558: return(0);
559: }
561: /*@
562: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
564: Not Collective
566: Input Arguments:
567: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
569: Output Arguments:
570: + hmin - minimum time step
571: - hmax - maximum time step
573: Level: intermediate
575: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
576: @*/
577: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
578: {
584: if (hmin) *hmin = adapt->dt_min;
585: if (hmax) *hmax = adapt->dt_max;
586: return(0);
587: }
589: /*
590: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
592: Collective on TSAdapt
594: Input Parameter:
595: . adapt - the TSAdapt context
597: Options Database Keys:
598: + -ts_adapt_type <type> - algorithm to use for adaptivity
599: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
600: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
601: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
602: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
603: . -ts_adapt_dt_min <min> - minimum timestep to use
604: . -ts_adapt_dt_max <max> - maximum timestep to use
605: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
606: - -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
608: Level: advanced
610: Notes:
611: This function is automatically called by TSSetFromOptions()
613: .keywords: TS, TSGetAdapt(), TSAdaptSetType(), TSAdaptSetStepLimits()
615: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
616: TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
617: */
618: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
619: {
621: char type[256] = TSADAPTBASIC;
622: PetscReal safety,reject_safety,clip[2],hmin,hmax;
623: PetscBool set,flg;
624: PetscInt two;
628: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
629: * function can only be called from inside TSSetFromOptions() */
630: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
631: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
632: if (flg || !((PetscObject)adapt)->type_name) {
633: TSAdaptSetType(adapt,type);
634: }
636: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
637: if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}
639: safety = adapt->safety; reject_safety = adapt->reject_safety;
640: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
641: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
642: if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}
644: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
645: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
646: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
647: if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}
649: hmin = adapt->dt_min; hmax = adapt->dt_max;
650: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
651: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
652: if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}
654: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);
656: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
657: if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
659: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
660: if (set) {TSAdaptSetMonitor(adapt,flg);}
662: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
663: PetscOptionsTail();
664: return(0);
665: }
667: /*@
668: TSAdaptCandidatesClear - clear any previously set candidate schemes
670: Logically collective on TSAdapt
672: Input Argument:
673: . adapt - adaptive controller
675: Level: developer
677: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
678: @*/
679: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
680: {
685: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
686: return(0);
687: }
689: /*@C
690: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
692: Logically collective on TSAdapt
694: Input Arguments:
695: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
696: . name - name of the candidate scheme to add
697: . order - order of the candidate scheme
698: . stageorder - stage order of the candidate scheme
699: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
700: . cost - relative measure of the amount of work required for the candidate scheme
701: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
703: Note:
704: This routine is not available in Fortran.
706: Level: developer
708: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
709: @*/
710: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
711: {
712: PetscInt c;
716: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
717: if (inuse) {
718: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
719: adapt->candidates.inuse_set = PETSC_TRUE;
720: }
721: /* first slot if this is the current scheme, otherwise the next available slot */
722: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
724: adapt->candidates.name[c] = name;
725: adapt->candidates.order[c] = order;
726: adapt->candidates.stageorder[c] = stageorder;
727: adapt->candidates.ccfl[c] = ccfl;
728: adapt->candidates.cost[c] = cost;
729: adapt->candidates.n++;
730: return(0);
731: }
733: /*@C
734: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
736: Not Collective
738: Input Arguments:
739: . adapt - time step adaptivity context
741: Output Arguments:
742: + n - number of candidate schemes, always at least 1
743: . order - the order of each candidate scheme
744: . stageorder - the stage order of each candidate scheme
745: . ccfl - the CFL coefficient of each scheme
746: - cost - the relative cost of each scheme
748: Level: developer
750: Note:
751: The current scheme is always returned in the first slot
753: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
754: @*/
755: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
756: {
759: if (n) *n = adapt->candidates.n;
760: if (order) *order = adapt->candidates.order;
761: if (stageorder) *stageorder = adapt->candidates.stageorder;
762: if (ccfl) *ccfl = adapt->candidates.ccfl;
763: if (cost) *cost = adapt->candidates.cost;
764: return(0);
765: }
767: /*@C
768: TSAdaptChoose - choose which method and step size to use for the next step
770: Collective on TSAdapt
772: Input Arguments:
773: + adapt - adaptive contoller
774: - h - current step size
776: Output Arguments:
777: + next_sc - optional, scheme to use for the next step
778: . next_h - step size to use for the next step
779: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
781: Note:
782: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
783: being retried after an initial rejection.
785: Level: developer
787: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
788: @*/
789: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
790: {
792: PetscInt ncandidates = adapt->candidates.n;
793: PetscInt scheme = 0;
794: PetscReal wlte = -1.0;
795: PetscReal wltea = -1.0;
796: PetscReal wlter = -1.0;
804: if (next_sc) *next_sc = 0;
806: /* Do not mess with adaptivity while handling events*/
807: if (ts->event && ts->event->status != TSEVENT_NONE) {
808: *next_h = h;
809: *accept = PETSC_TRUE;
810: return(0);
811: }
813: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
814: 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);
815: if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
816: if (next_sc) *next_sc = scheme;
818: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
819: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
820: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
821: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
822: PetscReal a = (PetscReal)1.01; /* allow 1% step size increase */
823: if (t < tmax && tend > tmax) *next_h = hmax;
824: if (t < tmax && tend < tmax && h > hmax/2) *next_h = hmax/2;
825: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
826: }
828: if (adapt->monitor) {
829: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
830: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
831: if (wlte < 0) {
832: 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);
833: } else {
834: 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);
835: }
836: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
837: }
838: return(0);
839: }
841: /*@
842: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)
844: Collective on TSAdapt
846: Input Arguments:
847: + adapt - adaptive controller context
848: . ts - time stepper
849: . t - Current simulation time
850: - Y - Current solution vector
852: Output Arguments:
853: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
855: Level: developer
857: .seealso:
858: @*/
859: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
860: {
861: PetscErrorCode ierr;
862: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
869: if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
870: if (snesreason < 0) {
871: *accept = PETSC_FALSE;
872: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
873: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
874: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
875: if (adapt->monitor) {
876: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
877: 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);
878: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
879: }
880: }
881: } else {
882: *accept = PETSC_TRUE;
883: TSFunctionDomainError(ts,t,Y,accept);
884: if(*accept && adapt->checkstage) {
885: (*adapt->checkstage)(adapt,ts,t,Y,accept);
886: }
887: }
889: if(!(*accept) && !ts->reason) {
890: PetscReal dt,new_dt;
891: TSGetTimeStep(ts,&dt);
892: new_dt = dt * adapt->scale_solve_failed;
893: TSSetTimeStep(ts,new_dt);
894: adapt->timestepjustincreased += 4;
895: if (adapt->monitor) {
896: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
897: 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);
898: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
899: }
900: }
901: return(0);
902: }
904: /*@
905: TSAdaptCreate - create an adaptive controller context for time stepping
907: Collective on MPI_Comm
909: Input Parameter:
910: . comm - The communicator
912: Output Parameter:
913: . adapt - new TSAdapt object
915: Level: developer
917: Notes:
918: TSAdapt creation is handled by TS, so users should not need to call this function.
920: .keywords: TSAdapt, create
921: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
922: @*/
923: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
924: {
926: TSAdapt adapt;
930: *inadapt = NULL;
931: TSAdaptInitializePackage();
933: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
935: adapt->always_accept = PETSC_FALSE;
936: adapt->safety = 0.9;
937: adapt->reject_safety = 0.5;
938: adapt->clip[0] = 0.1;
939: adapt->clip[1] = 10.;
940: adapt->dt_min = 1e-20;
941: adapt->dt_max = 1e+20;
942: adapt->scale_solve_failed = 0.25;
943: adapt->wnormtype = NORM_2;
945: *inadapt = adapt;
946: return(0);
947: }