Actual source code: tsadapt.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

  2:  #include <petsc/private/tsimpl.h>

  4: PetscClassId TSADAPT_CLASSID;

  6: static PetscFunctionList TSAdaptList;
  7: static PetscBool         TSAdaptPackageInitialized;
  8: static PetscBool         TSAdaptRegisterAllCalled;

 10: PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt);
 11: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt);
 12: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt);
 13: PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt);
 14: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt);

 16: /*@C
 17:    TSAdaptRegister -  adds a TSAdapt implementation

 19:    Not Collective

 21:    Input Parameters:
 22: +  name_scheme - name of user-defined adaptivity scheme
 23: -  routine_create - routine to create method context

 25:    Notes:
 26:    TSAdaptRegister() may be called multiple times to add several user-defined families.

 28:    Sample usage:
 29: .vb
 30:    TSAdaptRegister("my_scheme",MySchemeCreate);
 31: .ve

 33:    Then, your scheme can be chosen with the procedural interface via
 34: $     TSAdaptSetType(ts,"my_scheme")
 35:    or at runtime via the option
 36: $     -ts_adapt_type my_scheme

 38:    Level: advanced

 40: .keywords: TSAdapt, register

 42: .seealso: TSAdaptRegisterAll()
 43: @*/
 44: PetscErrorCode  TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
 45: {

 49:   TSAdaptInitializePackage();
 50:   PetscFunctionListAdd(&TSAdaptList,sname,function);
 51:   return(0);
 52: }

 54: /*@C
 55:   TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt

 57:   Not Collective

 59:   Level: advanced

 61: .keywords: TSAdapt, register, all

 63: .seealso: TSAdaptRegisterDestroy()
 64: @*/
 65: PetscErrorCode  TSAdaptRegisterAll(void)
 66: {

 70:   if (TSAdaptRegisterAllCalled) return(0);
 71:   TSAdaptRegisterAllCalled = PETSC_TRUE;
 72:   TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
 73:   TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);
 74:   TSAdaptRegister(TSADAPTDSP,  TSAdaptCreate_DSP);
 75:   TSAdaptRegister(TSADAPTCFL,  TSAdaptCreate_CFL);
 76:   TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);
 77:   return(0);
 78: }

 80: /*@C
 81:   TSAdaptFinalizePackage - This function destroys everything in the TS package. It is
 82:   called from PetscFinalize().

 84:   Level: developer

 86: .keywords: Petsc, destroy, package
 87: .seealso: PetscFinalize()
 88: @*/
 89: PetscErrorCode  TSAdaptFinalizePackage(void)
 90: {

 94:   PetscFunctionListDestroy(&TSAdaptList);
 95:   TSAdaptPackageInitialized = PETSC_FALSE;
 96:   TSAdaptRegisterAllCalled  = PETSC_FALSE;
 97:   return(0);
 98: }

100: /*@C
101:   TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
102:   called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
103:   TSAdaptCreate() when using static libraries.

105:   Level: developer

107: .keywords: TSAdapt, initialize, package
108: .seealso: PetscInitialize()
109: @*/
110: PetscErrorCode  TSAdaptInitializePackage(void)
111: {

115:   if (TSAdaptPackageInitialized) return(0);
116:   TSAdaptPackageInitialized = PETSC_TRUE;
117:   PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
118:   TSAdaptRegisterAll();
119:   PetscRegisterFinalize(TSAdaptFinalizePackage);
120:   return(0);
121: }

123: /*@C
124:   TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE

126:   Logicially Collective on TSAdapt

128:   Input Parameter:
129: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
130: - type - either  TSADAPTBASIC or TSADAPTNONE

132:   Options Database:
133: . -ts_adapt_type <basic or dsp or none> - to set the adapter type

135:   Level: intermediate

137: .keywords: TSAdapt, create

139: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
140: @*/
141: PetscErrorCode  TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
142: {
143:   PetscBool      match;
144:   PetscErrorCode ierr,(*r)(TSAdapt);

149:   PetscObjectTypeCompare((PetscObject)adapt,type,&match);
150:   if (match) return(0);
151:   PetscFunctionListFind(TSAdaptList,type,&r);
152:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
153:   if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
154:   PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
155:   PetscObjectChangeTypeName((PetscObject)adapt,type);
156:   (*r)(adapt);
157:   return(0);
158: }

160: /*@C
161:   TSAdaptGetType - gets the TS adapter method type (as a string).

163:   Not Collective

165:   Input Parameter:
166: . adapt - The TS adapter, most likely obtained with TSGetAdapt()

168:   Output Parameter:
169: . type - The name of TS adapter method

171:   Level: intermediate

173: .keywords: TSAdapt, get, type
174: .seealso TSAdaptSetType()
175: @*/
176: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
177: {
181:   *type = ((PetscObject)adapt)->type_name;
182:   return(0);
183: }

185: PetscErrorCode  TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
186: {

191:   PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
192:   return(0);
193: }

195: /*@C
196:   TSAdaptLoad - Loads a TSAdapt that has been stored in binary  with TSAdaptView().

198:   Collective on PetscViewer

200:   Input Parameters:
201: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
202:            some related function before a call to TSAdaptLoad().
203: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
204:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

206:    Level: intermediate

208:   Notes:
209:    The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.

211:   Notes for advanced users:
212:   Most users should not need to know the details of the binary storage
213:   format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
214:   But for anyone who's interested, the standard binary matrix storage
215:   format is
216: .vb
217:      has not yet been determined
218: .ve

220: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
221: @*/
222: PetscErrorCode  TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
223: {
225:   PetscBool      isbinary;
226:   char           type[256];

231:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
232:   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

234:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
235:   TSAdaptSetType(adapt,type);
236:   if (adapt->ops->load) {
237:     (*adapt->ops->load)(adapt,viewer);
238:   }
239:   return(0);
240: }

242: PetscErrorCode  TSAdaptView(TSAdapt adapt,PetscViewer viewer)
243: {
245:   PetscBool      iascii,isbinary,isnone;

249:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
252:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
253:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
254:   if (iascii) {
255:     PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
256:     PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
257:     if (!isnone) {
258:       if (adapt->always_accept) {PetscViewerASCIIPrintf(viewer,"  always accepting steps\n");}
259:       PetscViewerASCIIPrintf(viewer,"  safety factor %g\n",(double)adapt->safety);
260:       PetscViewerASCIIPrintf(viewer,"  extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
261:       PetscViewerASCIIPrintf(viewer,"  clip fastest increase %g\n",(double)adapt->clip[1]);
262:       PetscViewerASCIIPrintf(viewer,"  clip fastest decrease %g\n",(double)adapt->clip[0]);
263:       PetscViewerASCIIPrintf(viewer,"  maximum allowed timestep %g\n",(double)adapt->dt_max);
264:       PetscViewerASCIIPrintf(viewer,"  minimum allowed timestep %g\n",(double)adapt->dt_min);
265:     }
266:     if (adapt->ops->view) {
267:       PetscViewerASCIIPushTab(viewer);
268:       (*adapt->ops->view)(adapt,viewer);
269:       PetscViewerASCIIPopTab(viewer);
270:     }
271:   } else if (isbinary) {
272:     char type[256];

274:     /* need to save FILE_CLASS_ID for adapt class */
275:     PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
276:     PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
277:   } else if (adapt->ops->view) {
278:     (*adapt->ops->view)(adapt,viewer);
279:   }
280:   return(0);
281: }

283: /*@
284:    TSAdaptReset - Resets a TSAdapt context.

286:    Collective on TS

288:    Input Parameter:
289: .  adapt - the TSAdapt context obtained from TSAdaptCreate()

291:    Level: developer

293: .seealso: TSAdaptCreate(), TSAdaptDestroy()
294: @*/
295: PetscErrorCode  TSAdaptReset(TSAdapt adapt)
296: {

301:   if (adapt->ops->reset) {(*adapt->ops->reset)(adapt);}
302:   return(0);
303: }

305: PetscErrorCode  TSAdaptDestroy(TSAdapt *adapt)
306: {

310:   if (!*adapt) return(0);
312:   if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return(0);}

314:   TSAdaptReset(*adapt);

316:   if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
317:   PetscViewerDestroy(&(*adapt)->monitor);
318:   PetscHeaderDestroy(adapt);
319:   return(0);
320: }

322: /*@
323:    TSAdaptSetMonitor - Monitor the choices made by the adaptive controller

325:    Collective on TSAdapt

327:    Input Arguments:
328: +  adapt - adaptive controller context
329: -  flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable

331:    Options Database Keys:
332: .  -ts_adapt_monitor - to turn on monitoring

334:    Level: intermediate

336: .seealso: TSAdaptChoose()
337: @*/
338: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
339: {

345:   if (flg) {
346:     if (!adapt->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);}
347:   } else {
348:     PetscViewerDestroy(&adapt->monitor);
349:   }
350:   return(0);
351: }

353: /*@C
354:    TSAdaptSetCheckStage - Set a callback to check convergence for a stage

356:    Logically collective on TSAdapt

358:    Input Arguments:
359: +  adapt - adaptive controller context
360: -  func - stage check function

362:    Arguments of func:
363: $  PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)

365: +  adapt - adaptive controller context
366: .  ts - time stepping context
367: -  accept - pending choice of whether to accept, can be modified by this routine

369:    Level: advanced

371: .seealso: TSAdaptChoose()
372: @*/
373: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*))
374: {

378:   adapt->checkstage = func;
379:   return(0);
380: }

382: /*@
383:    TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of
384:    any error or stability condition not meeting the prescribed goal.

386:    Logically collective on TSAdapt

388:    Input Arguments:
389: +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
390: -  flag - whether to always accept steps

392:    Options Database Keys:
393: .  -ts_adapt_always_accept - to always accept steps

395:    Level: intermediate

397: .seealso: TSAdapt, TSAdaptChoose()
398: @*/
399: PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag)
400: {
404:   adapt->always_accept = flag;
405:   return(0);
406: }

408: /*@
409:    TSAdaptSetSafety - Set safety factors

411:    Logically collective on TSAdapt

413:    Input Arguments:
414: +  adapt - adaptive controller context
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:    Options Database Keys:
419: +  -ts_adapt_safety <safety> - to set safety factor
420: -  -ts_adapt_reject_safety <reject_safety> - to set reject safety factor

422:    Level: intermediate

424: .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose()
425: @*/
426: PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety)
427: {
432:   if (safety != PETSC_DEFAULT && safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be non negative",(double)safety);
433:   if (safety != PETSC_DEFAULT && safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be less than one",(double)safety);
434:   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);
435:   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);
436:   if (safety != PETSC_DEFAULT) adapt->safety = safety;
437:   if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety;
438:   return(0);
439: }

441: /*@
442:    TSAdaptGetSafety - Get safety factors

444:    Not Collective

446:    Input Arguments:
447: .  adapt - adaptive controller context

449:    Ouput Arguments:
450: .  safety - safety factor relative to target error/stability goal
451: +  reject_safety - extra safety factor to apply if the last step was rejected

453:    Level: intermediate

455: .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose()
456: @*/
457: PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety)
458: {
463:   if (safety)        *safety        = adapt->safety;
464:   if (reject_safety) *reject_safety = adapt->reject_safety;
465:   return(0);
466: }

468: /*@
469:    TSAdaptSetClip - Sets the admissible decrease/increase factor in step size

471:    Logically collective on TSAdapt

473:    Input Arguments:
474: +  adapt - adaptive controller context
475: .  low - admissible decrease factor
476: -  high - admissible increase factor

478:    Options Database Keys:
479: .  -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors

481:    Level: intermediate

483: .seealso: TSAdaptChoose(), TSAdaptGetClip()
484: @*/
485: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
486: {
491:   if (low  != PETSC_DEFAULT && low  < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
492:   if (low  != PETSC_DEFAULT && low  > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
493:   if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high);
494:   if (low  != PETSC_DEFAULT) adapt->clip[0] = low;
495:   if (high != PETSC_DEFAULT) adapt->clip[1] = high;
496:   return(0);
497: }

499: /*@
500:    TSAdaptGetClip - Gets the admissible decrease/increase factor in step size

502:    Not Collective

504:    Input Arguments:
505: .  adapt - adaptive controller context

507:    Ouput Arguments:
508: +  low - optional, admissible decrease factor
509: -  high - optional, admissible increase factor

511:    Level: intermediate

513: .seealso: TSAdaptChoose(), TSAdaptSetClip()
514: @*/
515: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
516: {
521:   if (low)  *low  = adapt->clip[0];
522:   if (high) *high = adapt->clip[1];
523:   return(0);
524: }

526: /*@
527:    TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller

529:    Logically collective on TSAdapt

531:    Input Arguments:
532: +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
533: .  hmin - minimum time step
534: -  hmax - maximum time step

536:    Options Database Keys:
537: +  -ts_adapt_dt_min <min> - to set minimum time step
538: -  -ts_adapt_dt_max <max> - to set maximum time step

540:    Level: intermediate

542: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
543: @*/
544: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
545: {

551:   if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
552:   if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
553:   if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
554:   if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
555:   hmin = adapt->dt_min;
556:   hmax = adapt->dt_max;
557:   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);
558:   return(0);
559: }

561: /*@
562:    TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller

564:    Not Collective

566:    Input Arguments:
567: .  adapt - time step adaptivity context, usually gotten with TSGetAdapt()

569:    Output Arguments:
570: +  hmin - minimum time step
571: -  hmax - maximum time step

573:    Level: intermediate

575: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
576: @*/
577: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
578: {

584:   if (hmin) *hmin = adapt->dt_min;
585:   if (hmax) *hmax = adapt->dt_max;
586:   return(0);
587: }

589: /*
590:    TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.

592:    Collective on TSAdapt

594:    Input Parameter:
595: .  adapt - the TSAdapt context

597:    Options Database Keys:
598: +  -ts_adapt_type <type> - algorithm to use for adaptivity
599: .  -ts_adapt_always_accept - always accept steps regardless of error/stability goals
600: .  -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
601: .  -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
602: .  -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
603: .  -ts_adapt_dt_min <min> - minimum timestep to use
604: .  -ts_adapt_dt_max <max> - maximum timestep to use
605: .  -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
606: -  -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates

608:    Level: advanced

610:    Notes:
611:    This function is automatically called by TSSetFromOptions()

613: .keywords: TS, TSGetAdapt(), TSAdaptSetType(), TSAdaptSetStepLimits()

615: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
616:           TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
617: */
618: PetscErrorCode  TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
619: {
621:   char           type[256] = TSADAPTBASIC;
622:   PetscReal      safety,reject_safety,clip[2],hmin,hmax;
623:   PetscBool      set,flg;
624:   PetscInt       two;

628:   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
629:    * function can only be called from inside TSSetFromOptions()  */
630:   PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");
631:   PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
632:   if (flg || !((PetscObject)adapt)->type_name) {
633:     TSAdaptSetType(adapt,type);
634:   }

636:   PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);
637:   if (set) {TSAdaptSetAlwaysAccept(adapt,flg);}

639:   safety = adapt->safety; reject_safety = adapt->reject_safety;
640:   PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);
641:   PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);
642:   if (set || flg) {TSAdaptSetSafety(adapt,safety,reject_safety);}

644:   two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1];
645:   PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);
646:   if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip");
647:   if (set) {TSAdaptSetClip(adapt,clip[0],clip[1]);}

649:   hmin = adapt->dt_min; hmax = adapt->dt_max;
650:   PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
651:   PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
652:   if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}

654:   PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);

656:   PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);
657:   if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported");

659:   PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
660:   if (set) {TSAdaptSetMonitor(adapt,flg);}

662:   if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
663:   PetscOptionsTail();
664:   return(0);
665: }

667: /*@
668:    TSAdaptCandidatesClear - clear any previously set candidate schemes

670:    Logically collective on TSAdapt

672:    Input Argument:
673: .  adapt - adaptive controller

675:    Level: developer

677: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
678: @*/
679: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
680: {

685:   PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
686:   return(0);
687: }

689: /*@C
690:    TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from

692:    Logically collective on TSAdapt

694:    Input Arguments:
695: +  adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
696: .  name - name of the candidate scheme to add
697: .  order - order of the candidate scheme
698: .  stageorder - stage order of the candidate scheme
699: .  ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
700: .  cost - relative measure of the amount of work required for the candidate scheme
701: -  inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme

703:    Note:
704:    This routine is not available in Fortran.

706:    Level: developer

708: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
709: @*/
710: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
711: {
712:   PetscInt c;

716:   if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
717:   if (inuse) {
718:     if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
719:     adapt->candidates.inuse_set = PETSC_TRUE;
720:   }
721:   /* first slot if this is the current scheme, otherwise the next available slot */
722:   c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;

724:   adapt->candidates.name[c]       = name;
725:   adapt->candidates.order[c]      = order;
726:   adapt->candidates.stageorder[c] = stageorder;
727:   adapt->candidates.ccfl[c]       = ccfl;
728:   adapt->candidates.cost[c]       = cost;
729:   adapt->candidates.n++;
730:   return(0);
731: }

733: /*@C
734:    TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost

736:    Not Collective

738:    Input Arguments:
739: .  adapt - time step adaptivity context

741:    Output Arguments:
742: +  n - number of candidate schemes, always at least 1
743: .  order - the order of each candidate scheme
744: .  stageorder - the stage order of each candidate scheme
745: .  ccfl - the CFL coefficient of each scheme
746: -  cost - the relative cost of each scheme

748:    Level: developer

750:    Note:
751:    The current scheme is always returned in the first slot

753: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
754: @*/
755: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
756: {
759:   if (n) *n = adapt->candidates.n;
760:   if (order) *order = adapt->candidates.order;
761:   if (stageorder) *stageorder = adapt->candidates.stageorder;
762:   if (ccfl) *ccfl = adapt->candidates.ccfl;
763:   if (cost) *cost = adapt->candidates.cost;
764:   return(0);
765: }

767: /*@C
768:    TSAdaptChoose - choose which method and step size to use for the next step

770:    Collective on TSAdapt

772:    Input Arguments:
773: +  adapt - adaptive contoller
774: -  h - current step size

776:    Output Arguments:
777: +  next_sc - optional, scheme to use for the next step
778: .  next_h - step size to use for the next step
779: -  accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size

781:    Note:
782:    The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
783:    being retried after an initial rejection.

785:    Level: developer

787: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
788: @*/
789: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
790: {
792:   PetscInt       ncandidates = adapt->candidates.n;
793:   PetscInt       scheme = 0;
794:   PetscReal      wlte = -1.0;
795:   PetscReal      wltea = -1.0;
796:   PetscReal      wlter = -1.0;

804:   if (next_sc) *next_sc = 0;

806:   /* Do not mess with adaptivity while handling events*/
807:   if (ts->event && ts->event->status != TSEVENT_NONE) {
808:     *next_h = h;
809:     *accept = PETSC_TRUE;
810:     return(0);
811:   }

813:   (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
814:   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);
815:   if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
816:   if (next_sc) *next_sc = scheme;

818:   if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
819:     /* Increase/reduce step size if end time of next step is close to or overshoots max time */
820:     PetscReal t = ts->ptime + ts->time_step, h = *next_h;
821:     PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
822:     PetscReal a = (PetscReal)1.01; /* allow 1% step size increase */
823:     if (t < tmax && tend > tmax) *next_h = hmax;
824:     if (t < tmax && tend < tmax && h > hmax/2) *next_h = hmax/2;
825:     if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
826:   }

828:   if (adapt->monitor) {
829:     const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
830:     PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
831:     if (wlte < 0) {
832:       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);
833:     } else {
834:       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);
835:     }
836:     PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
837:   }
838:   return(0);
839: }

841: /*@
842:    TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)

844:    Collective on TSAdapt

846:    Input Arguments:
847: +  adapt - adaptive controller context
848: .  ts - time stepper
849: .  t - Current simulation time
850: -  Y - Current solution vector

852:    Output Arguments:
853: .  accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject

855:    Level: developer

857: .seealso:
858: @*/
859: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
860: {
861:   PetscErrorCode      ierr;
862:   SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;


869:   if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
870:   if (snesreason < 0) {
871:     *accept = PETSC_FALSE;
872:     if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
873:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
874:       PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
875:       if (adapt->monitor) {
876:         PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
877:         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);
878:         PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
879:       }
880:     }
881:   } else {
882:     *accept = PETSC_TRUE;
883:     TSFunctionDomainError(ts,t,Y,accept);
884:     if(*accept && adapt->checkstage) {
885:       (*adapt->checkstage)(adapt,ts,t,Y,accept);
886:     }
887:   }

889:   if(!(*accept) && !ts->reason) {
890:     PetscReal dt,new_dt;
891:     TSGetTimeStep(ts,&dt);
892:     new_dt = dt * adapt->scale_solve_failed;
893:     TSSetTimeStep(ts,new_dt);
894:     adapt->timestepjustincreased += 4;
895:     if (adapt->monitor) {
896:       PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
897:       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);
898:       PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
899:     }
900:   }
901:   return(0);
902: }

904: /*@
905:   TSAdaptCreate - create an adaptive controller context for time stepping

907:   Collective on MPI_Comm

909:   Input Parameter:
910: . comm - The communicator

912:   Output Parameter:
913: . adapt - new TSAdapt object

915:   Level: developer

917:   Notes:
918:   TSAdapt creation is handled by TS, so users should not need to call this function.

920: .keywords: TSAdapt, create
921: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
922: @*/
923: PetscErrorCode  TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
924: {
926:   TSAdapt        adapt;

930:   *inadapt = NULL;
931:   TSAdaptInitializePackage();

933:   PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);

935:   adapt->always_accept      = PETSC_FALSE;
936:   adapt->safety             = 0.9;
937:   adapt->reject_safety      = 0.5;
938:   adapt->clip[0]            = 0.1;
939:   adapt->clip[1]            = 10.;
940:   adapt->dt_min             = 1e-20;
941:   adapt->dt_max             = 1e+20;
942:   adapt->scale_solve_failed = 0.25;
943:   adapt->wnormtype          = NORM_2;

945:   *inadapt = adapt;
946:   return(0);
947: }