Actual source code: tsadapt.c

petsc-3.11.4 2019-09-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);
 15: PETSC_EXTERN PetscErrorCode TSAdaptCreate_History(TSAdapt);

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

 20:    Not Collective

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

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

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

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

 39:    Level: advanced

 41: .keywords: TSAdapt, register

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

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

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

 58:   Not Collective

 60:   Level: advanced

 62: .keywords: TSAdapt, register, all

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

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

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

 86:   Level: developer

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

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

102: /*@C
103:   TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
104:   called from TSInitializePackage().

106:   Level: developer

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

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

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

127:   Logicially Collective on TSAdapt

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

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

136:   Level: intermediate

138: .keywords: TSAdapt, create

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

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

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

164:   Not Collective

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

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

172:   Level: intermediate

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

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

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

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

199:   Collective on PetscViewer

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

207:    Level: intermediate

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

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

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

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

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

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

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

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

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

287:    Collective on TS

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

292:    Level: developer

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

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

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

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

315:   TSAdaptReset(*adapt);

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

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

326:    Collective on TSAdapt

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

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

335:    Level: intermediate

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

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

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

357:    Logically collective on TSAdapt

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

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

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

370:    Level: advanced

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

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

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

387:    Logically collective on TSAdapt

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

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

396:    Level: intermediate

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

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

412:    Logically collective on TSAdapt

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

419:    Options Database Keys:
420: +  -ts_adapt_safety <safety> - to set safety factor
421: -  -ts_adapt_reject_safety <reject_safety> - to set reject safety factor

423:    Level: intermediate

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

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

445:    Not Collective

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

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

454:    Level: intermediate

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

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

472:    Logically collective on TSAdapt

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

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

482:    Level: intermediate

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

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

503:    Not Collective

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

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

512:    Level: intermediate

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

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

530:    Logically collective on TSAdapt

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

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

541:    Level: intermediate

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

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

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

565:    Not Collective

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

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

574:    Level: intermediate

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

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

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

593:    Collective on TSAdapt

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

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

610:    Level: advanced

612:    Notes:
613:    This function is automatically called by TSSetFromOptions()

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

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

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

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

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

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

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

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

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

661:   PetscOptionsInt("-ts_adapt_time_step_increase_delay","Number of timesteps to delay increasing the time step after it has been decreased due to failed solver","TSAdaptSetTimeStepIncreaseDelay",adapt->timestepjustdecreased_delay,&adapt->timestepjustdecreased_delay,NULL);
662: 
663:   PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
664:   if (set) {TSAdaptSetMonitor(adapt,flg);}

666:   if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
667:   PetscOptionsTail();
668:   return(0);
669: }

671: /*@
672:    TSAdaptCandidatesClear - clear any previously set candidate schemes

674:    Logically collective on TSAdapt

676:    Input Argument:
677: .  adapt - adaptive controller

679:    Level: developer

681: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
682: @*/
683: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
684: {

689:   PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
690:   return(0);
691: }

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

696:    Logically collective on TSAdapt

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

707:    Note:
708:    This routine is not available in Fortran.

710:    Level: developer

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

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

728:   adapt->candidates.name[c]       = name;
729:   adapt->candidates.order[c]      = order;
730:   adapt->candidates.stageorder[c] = stageorder;
731:   adapt->candidates.ccfl[c]       = ccfl;
732:   adapt->candidates.cost[c]       = cost;
733:   adapt->candidates.n++;
734:   return(0);
735: }

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

740:    Not Collective

742:    Input Arguments:
743: .  adapt - time step adaptivity context

745:    Output Arguments:
746: +  n - number of candidate schemes, always at least 1
747: .  order - the order of each candidate scheme
748: .  stageorder - the stage order of each candidate scheme
749: .  ccfl - the CFL coefficient of each scheme
750: -  cost - the relative cost of each scheme

752:    Level: developer

754:    Note:
755:    The current scheme is always returned in the first slot

757: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
758: @*/
759: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
760: {
763:   if (n) *n = adapt->candidates.n;
764:   if (order) *order = adapt->candidates.order;
765:   if (stageorder) *stageorder = adapt->candidates.stageorder;
766:   if (ccfl) *ccfl = adapt->candidates.ccfl;
767:   if (cost) *cost = adapt->candidates.cost;
768:   return(0);
769: }

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

774:    Collective on TSAdapt

776:    Input Arguments:
777: +  adapt - adaptive contoller
778: -  h - current step size

780:    Output Arguments:
781: +  next_sc - optional, scheme to use for the next step
782: .  next_h - step size to use for the next step
783: -  accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size

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

789:    Level: developer

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

808:   if (next_sc) *next_sc = 0;

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

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

822:   if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
823:     /* Increase/reduce step size if end time of next step is close to or overshoots max time */
824:     PetscReal t = ts->ptime + ts->time_step, h = *next_h;
825:     PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
826:     PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
827:     PetscReal b = adapt->matchstepfac[1];
828:     if (t < tmax && tend > tmax) *next_h = hmax;
829:     if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
830:     if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
831:   }

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

846: /*@
847:    TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver
848:                                      before increasing the time step.

850:    Logicially Collective on TSAdapt

852:    Input Arguments:
853: +  adapt - adaptive controller context
854: -  cnt - the number of timesteps

856:    Options Database Key:
857: .  -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase

859:    Notes: This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0.
860:           The successful use of this option is problem dependent

862:    Developer Note: there is no theory to support this option

864:    Level: advanced

866: .seealso:
867: @*/
868: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
869: {
871:   adapt->timestepjustdecreased_delay = cnt;
872:   return(0);
873: }

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

879:    Collective on TSAdapt

881:    Input Arguments:
882: +  adapt - adaptive controller context
883: .  ts - time stepper
884: .  t - Current simulation time
885: -  Y - Current solution vector

887:    Output Arguments:
888: .  accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject

890:    Level: developer

892: .seealso:
893: @*/
894: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
895: {
896:   PetscErrorCode      ierr;
897:   SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;


904:   if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
905:   if (snesreason < 0) {
906:     *accept = PETSC_FALSE;
907:     if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
908:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
909:       PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
910:       if (adapt->monitor) {
911:         PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
912:         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);
913:         PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
914:       }
915:     }
916:   } else {
917:     *accept = PETSC_TRUE;
918:     TSFunctionDomainError(ts,t,Y,accept);
919:     if(*accept && adapt->checkstage) {
920:       (*adapt->checkstage)(adapt,ts,t,Y,accept);
921:     }
922:   }

924:   if(!(*accept) && !ts->reason) {
925:     PetscReal dt,new_dt;
926:     TSGetTimeStep(ts,&dt);
927:     new_dt = dt * adapt->scale_solve_failed;
928:     TSSetTimeStep(ts,new_dt);
929:     adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
930:     if (adapt->monitor) {
931:       PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
932:       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);
933:       PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
934:     }
935:   }
936:   return(0);
937: }

939: /*@
940:   TSAdaptCreate - create an adaptive controller context for time stepping

942:   Collective on MPI_Comm

944:   Input Parameter:
945: . comm - The communicator

947:   Output Parameter:
948: . adapt - new TSAdapt object

950:   Level: developer

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

955: .keywords: TSAdapt, create
956: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
957: @*/
958: PetscErrorCode  TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
959: {
961:   TSAdapt        adapt;

965:   *inadapt = NULL;
966:   TSAdaptInitializePackage();

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

970:   adapt->always_accept      = PETSC_FALSE;
971:   adapt->safety             = 0.9;
972:   adapt->reject_safety      = 0.5;
973:   adapt->clip[0]            = 0.1;
974:   adapt->clip[1]            = 10.;
975:   adapt->dt_min             = 1e-20;
976:   adapt->dt_max             = 1e+20;
977:   adapt->scale_solve_failed = 0.25;
978:   /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
979:      to prevent from situations were unreasonably small time steps are taken in order to match the final time */
980:   adapt->matchstepfac[0]    = 0.01; /* allow 1% step size increase in the last step */
981:   adapt->matchstepfac[1]    = 2.0;  /* halve last step if it is greater than what remains divided this factor */
982:   adapt->wnormtype          = NORM_2;
983:   adapt->timestepjustdecreased_delay = 0;

985:   *inadapt = adapt;
986:   return(0);
987: }