Actual source code: tsadapt.c
2: #include <petsc/private/tsimpl.h>
4: PetscClassId TSADAPT_CLASSID;
6: static PetscFunctionList TSAdaptList;
7: static PetscBool TSAdaptPackageInitialized;
8: static PetscBool TSAdaptRegisterAllCalled;
10: PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt);
11: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt);
12: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt);
13: PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt);
14: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt);
15: PETSC_EXTERN PetscErrorCode TSAdaptCreate_History(TSAdapt);
17: /*@C
18: TSAdaptRegister - adds a TSAdapt implementation
20: Not Collective
22: Input Parameters:
23: + name_scheme - name of user-defined adaptivity scheme
24: - routine_create - routine to create method context
26: Notes:
27: TSAdaptRegister() may be called multiple times to add several user-defined families.
29: Sample usage:
30: .vb
31: TSAdaptRegister("my_scheme",MySchemeCreate);
32: .ve
34: Then, your scheme can be chosen with the procedural interface via
35: $ TSAdaptSetType(ts,"my_scheme")
36: or at runtime via the option
37: $ -ts_adapt_type my_scheme
39: Level: advanced
41: .seealso: TSAdaptRegisterAll()
42: @*/
43: PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
44: {
45: TSAdaptInitializePackage();
46: PetscFunctionListAdd(&TSAdaptList,sname,function);
47: return 0;
48: }
50: /*@C
51: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
53: Not Collective
55: Level: advanced
57: .seealso: TSAdaptRegisterDestroy()
58: @*/
59: PetscErrorCode TSAdaptRegisterAll(void)
60: {
61: if (TSAdaptRegisterAllCalled) return 0;
62: TSAdaptRegisterAllCalled = PETSC_TRUE;
63: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
64: TSAdaptRegister(TSADAPTBASIC, TSAdaptCreate_Basic);
65: TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP);
66: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
67: TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);
68: TSAdaptRegister(TSADAPTHISTORY,TSAdaptCreate_History);
69: return 0;
70: }
72: /*@C
73: TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
74: called from PetscFinalize().
76: Level: developer
78: .seealso: PetscFinalize()
79: @*/
80: PetscErrorCode TSAdaptFinalizePackage(void)
81: {
82: PetscFunctionListDestroy(&TSAdaptList);
83: TSAdaptPackageInitialized = PETSC_FALSE;
84: TSAdaptRegisterAllCalled = PETSC_FALSE;
85: return 0;
86: }
88: /*@C
89: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
90: called from TSInitializePackage().
92: Level: developer
94: .seealso: PetscInitialize()
95: @*/
96: PetscErrorCode TSAdaptInitializePackage(void)
97: {
98: if (TSAdaptPackageInitialized) return 0;
99: TSAdaptPackageInitialized = PETSC_TRUE;
100: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
101: TSAdaptRegisterAll();
102: PetscRegisterFinalize(TSAdaptFinalizePackage);
103: return 0;
104: }
106: /*@C
107: TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE
109: Logicially Collective on TSAdapt
111: Input Parameters:
112: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
113: - type - either TSADAPTBASIC or TSADAPTNONE
115: Options Database:
116: . -ts_adapt_type <basic or dsp or none> - to set the adapter type
118: Level: intermediate
120: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
121: @*/
122: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
123: {
124: PetscBool match;
125: PetscErrorCode (*r)(TSAdapt);
129: PetscObjectTypeCompare((PetscObject)adapt,type,&match);
130: if (match) return 0;
131: PetscFunctionListFind(TSAdaptList,type,&r);
133: if (adapt->ops->destroy) (*adapt->ops->destroy)(adapt);
134: PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
135: PetscObjectChangeTypeName((PetscObject)adapt,type);
136: (*r)(adapt);
137: return 0;
138: }
140: /*@C
141: TSAdaptGetType - gets the TS adapter method type (as a string).
143: Not Collective
145: Input Parameter:
146: . adapt - The TS adapter, most likely obtained with TSGetAdapt()
148: Output Parameter:
149: . type - The name of TS adapter method
151: Level: intermediate
153: .seealso TSAdaptSetType()
154: @*/
155: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
156: {
159: *type = ((PetscObject)adapt)->type_name;
160: return 0;
161: }
163: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
164: {
166: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
167: return 0;
168: }
170: /*@C
171: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
173: Collective on PetscViewer
175: Input Parameters:
176: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
177: some related function before a call to TSAdaptLoad().
178: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
179: HDF5 file viewer, obtained from PetscViewerHDF5Open()
181: Level: intermediate
183: Notes:
184: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
186: Notes for advanced users:
187: Most users should not need to know the details of the binary storage
188: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
189: But for anyone who's interested, the standard binary matrix storage
190: format is
191: .vb
192: has not yet been determined
193: .ve
195: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
196: @*/
197: PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
198: {
199: PetscBool isbinary;
200: char type[256];
204: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
207: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
208: TSAdaptSetType(adapt,type);
209: if (adapt->ops->load) (*adapt->ops->load)(adapt,viewer);
210: return 0;
211: }
213: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
214: {
215: PetscBool iascii,isbinary,isnone,isglee;
218: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);
221: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
222: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
223: if (iascii) {
224: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
225: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
226: PetscObjectTypeCompare((PetscObject)adapt,TSADAPTGLEE,&isglee);
227: if (!isnone) {
228: if (adapt->always_accept) PetscViewerASCIIPrintf(viewer," always accepting steps\n");
229: PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);
230: PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
231: PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);
232: PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);
233: PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);
234: PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);
235: PetscViewerASCIIPrintf(viewer," maximum solution absolute value to be ignored %g\n",(double)adapt->ignore_max);
236: }
237: if (isglee) {
238: if (adapt->glee_use_local) {
239: PetscViewerASCIIPrintf(viewer," GLEE uses local error control\n");
240: } else {
241: PetscViewerASCIIPrintf(viewer," GLEE uses global error control\n");
242: }
243: }
244: if (adapt->ops->view) {
245: PetscViewerASCIIPushTab(viewer);
246: (*adapt->ops->view)(adapt,viewer);
247: PetscViewerASCIIPopTab(viewer);
248: }
249: } else if (isbinary) {
250: char type[256];
252: /* need to save FILE_CLASS_ID for adapt class */
253: PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
254: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
255: } else if (adapt->ops->view) {
256: (*adapt->ops->view)(adapt,viewer);
257: }
258: return 0;
259: }
261: /*@
262: TSAdaptReset - Resets a TSAdapt context.
264: Collective on TS
266: Input Parameter:
267: . adapt - the TSAdapt context obtained from TSAdaptCreate()
269: Level: developer
271: .seealso: TSAdaptCreate(), TSAdaptDestroy()
272: @*/
273: PetscErrorCode TSAdaptReset(TSAdapt adapt)
274: {
276: if (adapt->ops->reset) (*adapt->ops->reset)(adapt);
277: return 0;
278: }
280: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
281: {
282: if (!*adapt) return 0;
284: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return 0;}
286: TSAdaptReset(*adapt);
288: if ((*adapt)->ops->destroy) (*(*adapt)->ops->destroy)(*adapt);
289: PetscViewerDestroy(&(*adapt)->monitor);
290: PetscHeaderDestroy(adapt);
291: return 0;
292: }
294: /*@
295: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
297: Collective on TSAdapt
299: Input Parameters:
300: + adapt - adaptive controller context
301: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
303: Options Database Keys:
304: . -ts_adapt_monitor - to turn on monitoring
306: Level: intermediate
308: .seealso: TSAdaptChoose()
309: @*/
310: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
311: {
314: if (flg) {
315: if (!adapt->monitor) PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);
316: } else {
317: PetscViewerDestroy(&adapt->monitor);
318: }
319: return 0;
320: }
322: /*@C
323: TSAdaptSetCheckStage - Set a callback to check convergence for a stage
325: Logically collective on TSAdapt
327: Input Parameters:
328: + adapt - adaptive controller context
329: - func - stage check function
331: Arguments of func:
332: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
334: + adapt - adaptive controller context
335: . ts - time stepping context
336: - accept - pending choice of whether to accept, can be modified by this routine
338: Level: advanced
340: .seealso: TSAdaptChoose()
341: @*/
342: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*))
343: {
345: adapt->checkstage = func;
346: return 0;
347: }
349: /*@
350: TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of
351: any error or stability condition not meeting the prescribed goal.
353: Logically collective on TSAdapt
355: Input Parameters:
356: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
357: - flag - whether to always accept steps
359: Options Database Keys:
360: . -ts_adapt_always_accept - to always accept steps
362: Level: intermediate
364: .seealso: TSAdapt, TSAdaptChoose()
365: @*/
366: PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)
367: {
370: adapt->always_accept = flag;
371: return 0;
372: }
374: /*@
375: TSAdaptSetSafety - Set safety factors
377: Logically collective on TSAdapt
379: Input Parameters:
380: + adapt - adaptive controller context
381: . safety - safety factor relative to target error/stability goal
382: - reject_safety - extra safety factor to apply if the last step was rejected
384: Options Database Keys:
385: + -ts_adapt_safety <safety> - to set safety factor
386: - -ts_adapt_reject_safety <reject_safety> - to set reject safety factor
388: Level: intermediate
390: .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose()
391: @*/
392: PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)
393: {
401: if (safety != PETSC_DEFAULT) adapt->safety = safety;
402: if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety;
403: return 0;
404: }
406: /*@
407: TSAdaptGetSafety - Get safety factors
409: Not Collective
411: Input Parameter:
412: . adapt - adaptive controller context
414: Output Parameters:
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: Level: intermediate
420: .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose()
421: @*/
422: PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety)
423: {
427: if (safety) *safety = adapt->safety;
428: if (reject_safety) *reject_safety = adapt->reject_safety;
429: return 0;
430: }
432: /*@
433: TSAdaptSetMaxIgnore - Set error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value). A negative value (default) of the threshold leads to considering all solution components.
435: Logically collective on TSAdapt
437: Input Parameters:
438: + adapt - adaptive controller context
439: - max_ignore - threshold for solution components that are ignored during error estimation
441: Options Database Keys:
442: . -ts_adapt_max_ignore <max_ignore> - to set the threshold
444: Level: intermediate
446: .seealso: TSAdapt, TSAdaptGetMaxIgnore(), TSAdaptChoose()
447: @*/
448: PetscErrorCode TSAdaptSetMaxIgnore(TSAdapt adapt,PetscReal max_ignore)
449: {
452: adapt->ignore_max = max_ignore;
453: return 0;
454: }
456: /*@
457: TSAdaptGetMaxIgnore - Get error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value).
459: Not Collective
461: Input Parameter:
462: . adapt - adaptive controller context
464: Output Parameter:
465: . max_ignore - threshold for solution components that are ignored during error estimation
467: Level: intermediate
469: .seealso: TSAdapt, TSAdaptSetMaxIgnore(), TSAdaptChoose()
470: @*/
471: PetscErrorCode TSAdaptGetMaxIgnore(TSAdapt adapt,PetscReal *max_ignore)
472: {
475: *max_ignore = adapt->ignore_max;
476: return 0;
477: }
479: /*@
480: TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
482: Logically collective on TSAdapt
484: Input Parameters:
485: + adapt - adaptive controller context
486: . low - admissible decrease factor
487: - high - admissible increase factor
489: Options Database Keys:
490: . -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors
492: Level: intermediate
494: .seealso: TSAdaptChoose(), TSAdaptGetClip(), TSAdaptSetScaleSolveFailed()
495: @*/
496: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
497: {
504: if (low != PETSC_DEFAULT) adapt->clip[0] = low;
505: if (high != PETSC_DEFAULT) adapt->clip[1] = high;
506: return 0;
507: }
509: /*@
510: TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
512: Not Collective
514: Input Parameter:
515: . adapt - adaptive controller context
517: Output Parameters:
518: + low - optional, admissible decrease factor
519: - high - optional, admissible increase factor
521: Level: intermediate
523: .seealso: TSAdaptChoose(), TSAdaptSetClip(), TSAdaptSetScaleSolveFailed()
524: @*/
525: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
526: {
530: if (low) *low = adapt->clip[0];
531: if (high) *high = adapt->clip[1];
532: return 0;
533: }
535: /*@
536: TSAdaptSetScaleSolveFailed - Scale step by this factor if solve fails
538: Logically collective on TSAdapt
540: Input Parameters:
541: + adapt - adaptive controller context
542: - scale - scale
544: Options Database Keys:
545: . -ts_adapt_scale_solve_failed <scale> - to set scale step by this factor if solve fails
547: Level: intermediate
549: .seealso: TSAdaptChoose(), TSAdaptGetScaleSolveFailed(), TSAdaptGetClip()
550: @*/
551: PetscErrorCode TSAdaptSetScaleSolveFailed(TSAdapt adapt,PetscReal scale)
552: {
557: if (scale != PETSC_DEFAULT) adapt->scale_solve_failed = scale;
558: return 0;
559: }
561: /*@
562: TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size
564: Not Collective
566: Input Parameter:
567: . adapt - adaptive controller context
569: Output Parameter:
570: . scale - scale factor
572: Level: intermediate
574: .seealso: TSAdaptChoose(), TSAdaptSetScaleSolveFailed(), TSAdaptSetClip()
575: @*/
576: PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal *scale)
577: {
580: if (scale) *scale = adapt->scale_solve_failed;
581: return 0;
582: }
584: /*@
585: TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
587: Logically collective on TSAdapt
589: Input Parameters:
590: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
591: . hmin - minimum time step
592: - hmax - maximum time step
594: Options Database Keys:
595: + -ts_adapt_dt_min <min> - to set minimum time step
596: - -ts_adapt_dt_max <max> - to set maximum time step
598: Level: intermediate
600: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
601: @*/
602: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
603: {
609: if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
610: if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
611: hmin = adapt->dt_min;
612: hmax = adapt->dt_max;
614: return 0;
615: }
617: /*@
618: TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
620: Not Collective
622: Input Parameter:
623: . adapt - time step adaptivity context, usually gotten with TSGetAdapt()
625: Output Parameters:
626: + hmin - minimum time step
627: - hmax - maximum time step
629: Level: intermediate
631: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
632: @*/
633: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
634: {
638: if (hmin) *hmin = adapt->dt_min;
639: if (hmax) *hmax = adapt->dt_max;
640: return 0;
641: }
643: /*
644: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
646: Collective on TSAdapt
648: Input Parameter:
649: . adapt - the TSAdapt context
651: Options Database Keys:
652: + -ts_adapt_type <type> - algorithm to use for adaptivity
653: . -ts_adapt_always_accept - always accept steps regardless of error/stability goals
654: . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
655: . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
656: . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
657: . -ts_adapt_dt_min <min> - minimum timestep to use
658: . -ts_adapt_dt_max <max> - maximum timestep to use
659: . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
660: . -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
661: - -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver
663: Level: advanced
665: Notes:
666: This function is automatically called by TSSetFromOptions()
668: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
669: TSAdaptSetClip(), TSAdaptSetScaleSolveFailed(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
670: */
671: PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
672: {
673: char type[256] = TSADAPTBASIC;
674: PetscReal safety,reject_safety,clip[2],scale,hmin,hmax;
675: PetscBool set,flg;
676: PetscInt two;
679: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
680: * function can only be called from inside TSSetFromOptions() */
681: PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
682: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
683: if (flg || !((PetscObject)adapt)->type_name) {
684: TSAdaptSetType(adapt,type);
685: }
687: PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
688: if (set) TSAdaptSetAlwaysAccept(adapt,flg);
690: safety = adapt->safety; reject_safety = adapt->reject_safety;
691: PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
692: PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
693: if (set || flg) TSAdaptSetSafety(adapt,safety,reject_safety);
695: two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
696: PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
698: if (set) TSAdaptSetClip(adapt,clip[0],clip[1]);
700: hmin = adapt->dt_min; hmax = adapt->dt_max;
701: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
702: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
703: if (set || flg) TSAdaptSetStepLimits(adapt,hmin,hmax);
705: PetscOptionsReal("-ts_adapt_max_ignore","Adaptor ignores (absolute) solution values smaller than this value","",adapt->ignore_max,&adapt->ignore_max,&set);
706: PetscOptionsBool("-ts_adapt_glee_use_local","GLEE adaptor uses local error estimation for step control","",adapt->glee_use_local,&adapt->glee_use_local,&set);
708: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","TSAdaptSetScaleSolveFailed",adapt->scale_solve_failed,&scale,&set);
709: if (set) TSAdaptSetScaleSolveFailed(adapt,scale);
711: PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
714: 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);
716: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
717: if (set) TSAdaptSetMonitor(adapt,flg);
719: if (adapt->ops->setfromoptions) (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);
720: PetscOptionsTail();
721: return 0;
722: }
724: /*@
725: TSAdaptCandidatesClear - clear any previously set candidate schemes
727: Logically collective on TSAdapt
729: Input Parameter:
730: . adapt - adaptive controller
732: Level: developer
734: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
735: @*/
736: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
737: {
739: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
740: return 0;
741: }
743: /*@C
744: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
746: Logically collective on TSAdapt
748: Input Parameters:
749: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
750: . name - name of the candidate scheme to add
751: . order - order of the candidate scheme
752: . stageorder - stage order of the candidate scheme
753: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
754: . cost - relative measure of the amount of work required for the candidate scheme
755: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
757: Note:
758: This routine is not available in Fortran.
760: Level: developer
762: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
763: @*/
764: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
765: {
766: PetscInt c;
770: if (inuse) {
772: adapt->candidates.inuse_set = PETSC_TRUE;
773: }
774: /* first slot if this is the current scheme, otherwise the next available slot */
775: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
777: adapt->candidates.name[c] = name;
778: adapt->candidates.order[c] = order;
779: adapt->candidates.stageorder[c] = stageorder;
780: adapt->candidates.ccfl[c] = ccfl;
781: adapt->candidates.cost[c] = cost;
782: adapt->candidates.n++;
783: return 0;
784: }
786: /*@C
787: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
789: Not Collective
791: Input Parameter:
792: . adapt - time step adaptivity context
794: Output Parameters:
795: + n - number of candidate schemes, always at least 1
796: . order - the order of each candidate scheme
797: . stageorder - the stage order of each candidate scheme
798: . ccfl - the CFL coefficient of each scheme
799: - cost - the relative cost of each scheme
801: Level: developer
803: Note:
804: The current scheme is always returned in the first slot
806: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
807: @*/
808: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
809: {
811: if (n) *n = adapt->candidates.n;
812: if (order) *order = adapt->candidates.order;
813: if (stageorder) *stageorder = adapt->candidates.stageorder;
814: if (ccfl) *ccfl = adapt->candidates.ccfl;
815: if (cost) *cost = adapt->candidates.cost;
816: return 0;
817: }
819: /*@C
820: TSAdaptChoose - choose which method and step size to use for the next step
822: Collective on TSAdapt
824: Input Parameters:
825: + adapt - adaptive contoller
826: . ts - time stepper
827: - h - current step size
829: Output Parameters:
830: + next_sc - optional, scheme to use for the next step
831: . next_h - step size to use for the next step
832: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
834: Note:
835: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
836: being retried after an initial rejection.
838: Level: developer
840: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
841: @*/
842: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
843: {
844: PetscInt ncandidates = adapt->candidates.n;
845: PetscInt scheme = 0;
846: PetscReal wlte = -1.0;
847: PetscReal wltea = -1.0;
848: PetscReal wlter = -1.0;
855: if (next_sc) *next_sc = 0;
857: /* Do not mess with adaptivity while handling events*/
858: if (ts->event && ts->event->status != TSEVENT_NONE) {
859: *next_h = h;
860: *accept = PETSC_TRUE;
861: return 0;
862: }
864: (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
867: if (next_sc) *next_sc = scheme;
869: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
870: /* Increase/reduce step size if end time of next step is close to or overshoots max time */
871: PetscReal t = ts->ptime + ts->time_step, h = *next_h;
872: PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
873: PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
874: PetscReal b = adapt->matchstepfac[1];
875: if (t < tmax && tend > tmax) *next_h = hmax;
876: if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
877: if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
878: }
880: if (adapt->monitor) {
881: const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
882: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
883: if (wlte < 0) {
884: 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);
885: } else {
886: 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);
887: }
888: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
889: }
890: return 0;
891: }
893: /*@
894: TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
895: before increasing the time step.
897: Logicially Collective on TSAdapt
899: Input Parameters:
900: + adapt - adaptive controller context
901: - cnt - the number of timesteps
903: Options Database Key:
904: . -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase
906: Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
907: The successful use of this option is problem dependent
909: Developer Note: there is no theory to support this option
911: Level: advanced
913: .seealso:
914: @*/
915: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
916: {
917: adapt->timestepjustdecreased_delay = cnt;
918: return 0;
919: }
921: /*@
922: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible)
924: Collective on TSAdapt
926: Input Parameters:
927: + adapt - adaptive controller context
928: . ts - time stepper
929: . t - Current simulation time
930: - Y - Current solution vector
932: Output Parameter:
933: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
935: Level: developer
937: .seealso:
938: @*/
939: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
940: {
941: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
947: if (ts->snes) SNESGetConvergedReason(ts->snes,&snesreason);
948: if (snesreason < 0) {
949: *accept = PETSC_FALSE;
950: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
951: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
952: PetscInfo(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
953: if (adapt->monitor) {
954: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
955: 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);
956: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
957: }
958: }
959: } else {
960: *accept = PETSC_TRUE;
961: TSFunctionDomainError(ts,t,Y,accept);
962: if (*accept && adapt->checkstage) {
963: (*adapt->checkstage)(adapt,ts,t,Y,accept);
964: if (!*accept) {
965: PetscInfo(ts,"Step=%D, solution rejected by user function provided by TSSetFunctionDomainError()\n",ts->steps);
966: if (adapt->monitor) {
967: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
968: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected by user function provided by TSSetFunctionDomainError()\n",((PetscObject)adapt)->type_name,ts->steps);
969: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
970: }
971: }
972: }
973: }
975: if (!(*accept) && !ts->reason) {
976: PetscReal dt,new_dt;
977: TSGetTimeStep(ts,&dt);
978: new_dt = dt * adapt->scale_solve_failed;
979: TSSetTimeStep(ts,new_dt);
980: adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
981: if (adapt->monitor) {
982: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
983: 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);
984: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
985: }
986: }
987: return 0;
988: }
990: /*@
991: TSAdaptCreate - create an adaptive controller context for time stepping
993: Collective
995: Input Parameter:
996: . comm - The communicator
998: Output Parameter:
999: . adapt - new TSAdapt object
1001: Level: developer
1003: Notes:
1004: TSAdapt creation is handled by TS, so users should not need to call this function.
1006: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1007: @*/
1008: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1009: {
1010: TSAdapt adapt;
1013: *inadapt = NULL;
1014: TSAdaptInitializePackage();
1016: PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
1018: adapt->always_accept = PETSC_FALSE;
1019: adapt->safety = 0.9;
1020: adapt->reject_safety = 0.5;
1021: adapt->clip[0] = 0.1;
1022: adapt->clip[1] = 10.;
1023: adapt->dt_min = 1e-20;
1024: adapt->dt_max = 1e+20;
1025: adapt->ignore_max = -1.0;
1026: adapt->glee_use_local = PETSC_TRUE;
1027: adapt->scale_solve_failed = 0.25;
1028: /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1029: to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1030: adapt->matchstepfac[0] = 0.01; /* allow 1% step size increase in the last step */
1031: adapt->matchstepfac[1] = 2.0; /* halve last step if it is greater than what remains divided this factor */
1032: adapt->wnormtype = NORM_2;
1033: adapt->timestepjustdecreased_delay = 0;
1035: *inadapt = adapt;
1036: return 0;
1037: }