Actual source code: tsadapt.c
petsc-3.11.4 2019-09-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);
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: .keywords: TSAdapt, register
43: .seealso: TSAdaptRegisterAll()
44: @*/
45: PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
46: {
50: TSAdaptInitializePackage();
51: PetscFunctionListAdd(&TSAdaptList,sname,function);
52: return(0);
53: }
55: /*@C
56: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
58: Not Collective
60: Level: advanced
62: .keywords: TSAdapt, register, all
64: .seealso: TSAdaptRegisterDestroy()
65: @*/
66: PetscErrorCode TSAdaptRegisterAll(void)
67: {
71: if (TSAdaptRegisterAllCalled) return(0);
72: TSAdaptRegisterAllCalled = PETSC_TRUE;
73: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
74: TSAdaptRegister(TSADAPTBASIC, TSAdaptCreate_Basic);
75: TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP);
76: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
77: TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);
78: TSAdaptRegister(TSADAPTHISTORY,TSAdaptCreate_History);
79: return(0);
80: }
82: /*@C
83: TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
84: called from PetscFinalize().
86: Level: developer
88: .keywords: Petsc, destroy, package
89: .seealso: PetscFinalize()
90: @*/
91: PetscErrorCode TSAdaptFinalizePackage(void)
92: {
96: PetscFunctionListDestroy(&TSAdaptList);
97: TSAdaptPackageInitialized = PETSC_FALSE;
98: TSAdaptRegisterAllCalled = PETSC_FALSE;
99: return(0);
100: }
102: /*@C
103: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
104: called from TSInitializePackage().
106: Level: developer
108: .keywords: TSAdapt, initialize, package
109: .seealso: PetscInitialize()
110: @*/
111: PetscErrorCode TSAdaptInitializePackage(void)
112: {
116: if (TSAdaptPackageInitialized) return(0);
117: TSAdaptPackageInitialized = PETSC_TRUE;
118: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
119: TSAdaptRegisterAll();
120: PetscRegisterFinalize(TSAdaptFinalizePackage);
121: return(0);
122: }
124: /*@C
125: TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE
127: Logicially Collective on TSAdapt
129: Input Parameter:
130: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
131: - type - either TSADAPTBASIC or TSADAPTNONE
133: Options Database:
134: . -ts_adapt_type <basic or dsp or none> - to set the adapter type
136: Level: intermediate
138: .keywords: TSAdapt, create
140: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
141: @*/
142: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
143: {
144: PetscBool match;
145: PetscErrorCode ierr,(*r)(TSAdapt);
150: PetscObjectTypeCompare((PetscObject)adapt,type,&match);
151: if (match) return(0);
152: PetscFunctionListFind(TSAdaptList,type,&r);
153: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
154: if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
155: PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
156: PetscObjectChangeTypeName((PetscObject)adapt,type);
157: (*r)(adapt);
158: return(0);
159: }
161: /*@C
162: TSAdaptGetType - gets the TS adapter method type (as a string).
164: Not Collective
166: Input Parameter:
167: . adapt - The TS adapter, most likely obtained with TSGetAdapt()
169: Output Parameter:
170: . type - The name of TS adapter method
172: Level: intermediate
174: .keywords: TSAdapt, get, type
175: .seealso TSAdaptSetType()
176: @*/
177: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
178: {
182: *type = ((PetscObject)adapt)->type_name;
183: return(0);
184: }
186: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
187: {
192: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
193: return(0);
194: }
196: /*@C
197: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
199: Collective on PetscViewer
201: Input Parameters:
202: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
203: some related function before a call to TSAdaptLoad().
204: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
205: HDF5 file viewer, obtained from PetscViewerHDF5Open()
207: Level: intermediate
209: Notes:
210: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
212: Notes for advanced users:
213: Most users should not need to know the details of the binary storage
214: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
215: But for anyone who's interested, the standard binary matrix storage
216: format is
217: .vb
218: has not yet been determined
219: .ve
221: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
222: @*/
223: PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
224: {
226: PetscBool isbinary;
227: char type[256];
232: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
233: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
235: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
236: TSAdaptSetType(adapt,type);
237: if (adapt->ops->load) {
238: (*adapt->ops->load)(adapt,viewer);
239: }
240: return(0);
241: }
243: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
244: {
246: PetscBool iascii,isbinary,isnone;
250: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
253: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
254: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
255: if (iascii) {
256: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
257: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
258: if (!isnone) {
259: if (adapt->always_accept) {PetscViewerASCIIPrintf(viewer," always accepting steps\n");}
260: PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);
261: PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
262: PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);
263: PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);
264: PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);
265: PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);
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: TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
472: Logically collective on TSAdapt
474: Input Arguments:
475: + adapt - adaptive controller context
476: . low - admissible decrease factor
477: - high - admissible increase factor
479: Options Database Keys:
480: . -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors
482: Level: intermediate
484: .seealso: TSAdaptChoose(), TSAdaptGetClip()
485: @*/
486: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
487: {
492: if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
493: if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
494: if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high);
495: if (low != PETSC_DEFAULT) adapt->clip[0] = low;
496: if (high != PETSC_DEFAULT) adapt->clip[1] = high;
497: return(0);
498: }
500: /*@
501: TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
503: Not Collective
505: Input Arguments:
506: . adapt - adaptive controller context
508: Ouput Arguments:
509: + low - optional, admissible decrease factor
510: - high - optional, admissible increase factor
512: Level: intermediate
514: .seealso: TSAdaptChoose(), TSAdaptSetClip()
515: @*/
516: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
517: {
522: if (low) *low = adapt->clip[0];
523: if (high) *high = adapt->clip[1];
524: return(0);
525: }
527: /*@
528: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
530: Logically collective on TSAdapt
532: Input Arguments:
533: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
534: . hmin - minimum time step
535: - hmax - maximum time step
537: Options Database Keys:
538: + -ts_adapt_dt_min <min> - to set minimum time step
539: - -ts_adapt_dt_max <max> - to set maximum time step
541: Level: intermediate
543: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
544: @*/
545: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
546: {
552: if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
553: if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
554: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
555: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
556: hmin = adapt->dt_min;
557: hmax = adapt->dt_max;
558: 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);
559: return(0);
560: }
562: /*@
563: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
565: Not Collective
567: Input Arguments:
568: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
570: Output Arguments:
571: + hmin - minimum time step
572: - hmax - maximum time step
574: Level: intermediate
576: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
577: @*/
578: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
579: {
585: if (hmin) *hmin = adapt->dt_min;
586: if (hmax) *hmax = adapt->dt_max;
587: return(0);
588: }
590: /*
591: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
593: Collective on TSAdapt
595: Input Parameter:
596: . adapt - the TSAdapt context
598: Options Database Keys:
599: + -ts_adapt_type <type> - algorithm to use for adaptivity
600: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
601: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
602: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
603: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
604: . -ts_adapt_dt_min <min> - minimum timestep to use
605: . -ts_adapt_dt_max <max> - maximum timestep to use
606: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
607: . -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
608: - -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver
610: Level: advanced
612: Notes:
613: This function is automatically called by TSSetFromOptions()
615: .keywords: TS, TSGetAdapt(), TSAdaptSetType(), TSAdaptSetStepLimits()
617: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
618: TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
619: */
620: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
621: {
623: char type[256] = TSADAPTBASIC;
624: PetscReal safety,reject_safety,clip[2],hmin,hmax;
625: PetscBool set,flg;
626: PetscInt two;
630: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
631: * function can only be called from inside TSSetFromOptions() */
632: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
633: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
634: if (flg || !((PetscObject)adapt)->type_name) {
635: TSAdaptSetType(adapt,type);
636: }
638: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
639: if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}
641: safety = adapt->safety; reject_safety = adapt->reject_safety;
642: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
643: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
644: if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}
646: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
647: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
648: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
649: if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}
651: hmin = adapt->dt_min; hmax = adapt->dt_max;
652: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
653: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
654: if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}
656: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);
658: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
659: if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
661: 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);
662:
663: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
664: if (set) {TSAdaptSetMonitor(adapt,flg);}
666: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
667: PetscOptionsTail();
668: return(0);
669: }
671: /*@
672: TSAdaptCandidatesClear - clear any previously set candidate schemes
674: Logically collective on TSAdapt
676: Input Argument:
677: . adapt - adaptive controller
679: Level: developer
681: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
682: @*/
683: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
684: {
689: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
690: return(0);
691: }
693: /*@C
694: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
696: Logically collective on TSAdapt
698: Input Arguments:
699: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
700: . name - name of the candidate scheme to add
701: . order - order of the candidate scheme
702: . stageorder - stage order of the candidate scheme
703: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
704: . cost - relative measure of the amount of work required for the candidate scheme
705: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
707: Note:
708: This routine is not available in Fortran.
710: Level: developer
712: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
713: @*/
714: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
715: {
716: PetscInt c;
720: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
721: if (inuse) {
722: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
723: adapt->candidates.inuse_set = PETSC_TRUE;
724: }
725: /* first slot if this is the current scheme, otherwise the next available slot */
726: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
728: adapt->candidates.name[c] = name;
729: adapt->candidates.order[c] = order;
730: adapt->candidates.stageorder[c] = stageorder;
731: adapt->candidates.ccfl[c] = ccfl;
732: adapt->candidates.cost[c] = cost;
733: adapt->candidates.n++;
734: return(0);
735: }
737: /*@C
738: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
740: Not Collective
742: Input Arguments:
743: . adapt - time step adaptivity context
745: Output Arguments:
746: + n - number of candidate schemes, always at least 1
747: . order - the order of each candidate scheme
748: . stageorder - the stage order of each candidate scheme
749: . ccfl - the CFL coefficient of each scheme
750: - cost - the relative cost of each scheme
752: Level: developer
754: Note:
755: The current scheme is always returned in the first slot
757: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
758: @*/
759: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
760: {
763: if (n) *n = adapt->candidates.n;
764: if (order) *order = adapt->candidates.order;
765: if (stageorder) *stageorder = adapt->candidates.stageorder;
766: if (ccfl) *ccfl = adapt->candidates.ccfl;
767: if (cost) *cost = adapt->candidates.cost;
768: return(0);
769: }
771: /*@C
772: TSAdaptChoose - choose which method and step size to use for the next step
774: Collective on TSAdapt
776: Input Arguments:
777: + adapt - adaptive contoller
778: - h - current step size
780: Output Arguments:
781: + next_sc - optional, scheme to use for the next step
782: . next_h - step size to use for the next step
783: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
785: Note:
786: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
787: being retried after an initial rejection.
789: Level: developer
791: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
792: @*/
793: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
794: {
796: PetscInt ncandidates = adapt->candidates.n;
797: PetscInt scheme = 0;
798: PetscReal wlte = -1.0;
799: PetscReal wltea = -1.0;
800: PetscReal wlter = -1.0;
808: if (next_sc) *next_sc = 0;
810: /* Do not mess with adaptivity while handling events*/
811: if (ts->event && ts->event->status != TSEVENT_NONE) {
812: *next_h = h;
813: *accept = PETSC_TRUE;
814: return(0);
815: }
817: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
818: 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);
819: if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
820: if (next_sc) *next_sc = scheme;
822: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
823: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
824: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
825: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
826: PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
827: PetscReal b = adapt->matchstepfac[1];
828: if (t < tmax && tend > tmax) *next_h = hmax;
829: if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
830: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
831: }
833: if (adapt->monitor) {
834: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
835: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
836: if (wlte < 0) {
837: 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);
838: } else {
839: 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);
840: }
841: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
842: }
843: return(0);
844: }
846: /*@
847: TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
848: before increasing the time step.
850: Logicially Collective on TSAdapt
852: Input Arguments:
853: + adapt - adaptive controller context
854: - cnt - the number of timesteps
856: Options Database Key:
857: . -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase
859: Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
860: The successful use of this option is problem dependent
862: Developer Note: there is no theory to support this option
864: Level: advanced
866: .seealso:
867: @*/
868: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
869: {
871: adapt->timestepjustdecreased_delay = cnt;
872: return(0);
873: }
875:
876: /*@
877: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)
879: Collective on TSAdapt
881: Input Arguments:
882: + adapt - adaptive controller context
883: . ts - time stepper
884: . t - Current simulation time
885: - Y - Current solution vector
887: Output Arguments:
888: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
890: Level: developer
892: .seealso:
893: @*/
894: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
895: {
896: PetscErrorCode ierr;
897: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
904: if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
905: if (snesreason < 0) {
906: *accept = PETSC_FALSE;
907: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
908: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
909: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
910: if (adapt->monitor) {
911: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
912: 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);
913: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
914: }
915: }
916: } else {
917: *accept = PETSC_TRUE;
918: TSFunctionDomainError(ts,t,Y,accept);
919: if(*accept && adapt->checkstage) {
920: (*adapt->checkstage)(adapt,ts,t,Y,accept);
921: }
922: }
924: if(!(*accept) && !ts->reason) {
925: PetscReal dt,new_dt;
926: TSGetTimeStep(ts,&dt);
927: new_dt = dt * adapt->scale_solve_failed;
928: TSSetTimeStep(ts,new_dt);
929: adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
930: if (adapt->monitor) {
931: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
932: 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);
933: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
934: }
935: }
936: return(0);
937: }
939: /*@
940: TSAdaptCreate - create an adaptive controller context for time stepping
942: Collective on MPI_Comm
944: Input Parameter:
945: . comm - The communicator
947: Output Parameter:
948: . adapt - new TSAdapt object
950: Level: developer
952: Notes:
953: TSAdapt creation is handled by TS, so users should not need to call this function.
955: .keywords: TSAdapt, create
956: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
957: @*/
958: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
959: {
961: TSAdapt adapt;
965: *inadapt = NULL;
966: TSAdaptInitializePackage();
968: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
970: adapt->always_accept = PETSC_FALSE;
971: adapt->safety = 0.9;
972: adapt->reject_safety = 0.5;
973: adapt->clip[0] = 0.1;
974: adapt->clip[1] = 10.;
975: adapt->dt_min = 1e-20;
976: adapt->dt_max = 1e+20;
977: adapt->scale_solve_failed = 0.25;
978: /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
979: to prevent from situations were unreasonably small time steps are taken in order to match the final time */
980: adapt->matchstepfac[0] = 0.01; /* allow 1% step size increase in the last step */
981: adapt->matchstepfac[1] = 2.0; /* halve last step if it is greater than what remains divided this factor */
982: adapt->wnormtype = NORM_2;
983: adapt->timestepjustdecreased_delay = 0;
985: *inadapt = adapt;
986: return(0);
987: }