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