Actual source code: tsadapt.c

petsc-3.8.2 2017-11-09
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:   PetscFunctionListAdd(&TSAdaptList,sname,function);
 50:   return(0);
 51: }

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

 56:   Not Collective

 58:   Level: advanced

 60: .keywords: TSAdapt, register, all

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

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

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

 83:   Level: developer

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

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

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

104:   Level: developer

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

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

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

125:   Logicially Collective on TSAdapt

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

131:   Options Database:
132: .  -ts_adapt_type basic or none - to setting the adapter type

134:   Level: intermediate

136: .keywords: TSAdapt, create

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

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

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

162:   Not Collective

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

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

170:   Level: intermediate

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

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

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

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

197:   Collective on PetscViewer

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

205:    Level: intermediate

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

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

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

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

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

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

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

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

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

285:    Collective on TS

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

290:    Level: developer

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

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

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

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

313:   TSAdaptReset(*adapt);

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

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

324:    Collective on TSAdapt

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

330:    Options Database Keys:
331: .  -ts_adapt_monitor

333:    Level: intermediate

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

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

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

355:    Logically collective on TSAdapt

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

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

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

368:    Level: advanced

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

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

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

385:    Logically collective on TSAdapt

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

391:    Options Database Keys:
392: .  -ts_adapt_always_accept

394:    Level: intermediate

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

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

410:    Logically collective on TSAdapt

412:    Input Arguments:
413: +  adapt - adaptive controller context
414: .  safety - safety factor relative to target error/stability goal
415: +  reject_safety - extra safety factor to apply if the last step was rejected

417:    Options Database Keys:
418: +  -ts_adapt_safety
419: -  -ts_adapt_reject_safety

421:    Level: intermediate

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

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

443:    Not Collective

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

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

452:    Level: intermediate

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

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

470:    Logically collective on TSAdapt

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

477:    Options Database Keys:
478: .  -ts_adapt_clip

480:    Level: intermediate

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

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

501:    Not Collective

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

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

510:    Level: intermediate

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

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

528:    Logically collective on TSAdapt

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

535:    Options Database Keys:
536: +  -ts_adapt_dt_min - minimum time step
537: -  -ts_adapt_dt_max - maximum time step

539:    Level: intermediate

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

550:   if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
551:   if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
552:   if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
553:   if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
554:   hmin = adapt->dt_min;
555:   hmax = adapt->dt_max;
556:   if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must geather than minimum time step %g",(double)hmax,(double)hmin);
557:   return(0);
558: }

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

563:    Not Collective

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

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

572:    Level: intermediate

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

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

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

591:    Collective on TSAdapt

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

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

607:    Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

669:    Logically collective on TSAdapt

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

674:    Level: developer

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

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

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

691:    Logically collective on TSAdapt

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

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

705:    Level: developer

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

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

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

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

735:    Not Collective

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

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

747:    Level: developer

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

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

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

769:    Collective on TSAdapt

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

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

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

784:    Level: developer

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

803:   if (next_sc) *next_sc = 0;

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

812:   (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);
813:   if (scheme < 0 || (ncandidates > 0 && scheme >= ncandidates)) SETERRQ2(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",scheme,ncandidates-1);
814:   if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h);
815:   if (next_sc) *next_sc = scheme;

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

827:   if (adapt->monitor) {
828:     const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
829:     PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
830:     if (wlte < 0) {
831:       PetscViewerASCIIPrintf(adapt->monitor,"    TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h);
832:     } else {
833:       PetscViewerASCIIPrintf(adapt->monitor,"    TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e wlte=%5.3g  wltea=%5.3g wlter=%5.3g\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h,(double)wlte,(double)wltea,(double)wlter);
834:     }
835:     PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
836:   }
837:   return(0);
838: }

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

843:    Collective on TSAdapt

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

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

854:    Level: developer

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


868:   if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
869:   if (snesreason < 0) {
870:     *accept = PETSC_FALSE;
871:     if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
872:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
873:       PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
874:       if (adapt->monitor) {
875:         PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
876:         PetscViewerASCIIPrintf(adapt->monitor,"    TSAdapt %s step %3D stage rejected t=%-11g+%10.3e, nonlinear solve failures %D greater than current TS allowed\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)ts->time_step,ts->num_snes_failures);
877:         PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
878:       }
879:     }
880:   } else {
881:     *accept = PETSC_TRUE;
882:     TSFunctionDomainError(ts,t,Y,accept);
883:     if(*accept && adapt->checkstage) {
884:       (*adapt->checkstage)(adapt,ts,t,Y,accept);
885:     }
886:   }

888:   if(!(*accept) && !ts->reason) {
889:     PetscReal dt,new_dt;
890:     TSGetTimeStep(ts,&dt);
891:     new_dt = dt * adapt->scale_solve_failed;
892:     TSSetTimeStep(ts,new_dt);
893:     adapt->timestepjustincreased += 4;
894:     if (adapt->monitor) {
895:       PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
896:       PetscViewerASCIIPrintf(adapt->monitor,"    TSAdapt %s step %3D stage rejected (%s) t=%-11g+%10.3e retrying with dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,SNESConvergedReasons[snesreason],(double)ts->ptime,(double)dt,(double)new_dt);
897:       PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
898:     }
899:   }
900:   return(0);
901: }

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

906:   Collective on MPI_Comm

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

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

914:   Level: developer

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

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

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

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

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

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