1 
2 #include <petsc/private/tsimpl.h> /*I  "petscts.h" I*/
3 
4 PetscClassId TSADAPT_CLASSID;
5 
6 static PetscFunctionList TSAdaptList;
7 static PetscBool         TSAdaptPackageInitialized;
8 static PetscBool         TSAdaptRegisterAllCalled;
9 
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);
16 
17 /*@C
18    TSAdaptRegister -  adds a TSAdapt implementation
19 
20    Not Collective
21 
22    Input Parameters:
23 +  name_scheme - name of user-defined adaptivity scheme
24 -  routine_create - routine to create method context
25 
26    Notes:
27    TSAdaptRegister() may be called multiple times to add several user-defined families.
28 
29    Sample usage:
30 .vb
31    TSAdaptRegister("my_scheme",MySchemeCreate);
32 .ve
33 
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
38 
39    Level: advanced
40 
41 .seealso: TSAdaptRegisterAll()
42 @*/
TSAdaptRegister(const char sname[],PetscErrorCode (* function)(TSAdapt))43 PetscErrorCode  TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
44 {
45   PetscErrorCode ierr;
46 
47   PetscFunctionBegin;
48   ierr = TSAdaptInitializePackage();CHKERRQ(ierr);
49   ierr = PetscFunctionListAdd(&TSAdaptList,sname,function);CHKERRQ(ierr);
50   PetscFunctionReturn(0);
51 }
52 
53 /*@C
54   TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
55 
56   Not Collective
57 
58   Level: advanced
59 
60 .seealso: TSAdaptRegisterDestroy()
61 @*/
TSAdaptRegisterAll(void)62 PetscErrorCode  TSAdaptRegisterAll(void)
63 {
64   PetscErrorCode ierr;
65 
66   PetscFunctionBegin;
67   if (TSAdaptRegisterAllCalled) PetscFunctionReturn(0);
68   TSAdaptRegisterAllCalled = PETSC_TRUE;
69   ierr = TSAdaptRegister(TSADAPTNONE,   TSAdaptCreate_None);CHKERRQ(ierr);
70   ierr = TSAdaptRegister(TSADAPTBASIC,  TSAdaptCreate_Basic);CHKERRQ(ierr);
71   ierr = TSAdaptRegister(TSADAPTDSP,    TSAdaptCreate_DSP);CHKERRQ(ierr);
72   ierr = TSAdaptRegister(TSADAPTCFL,    TSAdaptCreate_CFL);CHKERRQ(ierr);
73   ierr = TSAdaptRegister(TSADAPTGLEE,   TSAdaptCreate_GLEE);CHKERRQ(ierr);
74   ierr = TSAdaptRegister(TSADAPTHISTORY,TSAdaptCreate_History);CHKERRQ(ierr);
75   PetscFunctionReturn(0);
76 }
77 
78 /*@C
79   TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
80   called from PetscFinalize().
81 
82   Level: developer
83 
84 .seealso: PetscFinalize()
85 @*/
TSAdaptFinalizePackage(void)86 PetscErrorCode  TSAdaptFinalizePackage(void)
87 {
88   PetscErrorCode ierr;
89 
90   PetscFunctionBegin;
91   ierr = PetscFunctionListDestroy(&TSAdaptList);CHKERRQ(ierr);
92   TSAdaptPackageInitialized = PETSC_FALSE;
93   TSAdaptRegisterAllCalled  = PETSC_FALSE;
94   PetscFunctionReturn(0);
95 }
96 
97 /*@C
98   TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
99   called from TSInitializePackage().
100 
101   Level: developer
102 
103 .seealso: PetscInitialize()
104 @*/
TSAdaptInitializePackage(void)105 PetscErrorCode  TSAdaptInitializePackage(void)
106 {
107   PetscErrorCode ierr;
108 
109   PetscFunctionBegin;
110   if (TSAdaptPackageInitialized) PetscFunctionReturn(0);
111   TSAdaptPackageInitialized = PETSC_TRUE;
112   ierr = PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);CHKERRQ(ierr);
113   ierr = TSAdaptRegisterAll();CHKERRQ(ierr);
114   ierr = PetscRegisterFinalize(TSAdaptFinalizePackage);CHKERRQ(ierr);
115   PetscFunctionReturn(0);
116 }
117 
118 /*@C
119   TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE
120 
121   Logicially Collective on TSAdapt
122 
123   Input Parameter:
124 + adapt - the TS adapter, most likely obtained with TSGetAdapt()
125 - type - either  TSADAPTBASIC or TSADAPTNONE
126 
127   Options Database:
128 . -ts_adapt_type <basic or dsp or none> - to set the adapter type
129 
130   Level: intermediate
131 
132 .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
133 @*/
TSAdaptSetType(TSAdapt adapt,TSAdaptType type)134 PetscErrorCode  TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
135 {
136   PetscBool      match;
137   PetscErrorCode ierr,(*r)(TSAdapt);
138 
139   PetscFunctionBegin;
140   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
141   PetscValidCharPointer(type,2);
142   ierr = PetscObjectTypeCompare((PetscObject)adapt,type,&match);CHKERRQ(ierr);
143   if (match) PetscFunctionReturn(0);
144   ierr = PetscFunctionListFind(TSAdaptList,type,&r);CHKERRQ(ierr);
145   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
146   if (adapt->ops->destroy) {ierr = (*adapt->ops->destroy)(adapt);CHKERRQ(ierr);}
147   ierr = PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));CHKERRQ(ierr);
148   ierr = PetscObjectChangeTypeName((PetscObject)adapt,type);CHKERRQ(ierr);
149   ierr = (*r)(adapt);CHKERRQ(ierr);
150   PetscFunctionReturn(0);
151 }
152 
153 /*@C
154   TSAdaptGetType - gets the TS adapter method type (as a string).
155 
156   Not Collective
157 
158   Input Parameter:
159 . adapt - The TS adapter, most likely obtained with TSGetAdapt()
160 
161   Output Parameter:
162 . type - The name of TS adapter method
163 
164   Level: intermediate
165 
166 .seealso TSAdaptSetType()
167 @*/
TSAdaptGetType(TSAdapt adapt,TSAdaptType * type)168 PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
169 {
170   PetscFunctionBegin;
171   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
172   PetscValidPointer(type,2);
173   *type = ((PetscObject)adapt)->type_name;
174   PetscFunctionReturn(0);
175 }
176 
TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])177 PetscErrorCode  TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
178 {
179   PetscErrorCode ierr;
180 
181   PetscFunctionBegin;
182   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
183   ierr = PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);CHKERRQ(ierr);
184   PetscFunctionReturn(0);
185 }
186 
187 /*@C
188   TSAdaptLoad - Loads a TSAdapt that has been stored in binary  with TSAdaptView().
189 
190   Collective on PetscViewer
191 
192   Input Parameters:
193 + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
194            some related function before a call to TSAdaptLoad().
195 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
196            HDF5 file viewer, obtained from PetscViewerHDF5Open()
197 
198    Level: intermediate
199 
200   Notes:
201    The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
202 
203   Notes for advanced users:
204   Most users should not need to know the details of the binary storage
205   format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
206   But for anyone who's interested, the standard binary matrix storage
207   format is
208 .vb
209      has not yet been determined
210 .ve
211 
212 .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
213 @*/
TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)214 PetscErrorCode  TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
215 {
216   PetscErrorCode ierr;
217   PetscBool      isbinary;
218   char           type[256];
219 
220   PetscFunctionBegin;
221   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
222   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
223   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
224   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
225 
226   ierr = PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);CHKERRQ(ierr);
227   ierr = TSAdaptSetType(adapt,type);CHKERRQ(ierr);
228   if (adapt->ops->load) {
229     ierr = (*adapt->ops->load)(adapt,viewer);CHKERRQ(ierr);
230   }
231   PetscFunctionReturn(0);
232 }
233 
TSAdaptView(TSAdapt adapt,PetscViewer viewer)234 PetscErrorCode  TSAdaptView(TSAdapt adapt,PetscViewer viewer)
235 {
236   PetscErrorCode ierr;
237   PetscBool      iascii,isbinary,isnone,isglee;
238 
239   PetscFunctionBegin;
240   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
241   if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);CHKERRQ(ierr);}
242   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
243   PetscCheckSameComm(adapt,1,viewer,2);
244   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
245   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
246   if (iascii) {
247     ierr = PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);CHKERRQ(ierr);
248     ierr = PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);CHKERRQ(ierr);
249     ierr = PetscObjectTypeCompare((PetscObject)adapt,TSADAPTGLEE,&isglee);CHKERRQ(ierr);
250     if (!isnone) {
251       if (adapt->always_accept) {ierr = PetscViewerASCIIPrintf(viewer,"  always accepting steps\n");CHKERRQ(ierr);}
252       ierr = PetscViewerASCIIPrintf(viewer,"  safety factor %g\n",(double)adapt->safety);CHKERRQ(ierr);
253       ierr = PetscViewerASCIIPrintf(viewer,"  extra safety factor after step rejection %g\n",(double)adapt->reject_safety);CHKERRQ(ierr);
254       ierr = PetscViewerASCIIPrintf(viewer,"  clip fastest increase %g\n",(double)adapt->clip[1]);CHKERRQ(ierr);
255       ierr = PetscViewerASCIIPrintf(viewer,"  clip fastest decrease %g\n",(double)adapt->clip[0]);CHKERRQ(ierr);
256       ierr = PetscViewerASCIIPrintf(viewer,"  maximum allowed timestep %g\n",(double)adapt->dt_max);CHKERRQ(ierr);
257       ierr = PetscViewerASCIIPrintf(viewer,"  minimum allowed timestep %g\n",(double)adapt->dt_min);CHKERRQ(ierr);
258       ierr = PetscViewerASCIIPrintf(viewer,"  maximum solution absolute value to be ignored %g\n",(double)adapt->ignore_max);CHKERRQ(ierr);
259     }
260     if (isglee) {
261       if (adapt->glee_use_local) {
262         ierr = PetscViewerASCIIPrintf(viewer,"  GLEE uses local error control\n");CHKERRQ(ierr);
263       } else {
264         ierr = PetscViewerASCIIPrintf(viewer,"  GLEE uses global error control\n");CHKERRQ(ierr);
265       }
266     }
267     if (adapt->ops->view) {
268       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
269       ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr);
270       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
271     }
272   } else if (isbinary) {
273     char type[256];
274 
275     /* need to save FILE_CLASS_ID for adapt class */
276     ierr = PetscStrncpy(type,((PetscObject)adapt)->type_name,256);CHKERRQ(ierr);
277     ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr);
278   } else if (adapt->ops->view) {
279     ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr);
280   }
281   PetscFunctionReturn(0);
282 }
283 
284 /*@
285    TSAdaptReset - Resets a TSAdapt context.
286 
287    Collective on TS
288 
289    Input Parameter:
290 .  adapt - the TSAdapt context obtained from TSAdaptCreate()
291 
292    Level: developer
293 
294 .seealso: TSAdaptCreate(), TSAdaptDestroy()
295 @*/
TSAdaptReset(TSAdapt adapt)296 PetscErrorCode  TSAdaptReset(TSAdapt adapt)
297 {
298   PetscErrorCode ierr;
299 
300   PetscFunctionBegin;
301   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
302   if (adapt->ops->reset) {ierr = (*adapt->ops->reset)(adapt);CHKERRQ(ierr);}
303   PetscFunctionReturn(0);
304 }
305 
TSAdaptDestroy(TSAdapt * adapt)306 PetscErrorCode  TSAdaptDestroy(TSAdapt *adapt)
307 {
308   PetscErrorCode ierr;
309 
310   PetscFunctionBegin;
311   if (!*adapt) PetscFunctionReturn(0);
312   PetscValidHeaderSpecific(*adapt,TSADAPT_CLASSID,1);
313   if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; PetscFunctionReturn(0);}
314 
315   ierr = TSAdaptReset(*adapt);CHKERRQ(ierr);
316 
317   if ((*adapt)->ops->destroy) {ierr = (*(*adapt)->ops->destroy)(*adapt);CHKERRQ(ierr);}
318   ierr = PetscViewerDestroy(&(*adapt)->monitor);CHKERRQ(ierr);
319   ierr = PetscHeaderDestroy(adapt);CHKERRQ(ierr);
320   PetscFunctionReturn(0);
321 }
322 
323 /*@
324    TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
325 
326    Collective on TSAdapt
327 
328    Input Arguments:
329 +  adapt - adaptive controller context
330 -  flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
331 
332    Options Database Keys:
333 .  -ts_adapt_monitor - to turn on monitoring
334 
335    Level: intermediate
336 
337 .seealso: TSAdaptChoose()
338 @*/
TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)339 PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
340 {
341   PetscErrorCode ierr;
342 
343   PetscFunctionBegin;
344   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
345   PetscValidLogicalCollectiveBool(adapt,flg,2);
346   if (flg) {
347     if (!adapt->monitor) {ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);CHKERRQ(ierr);}
348   } else {
349     ierr = PetscViewerDestroy(&adapt->monitor);CHKERRQ(ierr);
350   }
351   PetscFunctionReturn(0);
352 }
353 
354 /*@C
355    TSAdaptSetCheckStage - Set a callback to check convergence for a stage
356 
357    Logically collective on TSAdapt
358 
359    Input Arguments:
360 +  adapt - adaptive controller context
361 -  func - stage check function
362 
363    Arguments of func:
364 $  PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
365 
366 +  adapt - adaptive controller context
367 .  ts - time stepping context
368 -  accept - pending choice of whether to accept, can be modified by this routine
369 
370    Level: advanced
371 
372 .seealso: TSAdaptChoose()
373 @*/
TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (* func)(TSAdapt,TS,PetscReal,Vec,PetscBool *))374 PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*))
375 {
376 
377   PetscFunctionBegin;
378   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
379   adapt->checkstage = func;
380   PetscFunctionReturn(0);
381 }
382 
383 /*@
384    TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of
385    any error or stability condition not meeting the prescribed goal.
386 
387    Logically collective on TSAdapt
388 
389    Input Arguments:
390 +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
391 -  flag - whether to always accept steps
392 
393    Options Database Keys:
394 .  -ts_adapt_always_accept - to always accept steps
395 
396    Level: intermediate
397 
398 .seealso: TSAdapt, TSAdaptChoose()
399 @*/
TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)400 PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)
401 {
402   PetscFunctionBegin;
403   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
404   PetscValidLogicalCollectiveBool(adapt,flag,2);
405   adapt->always_accept = flag;
406   PetscFunctionReturn(0);
407 }
408 
409 /*@
410    TSAdaptSetSafety - Set safety factors
411 
412    Logically collective on TSAdapt
413 
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
418 
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
422 
423    Level: intermediate
424 
425 .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose()
426 @*/
TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)427 PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)
428 {
429   PetscFunctionBegin;
430   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
431   PetscValidLogicalCollectiveReal(adapt,safety,2);
432   PetscValidLogicalCollectiveReal(adapt,reject_safety,3);
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   PetscFunctionReturn(0);
440 }
441 
442 /*@
443    TSAdaptGetSafety - Get safety factors
444 
445    Not Collective
446 
447    Input Arguments:
448 .  adapt - adaptive controller context
449 
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
453 
454    Level: intermediate
455 
456 .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose()
457 @*/
TSAdaptGetSafety(TSAdapt adapt,PetscReal * safety,PetscReal * reject_safety)458 PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety)
459 {
460   PetscFunctionBegin;
461   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
462   if (safety)        PetscValidRealPointer(safety,2);
463   if (reject_safety) PetscValidRealPointer(reject_safety,3);
464   if (safety)        *safety        = adapt->safety;
465   if (reject_safety) *reject_safety = adapt->reject_safety;
466   PetscFunctionReturn(0);
467 }
468 
469 /*@
470    TSAdaptSetMaxIgnore - Set error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value). A negative value (default) of the threshold leads to considering all solution components.
471 
472    Logically collective on TSAdapt
473 
474    Input Arguments:
475 +  adapt - adaptive controller context
476 -  max_ignore - threshold for solution components that are ignored during error estimation
477 
478    Options Database Keys:
479 .  -ts_adapt_max_ignore <max_ignore> - to set the threshold
480 
481    Level: intermediate
482 
483 .seealso: TSAdapt, TSAdaptGetMaxIgnore(), TSAdaptChoose()
484 @*/
TSAdaptSetMaxIgnore(TSAdapt adapt,PetscReal max_ignore)485 PetscErrorCode TSAdaptSetMaxIgnore(TSAdapt adapt,PetscReal max_ignore)
486 {
487   PetscFunctionBegin;
488   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
489   PetscValidLogicalCollectiveReal(adapt,max_ignore,2);
490   adapt->ignore_max = max_ignore;
491   PetscFunctionReturn(0);
492 }
493 
494 /*@
495    TSAdaptGetMaxIgnore - Get error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value).
496 
497    Not Collective
498 
499    Input Arguments:
500 .  adapt - adaptive controller context
501 
502    Ouput Arguments:
503 .  max_ignore - threshold for solution components that are ignored during error estimation
504 
505    Level: intermediate
506 
507 .seealso: TSAdapt, TSAdaptSetMaxIgnore(), TSAdaptChoose()
508 @*/
TSAdaptGetMaxIgnore(TSAdapt adapt,PetscReal * max_ignore)509 PetscErrorCode TSAdaptGetMaxIgnore(TSAdapt adapt,PetscReal *max_ignore)
510 {
511   PetscFunctionBegin;
512   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
513   PetscValidRealPointer(max_ignore,2);
514   *max_ignore = adapt->ignore_max;
515   PetscFunctionReturn(0);
516 }
517 
518 
519 /*@
520    TSAdaptSetClip - Sets the admissible decrease/increase factor in step size
521 
522    Logically collective on TSAdapt
523 
524    Input Arguments:
525 +  adapt - adaptive controller context
526 .  low - admissible decrease factor
527 -  high - admissible increase factor
528 
529    Options Database Keys:
530 .  -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors
531 
532    Level: intermediate
533 
534 .seealso: TSAdaptChoose(), TSAdaptGetClip(), TSAdaptSetScaleSolveFailed()
535 @*/
TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)536 PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
537 {
538   PetscFunctionBegin;
539   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
540   PetscValidLogicalCollectiveReal(adapt,low,2);
541   PetscValidLogicalCollectiveReal(adapt,high,3);
542   if (low  != PETSC_DEFAULT && low  < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
543   if (low  != PETSC_DEFAULT && low  > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
544   if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be greater than one",(double)high);
545   if (low  != PETSC_DEFAULT) adapt->clip[0] = low;
546   if (high != PETSC_DEFAULT) adapt->clip[1] = high;
547   PetscFunctionReturn(0);
548 }
549 
550 /*@
551    TSAdaptGetClip - Gets the admissible decrease/increase factor in step size
552 
553    Not Collective
554 
555    Input Arguments:
556 .  adapt - adaptive controller context
557 
558    Ouput Arguments:
559 +  low - optional, admissible decrease factor
560 -  high - optional, admissible increase factor
561 
562    Level: intermediate
563 
564 .seealso: TSAdaptChoose(), TSAdaptSetClip(), TSAdaptSetScaleSolveFailed()
565 @*/
TSAdaptGetClip(TSAdapt adapt,PetscReal * low,PetscReal * high)566 PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
567 {
568   PetscFunctionBegin;
569   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
570   if (low)  PetscValidRealPointer(low,2);
571   if (high) PetscValidRealPointer(high,3);
572   if (low)  *low  = adapt->clip[0];
573   if (high) *high = adapt->clip[1];
574   PetscFunctionReturn(0);
575 }
576 
577 /*@
578    TSAdaptSetScaleSolveFailed - Scale step by this factor if solve fails
579 
580    Logically collective on TSAdapt
581 
582    Input Arguments:
583 +  adapt - adaptive controller context
584 -  scale - scale
585 
586    Options Database Keys:
587 .  -ts_adapt_scale_solve_failed <scale> - to set scale step by this factor if solve fails
588 
589    Level: intermediate
590 
591 .seealso: TSAdaptChoose(), TSAdaptGetScaleSolveFailed(), TSAdaptGetClip()
592 @*/
TSAdaptSetScaleSolveFailed(TSAdapt adapt,PetscReal scale)593 PetscErrorCode TSAdaptSetScaleSolveFailed(TSAdapt adapt,PetscReal scale)
594 {
595   PetscFunctionBegin;
596   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
597   PetscValidLogicalCollectiveReal(adapt,scale,2);
598   if (scale != PETSC_DEFAULT && scale <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be positive",(double)scale);
599   if (scale != PETSC_DEFAULT && scale  > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be less than one",(double)scale);
600   if (scale != PETSC_DEFAULT) adapt->scale_solve_failed = scale;
601   PetscFunctionReturn(0);
602 }
603 
604 /*@
605    TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size
606 
607    Not Collective
608 
609    Input Arguments:
610 .  adapt - adaptive controller context
611 
612    Ouput Arguments:
613 .  scale - scale factor
614 
615    Level: intermediate
616 
617 .seealso: TSAdaptChoose(), TSAdaptSetScaleSolveFailed(), TSAdaptSetClip()
618 @*/
TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal * scale)619 PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal *scale)
620 {
621   PetscFunctionBegin;
622   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
623   if (scale)  PetscValidRealPointer(scale,2);
624   if (scale)  *scale  = adapt->scale_solve_failed;
625   PetscFunctionReturn(0);
626 }
627 
628 /*@
629    TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller
630 
631    Logically collective on TSAdapt
632 
633    Input Arguments:
634 +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
635 .  hmin - minimum time step
636 -  hmax - maximum time step
637 
638    Options Database Keys:
639 +  -ts_adapt_dt_min <min> - to set minimum time step
640 -  -ts_adapt_dt_max <max> - to set maximum time step
641 
642    Level: intermediate
643 
644 .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
645 @*/
TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)646 PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
647 {
648 
649   PetscFunctionBegin;
650   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
651   PetscValidLogicalCollectiveReal(adapt,hmin,2);
652   PetscValidLogicalCollectiveReal(adapt,hmax,3);
653   if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
654   if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
655   if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
656   if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
657   hmin = adapt->dt_min;
658   hmax = adapt->dt_max;
659   if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must greater than minimum time step %g",(double)hmax,(double)hmin);
660   PetscFunctionReturn(0);
661 }
662 
663 /*@
664    TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller
665 
666    Not Collective
667 
668    Input Arguments:
669 .  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
670 
671    Output Arguments:
672 +  hmin - minimum time step
673 -  hmax - maximum time step
674 
675    Level: intermediate
676 
677 .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
678 @*/
TSAdaptGetStepLimits(TSAdapt adapt,PetscReal * hmin,PetscReal * hmax)679 PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
680 {
681 
682   PetscFunctionBegin;
683   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
684   if (hmin) PetscValidRealPointer(hmin,2);
685   if (hmax) PetscValidRealPointer(hmax,3);
686   if (hmin) *hmin = adapt->dt_min;
687   if (hmax) *hmax = adapt->dt_max;
688   PetscFunctionReturn(0);
689 }
690 
691 /*
692    TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
693 
694    Collective on TSAdapt
695 
696    Input Parameter:
697 .  adapt - the TSAdapt context
698 
699    Options Database Keys:
700 +  -ts_adapt_type <type> - algorithm to use for adaptivity
701 .  -ts_adapt_always_accept - always accept steps regardless of error/stability goals
702 .  -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
703 .  -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
704 .  -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
705 .  -ts_adapt_dt_min <min> - minimum timestep to use
706 .  -ts_adapt_dt_max <max> - maximum timestep to use
707 .  -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
708 .  -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
709 -  -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver
710 
711    Level: advanced
712 
713    Notes:
714    This function is automatically called by TSSetFromOptions()
715 
716 .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
717           TSAdaptSetClip(), TSAdaptSetScaleSolveFailed(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
718 */
TSAdaptSetFromOptions(PetscOptionItems * PetscOptionsObject,TSAdapt adapt)719 PetscErrorCode  TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
720 {
721   PetscErrorCode ierr;
722   char           type[256] = TSADAPTBASIC;
723   PetscReal      safety,reject_safety,clip[2],scale,hmin,hmax;
724   PetscBool      set,flg;
725   PetscInt       two;
726 
727   PetscFunctionBegin;
728   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
729   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
730    * function can only be called from inside TSSetFromOptions()  */
731   ierr = PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");CHKERRQ(ierr);
732   ierr = PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);CHKERRQ(ierr);
733   if (flg || !((PetscObject)adapt)->type_name) {
734     ierr = TSAdaptSetType(adapt,type);CHKERRQ(ierr);
735   }
736 
737   ierr = PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);CHKERRQ(ierr);
738   if (set) {ierr = TSAdaptSetAlwaysAccept(adapt,flg);CHKERRQ(ierr);}
739 
740   safety = adapt->safety; reject_safety = adapt->reject_safety;
741   ierr = PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);CHKERRQ(ierr);
742   ierr = PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);CHKERRQ(ierr);
743   if (set || flg) {ierr = TSAdaptSetSafety(adapt,safety,reject_safety);CHKERRQ(ierr);}
744 
745   two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
746   ierr = PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);CHKERRQ(ierr);
747   if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
748   if (set) {ierr = TSAdaptSetClip(adapt,clip[0],clip[1]);CHKERRQ(ierr);}
749 
750   hmin = adapt->dt_min; hmax = adapt->dt_max;
751   ierr = PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);CHKERRQ(ierr);
752   ierr = PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);CHKERRQ(ierr);
753   if (set || flg) {ierr = TSAdaptSetStepLimits(adapt,hmin,hmax);CHKERRQ(ierr);}
754 
755   ierr = PetscOptionsReal("-ts_adapt_max_ignore","Adaptor ignores (absolute) solution values smaller than this value","",adapt->ignore_max,&adapt->ignore_max,&set);CHKERRQ(ierr);
756   ierr = PetscOptionsBool("-ts_adapt_glee_use_local","GLEE adaptor uses local error estimation for step control","",adapt->glee_use_local,&adapt->glee_use_local,&set);CHKERRQ(ierr);
757 
758   ierr = PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","TSAdaptSetScaleSolveFailed",adapt->scale_solve_failed,&scale,&set);CHKERRQ(ierr);
759   if (set) {ierr = TSAdaptSetScaleSolveFailed(adapt,scale);CHKERRQ(ierr);}
760 
761   ierr = PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);CHKERRQ(ierr);
762   if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");
763 
764   ierr = 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);CHKERRQ(ierr);
765 
766   ierr = PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr);
767   if (set) {ierr = TSAdaptSetMonitor(adapt,flg);CHKERRQ(ierr);}
768 
769   if (adapt->ops->setfromoptions) {ierr = (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);CHKERRQ(ierr);}
770   ierr = PetscOptionsTail();CHKERRQ(ierr);
771   PetscFunctionReturn(0);
772 }
773 
774 /*@
775    TSAdaptCandidatesClear - clear any previously set candidate schemes
776 
777    Logically collective on TSAdapt
778 
779    Input Argument:
780 .  adapt - adaptive controller
781 
782    Level: developer
783 
784 .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
785 @*/
TSAdaptCandidatesClear(TSAdapt adapt)786 PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
787 {
788   PetscErrorCode ierr;
789 
790   PetscFunctionBegin;
791   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
792   ierr = PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));CHKERRQ(ierr);
793   PetscFunctionReturn(0);
794 }
795 
796 /*@C
797    TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
798 
799    Logically collective on TSAdapt
800 
801    Input Arguments:
802 +  adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
803 .  name - name of the candidate scheme to add
804 .  order - order of the candidate scheme
805 .  stageorder - stage order of the candidate scheme
806 .  ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
807 .  cost - relative measure of the amount of work required for the candidate scheme
808 -  inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
809 
810    Note:
811    This routine is not available in Fortran.
812 
813    Level: developer
814 
815 .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
816 @*/
TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)817 PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
818 {
819   PetscInt c;
820 
821   PetscFunctionBegin;
822   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
823   if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
824   if (inuse) {
825     if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
826     adapt->candidates.inuse_set = PETSC_TRUE;
827   }
828   /* first slot if this is the current scheme, otherwise the next available slot */
829   c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
830 
831   adapt->candidates.name[c]       = name;
832   adapt->candidates.order[c]      = order;
833   adapt->candidates.stageorder[c] = stageorder;
834   adapt->candidates.ccfl[c]       = ccfl;
835   adapt->candidates.cost[c]       = cost;
836   adapt->candidates.n++;
837   PetscFunctionReturn(0);
838 }
839 
840 /*@C
841    TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
842 
843    Not Collective
844 
845    Input Arguments:
846 .  adapt - time step adaptivity context
847 
848    Output Arguments:
849 +  n - number of candidate schemes, always at least 1
850 .  order - the order of each candidate scheme
851 .  stageorder - the stage order of each candidate scheme
852 .  ccfl - the CFL coefficient of each scheme
853 -  cost - the relative cost of each scheme
854 
855    Level: developer
856 
857    Note:
858    The current scheme is always returned in the first slot
859 
860 .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
861 @*/
TSAdaptCandidatesGet(TSAdapt adapt,PetscInt * n,const PetscInt ** order,const PetscInt ** stageorder,const PetscReal ** ccfl,const PetscReal ** cost)862 PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
863 {
864   PetscFunctionBegin;
865   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
866   if (n) *n = adapt->candidates.n;
867   if (order) *order = adapt->candidates.order;
868   if (stageorder) *stageorder = adapt->candidates.stageorder;
869   if (ccfl) *ccfl = adapt->candidates.ccfl;
870   if (cost) *cost = adapt->candidates.cost;
871   PetscFunctionReturn(0);
872 }
873 
874 /*@C
875    TSAdaptChoose - choose which method and step size to use for the next step
876 
877    Collective on TSAdapt
878 
879    Input Arguments:
880 +  adapt - adaptive contoller
881 -  h - current step size
882 
883    Output Arguments:
884 +  next_sc - optional, scheme to use for the next step
885 .  next_h - step size to use for the next step
886 -  accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
887 
888    Note:
889    The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
890    being retried after an initial rejection.
891 
892    Level: developer
893 
894 .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
895 @*/
TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept)896 PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
897 {
898   PetscErrorCode ierr;
899   PetscInt       ncandidates = adapt->candidates.n;
900   PetscInt       scheme = 0;
901   PetscReal      wlte = -1.0;
902   PetscReal      wltea = -1.0;
903   PetscReal      wlter = -1.0;
904 
905   PetscFunctionBegin;
906   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
907   PetscValidHeaderSpecific(ts,TS_CLASSID,2);
908   if (next_sc) PetscValidIntPointer(next_sc,4);
909   PetscValidPointer(next_h,5);
910   PetscValidIntPointer(accept,6);
911   if (next_sc) *next_sc = 0;
912 
913   /* Do not mess with adaptivity while handling events*/
914   if (ts->event && ts->event->status != TSEVENT_NONE) {
915     *next_h = h;
916     *accept = PETSC_TRUE;
917     PetscFunctionReturn(0);
918   }
919 
920   ierr = (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);CHKERRQ(ierr);
921   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);
922   if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
923   if (next_sc) *next_sc = scheme;
924 
925   if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
926     /* Increase/reduce step size if end time of next step is close to or overshoots max time */
927     PetscReal t = ts->ptime + ts->time_step, h = *next_h;
928     PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
929     PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
930     PetscReal b = adapt->matchstepfac[1];
931     if (t < tmax && tend > tmax) *next_h = hmax;
932     if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
933     if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
934   }
935 
936   if (adapt->monitor) {
937     const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
938     ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
939     if (wlte < 0) {
940       ierr = 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);CHKERRQ(ierr);
941     } else {
942       ierr = 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);CHKERRQ(ierr);
943     }
944     ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
945   }
946   PetscFunctionReturn(0);
947 }
948 
949 /*@
950    TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
951                                      before increasing the time step.
952 
953    Logicially Collective on TSAdapt
954 
955    Input Arguments:
956 +  adapt - adaptive controller context
957 -  cnt - the number of timesteps
958 
959    Options Database Key:
960 .  -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase
961 
962    Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
963           The successful use of this option is problem dependent
964 
965    Developer Note: there is no theory to support this option
966 
967    Level: advanced
968 
969 .seealso:
970 @*/
TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)971 PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
972 {
973   PetscFunctionBegin;
974   adapt->timestepjustdecreased_delay = cnt;
975   PetscFunctionReturn(0);
976 }
977 
978 /*@
979    TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible)
980 
981    Collective on TSAdapt
982 
983    Input Arguments:
984 +  adapt - adaptive controller context
985 .  ts - time stepper
986 .  t - Current simulation time
987 -  Y - Current solution vector
988 
989    Output Arguments:
990 .  accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
991 
992    Level: developer
993 
994 .seealso:
995 @*/
TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool * accept)996 PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
997 {
998   PetscErrorCode      ierr;
999   SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
1000 
1001   PetscFunctionBegin;
1002   PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1);
1003   PetscValidHeaderSpecific(ts,TS_CLASSID,2);
1004   PetscValidIntPointer(accept,3);
1005 
1006   if (ts->snes) {ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr);}
1007   if (snesreason < 0) {
1008     *accept = PETSC_FALSE;
1009     if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
1010       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
1011       ierr = PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);CHKERRQ(ierr);
1012       if (adapt->monitor) {
1013         ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
1014         ierr = 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);CHKERRQ(ierr);
1015         ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
1016       }
1017     }
1018   } else {
1019     *accept = PETSC_TRUE;
1020     ierr = TSFunctionDomainError(ts,t,Y,accept);CHKERRQ(ierr);
1021     if (*accept && adapt->checkstage) {
1022       ierr = (*adapt->checkstage)(adapt,ts,t,Y,accept);CHKERRQ(ierr);
1023       if (!*accept) {
1024         ierr = PetscInfo1(ts,"Step=%D, solution rejected by user function provided by TSSetFunctionDomainError()\n",ts->steps);CHKERRQ(ierr);
1025         if (adapt->monitor) {
1026           ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
1027           ierr = PetscViewerASCIIPrintf(adapt->monitor,"    TSAdapt %s step %3D stage rejected by user function provided by TSSetFunctionDomainError()\n",((PetscObject)adapt)->type_name,ts->steps);CHKERRQ(ierr);
1028           ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
1029         }
1030       }
1031     }
1032   }
1033 
1034   if (!(*accept) && !ts->reason) {
1035     PetscReal dt,new_dt;
1036     ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
1037     new_dt = dt * adapt->scale_solve_failed;
1038     ierr = TSSetTimeStep(ts,new_dt);CHKERRQ(ierr);
1039     adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
1040     if (adapt->monitor) {
1041       ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
1042       ierr = 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);CHKERRQ(ierr);
1043       ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr);
1044     }
1045   }
1046   PetscFunctionReturn(0);
1047 }
1048 
1049 /*@
1050   TSAdaptCreate - create an adaptive controller context for time stepping
1051 
1052   Collective
1053 
1054   Input Parameter:
1055 . comm - The communicator
1056 
1057   Output Parameter:
1058 . adapt - new TSAdapt object
1059 
1060   Level: developer
1061 
1062   Notes:
1063   TSAdapt creation is handled by TS, so users should not need to call this function.
1064 
1065 .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1066 @*/
TSAdaptCreate(MPI_Comm comm,TSAdapt * inadapt)1067 PetscErrorCode  TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1068 {
1069   PetscErrorCode ierr;
1070   TSAdapt        adapt;
1071 
1072   PetscFunctionBegin;
1073   PetscValidPointer(inadapt,1);
1074   *inadapt = NULL;
1075   ierr = TSAdaptInitializePackage();CHKERRQ(ierr);
1076 
1077   ierr = PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);CHKERRQ(ierr);
1078 
1079   adapt->always_accept      = PETSC_FALSE;
1080   adapt->safety             = 0.9;
1081   adapt->reject_safety      = 0.5;
1082   adapt->clip[0]            = 0.1;
1083   adapt->clip[1]            = 10.;
1084   adapt->dt_min             = 1e-20;
1085   adapt->dt_max             = 1e+20;
1086   adapt->ignore_max         = -1.0;
1087   adapt->glee_use_local     = PETSC_TRUE;
1088   adapt->scale_solve_failed = 0.25;
1089   /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1090      to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1091   adapt->matchstepfac[0]    = 0.01; /* allow 1% step size increase in the last step */
1092   adapt->matchstepfac[1]    = 2.0;  /* halve last step if it is greater than what remains divided this factor */
1093   adapt->wnormtype          = NORM_2;
1094   adapt->timestepjustdecreased_delay = 0;
1095 
1096   *inadapt = adapt;
1097   PetscFunctionReturn(0);
1098 }
1099