Actual source code: tsadapt.c
petsc-3.8.4 2018-03-24
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: 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: .keywords: TSAdapt, register, all
62: .seealso: TSAdaptRegisterDestroy()
63: @*/
64: PetscErrorCode TSAdaptRegisterAll(void)
65: {
69: if (TSAdaptRegisterAllCalled) return(0);
70: TSAdaptRegisterAllCalled = PETSC_TRUE;
71: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
72: TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);
73: TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP);
74: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
75: TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);
76: return(0);
77: }
79: /*@C
80: TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
81: called from PetscFinalize().
83: Level: developer
85: .keywords: Petsc, destroy, package
86: .seealso: PetscFinalize()
87: @*/
88: PetscErrorCode TSAdaptFinalizePackage(void)
89: {
93: PetscFunctionListDestroy(&TSAdaptList);
94: TSAdaptPackageInitialized = PETSC_FALSE;
95: TSAdaptRegisterAllCalled = PETSC_FALSE;
96: return(0);
97: }
99: /*@C
100: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
101: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
102: TSAdaptCreate() when using static libraries.
104: Level: developer
106: .keywords: TSAdapt, initialize, package
107: .seealso: PetscInitialize()
108: @*/
109: PetscErrorCode TSAdaptInitializePackage(void)
110: {
114: if (TSAdaptPackageInitialized) return(0);
115: TSAdaptPackageInitialized = PETSC_TRUE;
116: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
117: TSAdaptRegisterAll();
118: PetscRegisterFinalize(TSAdaptFinalizePackage);
119: return(0);
120: }
122: /*@C
123: TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE
125: Logicially Collective on TSAdapt
127: Input Parameter:
128: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
129: - type - either TSADAPTBASIC or TSADAPTNONE
131: Options Database:
132: . -ts_adapt_type basic or none - to setting the adapter type
134: Level: intermediate
136: .keywords: TSAdapt, create
138: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
139: @*/
140: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
141: {
142: PetscBool match;
143: PetscErrorCode ierr,(*r)(TSAdapt);
148: PetscObjectTypeCompare((PetscObject)adapt,type,&match);
149: if (match) return(0);
150: PetscFunctionListFind(TSAdaptList,type,&r);
151: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
152: if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
153: PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
154: PetscObjectChangeTypeName((PetscObject)adapt,type);
155: (*r)(adapt);
156: return(0);
157: }
159: /*@C
160: TSAdaptGetType - gets the TS adapter method type (as a string).
162: Not Collective
164: Input Parameter:
165: . adapt - The TS adapter, most likely obtained with TSGetAdapt()
167: Output Parameter:
168: . type - The name of TS adapter method
170: Level: intermediate
172: .keywords: TSAdapt, get, type
173: .seealso TSAdaptSetType()
174: @*/
175: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
176: {
180: *type = ((PetscObject)adapt)->type_name;
181: return(0);
182: }
184: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
185: {
190: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
191: return(0);
192: }
194: /*@C
195: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
197: Collective on PetscViewer
199: Input Parameters:
200: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
201: some related function before a call to TSAdaptLoad().
202: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
203: HDF5 file viewer, obtained from PetscViewerHDF5Open()
205: Level: intermediate
207: Notes:
208: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
210: Notes for advanced users:
211: Most users should not need to know the details of the binary storage
212: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
213: But for anyone who's interested, the standard binary matrix storage
214: format is
215: .vb
216: has not yet been determined
217: .ve
219: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
220: @*/
221: PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
222: {
224: PetscBool isbinary;
225: char type[256];
230: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
231: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
233: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
234: TSAdaptSetType(adapt,type);
235: if (adapt->ops->load) {
236: (*adapt->ops->load)(adapt,viewer);
237: }
238: return(0);
239: }
241: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
242: {
244: PetscBool iascii,isbinary,isnone;
248: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
251: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
252: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
253: if (iascii) {
254: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
255: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
256: if (!isnone) {
257: if (adapt->always_accept) {PetscViewerASCIIPrintf(viewer," always accepting steps\n");}
258: PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);
259: PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
260: PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);
261: PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);
262: PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);
263: PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);
264: }
265: if (adapt->ops->view) {
266: PetscViewerASCIIPushTab(viewer);
267: (*adapt->ops->view)(adapt,viewer);
268: PetscViewerASCIIPopTab(viewer);
269: }
270: } else if (isbinary) {
271: char type[256];
273: /* need to save FILE_CLASS_ID for adapt class */
274: PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
275: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
276: } else if (adapt->ops->view) {
277: (*adapt->ops->view)(adapt,viewer);
278: }
279: return(0);
280: }
282: /*@
283: TSAdaptReset - Resets a TSAdapt context.
285: Collective on TS
287: Input Parameter:
288: . adapt - the TSAdapt context obtained from TSAdaptCreate()
290: Level: developer
292: .seealso: TSAdaptCreate(), TSAdaptDestroy()
293: @*/
294: PetscErrorCode TSAdaptReset(TSAdapt adapt)
295: {
300: if (adapt->ops->reset) {(*adapt->ops->reset)(adapt);}
301: return(0);
302: }
304: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
305: {
309: if (!*adapt) return(0);
311: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return(0);}
313: TSAdaptReset(*adapt);
315: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
316: PetscViewerDestroy(&(*adapt)->monitor);
317: PetscHeaderDestroy(adapt);
318: return(0);
319: }
321: /*@
322: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
324: Collective on TSAdapt
326: Input Arguments:
327: + adapt - adaptive controller context
328: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
330: Options Database Keys:
331: . -ts_adapt_monitor
333: Level: intermediate
335: .seealso: TSAdaptChoose()
336: @*/
337: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
338: {
344: if (flg) {
345: if (!adapt->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);}
346: } else {
347: PetscViewerDestroy(&adapt->monitor);
348: }
349: return(0);
350: }
352: /*@C
353: TSAdaptSetCheckStage - Set a callback to check convergence for a stage
355: Logically collective on TSAdapt
357: Input Arguments:
358: + adapt - adaptive controller context
359: - func - stage check function
361: Arguments of func:
362: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
364: + adapt - adaptive controller context
365: . ts - time stepping context
366: - accept - pending choice of whether to accept, can be modified by this routine
368: Level: advanced
370: .seealso: TSAdaptChoose()
371: @*/
372: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*))
373: {
377: adapt->checkstage = func;
378: return(0);
379: }
381: /*@
382: TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of
383: any error or stability condition not meeting the prescribed goal.
385: Logically collective on TSAdapt
387: Input Arguments:
388: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
389: - flag - whether to always accept steps
391: Options Database Keys:
392: . -ts_adapt_always_accept
394: Level: intermediate
396: .seealso: TSAdapt, TSAdaptChoose()
397: @*/
398: PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)
399: {
403: adapt->always_accept = flag;
404: return(0);
405: }
407: /*@
408: TSAdaptSetSafety - Set safety factors
410: Logically collective on TSAdapt
412: Input Arguments:
413: + adapt - adaptive controller context
414: . safety - safety factor relative to target error/stability goal
415: + reject_safety - extra safety factor to apply if the last step was rejected
417: Options Database Keys:
418: + -ts_adapt_safety
419: - -ts_adapt_reject_safety
421: Level: intermediate
423: .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose()
424: @*/
425: PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)
426: {
431: if (safety != PETSC_DEFAULT && safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be non negative",(double)safety);
432: if (safety != PETSC_DEFAULT && safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be less than one",(double)safety);
433: 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);
434: 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);
435: if (safety != PETSC_DEFAULT) adapt->safety = safety;
436: if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety;
437: return(0);
438: }
440: /*@
441: TSAdaptGetSafety - Get safety factors
443: Not Collective
445: Input Arguments:
446: . adapt - adaptive controller context
448: Ouput Arguments:
449: . safety - safety factor relative to target error/stability goal
450: + reject_safety - extra safety factor to apply if the last step was rejected
452: Level: intermediate
454: .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose()
455: @*/
456: PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety)
457: {
462: if (safety) *safety = adapt->safety;
463: if (reject_safety) *reject_safety = adapt->reject_safety;
464: return(0);
465: }
467: /*@
468: TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
470: Logically collective on TSAdapt
472: Input Arguments:
473: + adapt - adaptive controller context
474: . low - admissible decrease factor
475: + high - admissible increase factor
477: Options Database Keys:
478: . -ts_adapt_clip
480: Level: intermediate
482: .seealso: TSAdaptChoose(), TSAdaptGetClip()
483: @*/
484: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
485: {
490: if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
491: if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
492: if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high);
493: if (low != PETSC_DEFAULT) adapt->clip[0] = low;
494: if (high != PETSC_DEFAULT) adapt->clip[1] = high;
495: return(0);
496: }
498: /*@
499: TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
501: Not Collective
503: Input Arguments:
504: . adapt - adaptive controller context
506: Ouput Arguments:
507: + low - optional, admissible decrease factor
508: - high - optional, admissible increase factor
510: Level: intermediate
512: .seealso: TSAdaptChoose(), TSAdaptSetClip()
513: @*/
514: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
515: {
520: if (low) *low = adapt->clip[0];
521: if (high) *high = adapt->clip[1];
522: return(0);
523: }
525: /*@
526: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
528: Logically collective on TSAdapt
530: Input Arguments:
531: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
532: . hmin - minimum time step
533: - hmax - maximum time step
535: Options Database Keys:
536: + -ts_adapt_dt_min - minimum time step
537: - -ts_adapt_dt_max - maximum time step
539: Level: intermediate
541: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
542: @*/
543: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
544: {
550: if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
551: if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
552: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
553: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
554: hmin = adapt->dt_min;
555: hmax = adapt->dt_max;
556: 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);
557: return(0);
558: }
560: /*@
561: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
563: Not Collective
565: Input Arguments:
566: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
568: Output Arguments:
569: + hmin - minimum time step
570: - hmax - maximum time step
572: Level: intermediate
574: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
575: @*/
576: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
577: {
583: if (hmin) *hmin = adapt->dt_min;
584: if (hmax) *hmax = adapt->dt_max;
585: return(0);
586: }
588: /*
589: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
591: Collective on TSAdapt
593: Input Parameter:
594: . adapt - the TSAdapt context
596: Options Database Keys:
597: + -ts_adapt_type <type> - algorithm to use for adaptivity
598: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
599: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
600: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
601: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
602: . -ts_adapt_dt_min <min> - minimum timestep to use
603: . -ts_adapt_dt_max <max> - maximum timestep to use
604: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
605: - -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
607: Level: advanced
609: Notes:
610: This function is automatically called by TSSetFromOptions()
612: .keywords: TS, TSGetAdapt(), TSAdaptSetType(), TSAdaptSetStepLimits()
614: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
615: TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
616: */
617: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
618: {
620: char type[256] = TSADAPTBASIC;
621: PetscReal safety,reject_safety,clip[2],hmin,hmax;
622: PetscBool set,flg;
623: PetscInt two;
627: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
628: * function can only be called from inside TSSetFromOptions() */
629: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
630: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
631: if (flg || !((PetscObject)adapt)->type_name) {
632: TSAdaptSetType(adapt,type);
633: }
635: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
636: if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}
638: safety = adapt->safety; reject_safety = adapt->reject_safety;
639: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
640: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
641: if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}
643: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
644: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
645: if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
646: if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}
648: hmin = adapt->dt_min; hmax = adapt->dt_max;
649: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
650: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
651: if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}
653: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);
655: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
656: if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
658: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
659: if (set) {TSAdaptSetMonitor(adapt,flg);}
661: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
662: PetscOptionsTail();
663: return(0);
664: }
666: /*@
667: TSAdaptCandidatesClear - clear any previously set candidate schemes
669: Logically collective on TSAdapt
671: Input Argument:
672: . adapt - adaptive controller
674: Level: developer
676: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
677: @*/
678: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
679: {
684: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
685: return(0);
686: }
688: /*@C
689: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
691: Logically collective on TSAdapt
693: Input Arguments:
694: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
695: . name - name of the candidate scheme to add
696: . order - order of the candidate scheme
697: . stageorder - stage order of the candidate scheme
698: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
699: . cost - relative measure of the amount of work required for the candidate scheme
700: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
702: Note:
703: This routine is not available in Fortran.
705: Level: developer
707: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
708: @*/
709: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
710: {
711: PetscInt c;
715: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
716: if (inuse) {
717: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
718: adapt->candidates.inuse_set = PETSC_TRUE;
719: }
720: /* first slot if this is the current scheme, otherwise the next available slot */
721: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
723: adapt->candidates.name[c] = name;
724: adapt->candidates.order[c] = order;
725: adapt->candidates.stageorder[c] = stageorder;
726: adapt->candidates.ccfl[c] = ccfl;
727: adapt->candidates.cost[c] = cost;
728: adapt->candidates.n++;
729: return(0);
730: }
732: /*@C
733: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
735: Not Collective
737: Input Arguments:
738: . adapt - time step adaptivity context
740: Output Arguments:
741: + n - number of candidate schemes, always at least 1
742: . order - the order of each candidate scheme
743: . stageorder - the stage order of each candidate scheme
744: . ccfl - the CFL coefficient of each scheme
745: - cost - the relative cost of each scheme
747: Level: developer
749: Note:
750: The current scheme is always returned in the first slot
752: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
753: @*/
754: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
755: {
758: if (n) *n = adapt->candidates.n;
759: if (order) *order = adapt->candidates.order;
760: if (stageorder) *stageorder = adapt->candidates.stageorder;
761: if (ccfl) *ccfl = adapt->candidates.ccfl;
762: if (cost) *cost = adapt->candidates.cost;
763: return(0);
764: }
766: /*@C
767: TSAdaptChoose - choose which method and step size to use for the next step
769: Collective on TSAdapt
771: Input Arguments:
772: + adapt - adaptive contoller
773: - h - current step size
775: Output Arguments:
776: + next_sc - optional, scheme to use for the next step
777: . next_h - step size to use for the next step
778: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
780: Note:
781: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
782: being retried after an initial rejection.
784: Level: developer
786: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
787: @*/
788: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
789: {
791: PetscInt ncandidates = adapt->candidates.n;
792: PetscInt scheme = 0;
793: PetscReal wlte = -1.0;
794: PetscReal wltea = -1.0;
795: PetscReal wlter = -1.0;
803: if (next_sc) *next_sc = 0;
805: /* Do not mess with adaptivity while handling events*/
806: if (ts->event && ts->event->status != TSEVENT_NONE) {
807: *next_h = h;
808: *accept = PETSC_TRUE;
809: return(0);
810: }
812: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
813: 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);
814: if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
815: if (next_sc) *next_sc = scheme;
817: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
818: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
819: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
820: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
821: PetscReal a = (PetscReal)1.01; /* allow 1% step size increase */
822: if (t < tmax && tend > tmax) *next_h = hmax;
823: if (t < tmax && tend < tmax && h > hmax/2) *next_h = hmax/2;
824: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
825: }
827: if (adapt->monitor) {
828: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
829: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
830: if (wlte < 0) {
831: 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);
832: } else {
833: 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);
834: }
835: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
836: }
837: return(0);
838: }
840: /*@
841: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)
843: Collective on TSAdapt
845: Input Arguments:
846: + adapt - adaptive controller context
847: . ts - time stepper
848: . t - Current simulation time
849: - Y - Current solution vector
851: Output Arguments:
852: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
854: Level: developer
856: .seealso:
857: @*/
858: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
859: {
860: PetscErrorCode ierr;
861: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
868: if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
869: if (snesreason < 0) {
870: *accept = PETSC_FALSE;
871: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
872: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
873: PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
874: if (adapt->monitor) {
875: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
876: 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);
877: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
878: }
879: }
880: } else {
881: *accept = PETSC_TRUE;
882: TSFunctionDomainError(ts,t,Y,accept);
883: if(*accept && adapt->checkstage) {
884: (*adapt->checkstage)(adapt,ts,t,Y,accept);
885: }
886: }
888: if(!(*accept) && !ts->reason) {
889: PetscReal dt,new_dt;
890: TSGetTimeStep(ts,&dt);
891: new_dt = dt * adapt->scale_solve_failed;
892: TSSetTimeStep(ts,new_dt);
893: adapt->timestepjustincreased += 4;
894: if (adapt->monitor) {
895: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
896: 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);
897: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
898: }
899: }
900: return(0);
901: }
903: /*@
904: TSAdaptCreate - create an adaptive controller context for time stepping
906: Collective on MPI_Comm
908: Input Parameter:
909: . comm - The communicator
911: Output Parameter:
912: . adapt - new TSAdapt object
914: Level: developer
916: Notes:
917: TSAdapt creation is handled by TS, so users should not need to call this function.
919: .keywords: TSAdapt, create
920: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
921: @*/
922: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
923: {
925: TSAdapt adapt;
929: *inadapt = NULL;
930: TSAdaptInitializePackage();
932: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
934: adapt->always_accept = PETSC_FALSE;
935: adapt->safety = 0.9;
936: adapt->reject_safety = 0.5;
937: adapt->clip[0] = 0.1;
938: adapt->clip[1] = 10.;
939: adapt->dt_min = 1e-20;
940: adapt->dt_max = 1e+20;
941: adapt->scale_solve_failed = 0.25;
942: adapt->wnormtype = NORM_2;
944: *inadapt = adapt;
945: return(0);
946: }