Actual source code: tsadapt.c

petsc-master 2019-08-17
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: .seealso: TSAdaptRegisterAll()
 42: @*/
 43: PetscErrorCode  TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
 44: {

 48:   TSAdaptInitializePackage();
 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: .seealso: TSAdaptRegisterDestroy()
 61: @*/
 62: PetscErrorCode  TSAdaptRegisterAll(void)
 63: {

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

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

 82:   Level: developer

 84: .seealso: PetscFinalize()
 85: @*/
 86: PetscErrorCode  TSAdaptFinalizePackage(void)
 87: {

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

 97: /*@C
 98:   TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
 99:   called from TSInitializePackage().

101:   Level: developer

103: .seealso: PetscInitialize()
104: @*/
105: PetscErrorCode  TSAdaptInitializePackage(void)
106: {

110:   if (TSAdaptPackageInitialized) return(0);
111:   TSAdaptPackageInitialized = PETSC_TRUE;
112:   PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
113:   TSAdaptRegisterAll();
114:   PetscRegisterFinalize(TSAdaptFinalizePackage);
115:   return(0);
116: }

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

121:   Logicially Collective on TSAdapt

123:   Input Parameter:
124: + adapt - the TS adapter, most likely obtained with TSGetAdapt()
125: - type - either  TSADAPTBASIC or TSADAPTNONE

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

130:   Level: intermediate

132: .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType()
133: @*/
134: PetscErrorCode  TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
135: {
136:   PetscBool      match;
137:   PetscErrorCode ierr,(*r)(TSAdapt);

142:   PetscObjectTypeCompare((PetscObject)adapt,type,&match);
143:   if (match) return(0);
144:   PetscFunctionListFind(TSAdaptList,type,&r);
145:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
146:   if (adapt->ops->destroy) {(*adapt->ops->destroy)(adapt);}
147:   PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));
148:   PetscObjectChangeTypeName((PetscObject)adapt,type);
149:   (*r)(adapt);
150:   return(0);
151: }

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

156:   Not Collective

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

161:   Output Parameter:
162: . type - The name of TS adapter method

164:   Level: intermediate

166: .seealso TSAdaptSetType()
167: @*/
168: PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type)
169: {
173:   *type = ((PetscObject)adapt)->type_name;
174:   return(0);
175: }

177: PetscErrorCode  TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
178: {

183:   PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
184:   return(0);
185: }

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

190:   Collective on PetscViewer

192:   Input Parameters:
193: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
194:            some related function before a call to TSAdaptLoad().
195: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
196:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

198:    Level: intermediate

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

203:   Notes for advanced users:
204:   Most users should not need to know the details of the binary storage
205:   format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
206:   But for anyone who's interested, the standard binary matrix storage
207:   format is
208: .vb
209:      has not yet been determined
210: .ve

212: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
213: @*/
214: PetscErrorCode  TSAdaptLoad(TSAdapt adapt,PetscViewer viewer)
215: {
217:   PetscBool      isbinary;
218:   char           type[256];

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

226:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
227:   TSAdaptSetType(adapt,type);
228:   if (adapt->ops->load) {
229:     (*adapt->ops->load)(adapt,viewer);
230:   }
231:   return(0);
232: }

234: PetscErrorCode  TSAdaptView(TSAdapt adapt,PetscViewer viewer)
235: {
237:   PetscBool      iascii,isbinary,isnone,isglee;

241:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);}
244:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
245:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
246:   if (iascii) {
247:     PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
248:     PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);
249:     PetscObjectTypeCompare((PetscObject)adapt,TSADAPTGLEE,&isglee);
250:     if (!isnone) {
251:       if (adapt->always_accept) {PetscViewerASCIIPrintf(viewer,"  always accepting steps\n");}
252:       PetscViewerASCIIPrintf(viewer,"  safety factor %g\n",(double)adapt->safety);
253:       PetscViewerASCIIPrintf(viewer,"  extra safety factor after step rejection %g\n",(double)adapt->reject_safety);
254:       PetscViewerASCIIPrintf(viewer,"  clip fastest increase %g\n",(double)adapt->clip[1]);
255:       PetscViewerASCIIPrintf(viewer,"  clip fastest decrease %g\n",(double)adapt->clip[0]);
256:       PetscViewerASCIIPrintf(viewer,"  maximum allowed timestep %g\n",(double)adapt->dt_max);
257:       PetscViewerASCIIPrintf(viewer,"  minimum allowed timestep %g\n",(double)adapt->dt_min);
258:       PetscViewerASCIIPrintf(viewer,"  maximum solution absolute value to be ignored %g\n",(double)adapt->ignore_max);
259:     }
260:     if (isglee) {
261:       if (adapt->glee_use_local) {
262:         PetscViewerASCIIPrintf(viewer,"  GLEE uses local error control\n");
263:       } else {
264:         PetscViewerASCIIPrintf(viewer,"  GLEE uses global error control\n");
265:       }
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:    TSAdaptSetMaxIgnore - Set error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value). A negative value (default) of the threshold leads to considering all solution components.

472:    Logically collective on TSAdapt

474:    Input Arguments:
475: +  adapt - adaptive controller context
476: -  max_ignore - threshold for solution components that are ignored during error estimation

478:    Options Database Keys:
479: .  -ts_adapt_max_ignore <max_ignore> - to set the threshold

481:    Level: intermediate

483: .seealso: TSAdapt, TSAdaptGetMaxIgnore(), TSAdaptChoose()
484: @*/
485: PetscErrorCode TSAdaptSetMaxIgnore(TSAdapt adapt,PetscReal max_ignore)
486: {
490:   adapt->ignore_max = max_ignore;
491:   return(0);
492: }

494: /*@
495:    TSAdaptGetMaxIgnore - Get error estimation threshold. Solution components below this threshold value will not be considered when computing error norms for time step adaptivity (in absolute value).

497:    Not Collective

499:    Input Arguments:
500: .  adapt - adaptive controller context

502:    Ouput Arguments:
503: .  max_ignore - threshold for solution components that are ignored during error estimation

505:    Level: intermediate

507: .seealso: TSAdapt, TSAdaptSetMaxIgnore(), TSAdaptChoose()
508: @*/
509: PetscErrorCode TSAdaptGetMaxIgnore(TSAdapt adapt,PetscReal *max_ignore)
510: {
514:   *max_ignore = adapt->ignore_max;
515:   return(0);
516: }


519: /*@
520:    TSAdaptSetClip - Sets the admissible decrease/increase factor in step size

522:    Logically collective on TSAdapt

524:    Input Arguments:
525: +  adapt - adaptive controller context
526: .  low - admissible decrease factor
527: -  high - admissible increase factor

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

532:    Level: intermediate

534: .seealso: TSAdaptChoose(), TSAdaptGetClip()
535: @*/
536: PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high)
537: {
542:   if (low  != PETSC_DEFAULT && low  < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low);
543:   if (low  != PETSC_DEFAULT && low  > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low);
544:   if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high);
545:   if (low  != PETSC_DEFAULT) adapt->clip[0] = low;
546:   if (high != PETSC_DEFAULT) adapt->clip[1] = high;
547:   return(0);
548: }

550: /*@
551:    TSAdaptGetClip - Gets the admissible decrease/increase factor in step size

553:    Not Collective

555:    Input Arguments:
556: .  adapt - adaptive controller context

558:    Ouput Arguments:
559: +  low - optional, admissible decrease factor
560: -  high - optional, admissible increase factor

562:    Level: intermediate

564: .seealso: TSAdaptChoose(), TSAdaptSetClip()
565: @*/
566: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
567: {
572:   if (low)  *low  = adapt->clip[0];
573:   if (high) *high = adapt->clip[1];
574:   return(0);
575: }

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

580:    Logically collective on TSAdapt

582:    Input Arguments:
583: +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
584: .  hmin - minimum time step
585: -  hmax - maximum time step

587:    Options Database Keys:
588: +  -ts_adapt_dt_min <min> - to set minimum time step
589: -  -ts_adapt_dt_max <max> - to set maximum time step

591:    Level: intermediate

593: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
594: @*/
595: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
596: {

602:   if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
603:   if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
604:   if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
605:   if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
606:   hmin = adapt->dt_min;
607:   hmax = adapt->dt_max;
608:   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);
609:   return(0);
610: }

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

615:    Not Collective

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

620:    Output Arguments:
621: +  hmin - minimum time step
622: -  hmax - maximum time step

624:    Level: intermediate

626: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
627: @*/
628: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
629: {

635:   if (hmin) *hmin = adapt->dt_min;
636:   if (hmax) *hmax = adapt->dt_max;
637:   return(0);
638: }

640: /*
641:    TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.

643:    Collective on TSAdapt

645:    Input Parameter:
646: .  adapt - the TSAdapt context

648:    Options Database Keys:
649: +  -ts_adapt_type <type> - algorithm to use for adaptivity
650: .  -ts_adapt_always_accept - always accept steps regardless of error/stability goals
651: .  -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
652: .  -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
653: .  -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
654: .  -ts_adapt_dt_min <min> - minimum timestep to use
655: .  -ts_adapt_dt_max <max> - maximum timestep to use
656: .  -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
657: .  -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
658: -  -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver

660:    Level: advanced

662:    Notes:
663:    This function is automatically called by TSSetFromOptions()

665: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
666:           TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
667: */
668: PetscErrorCode  TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
669: {
671:   char           type[256] = TSADAPTBASIC;
672:   PetscReal      safety,reject_safety,clip[2],hmin,hmax;
673:   PetscBool      set,flg;
674:   PetscInt       two;

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

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

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

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

699:   hmin = adapt->dt_min; hmax = adapt->dt_max;
700:   PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
701:   PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
702:   if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}

704:   PetscOptionsReal("-ts_adapt_max_ignore","Adaptor ignores (absolute) solution values smaller than this value","",adapt->ignore_max,&adapt->ignore_max,&set);
705:   PetscOptionsBool("-ts_adapt_glee_use_local","GLEE adaptor uses local error estimation for step control","",adapt->glee_use_local,&adapt->glee_use_local,&set);

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

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

712:   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);

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

717:   if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
718:   PetscOptionsTail();
719:   return(0);
720: }

722: /*@
723:    TSAdaptCandidatesClear - clear any previously set candidate schemes

725:    Logically collective on TSAdapt

727:    Input Argument:
728: .  adapt - adaptive controller

730:    Level: developer

732: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
733: @*/
734: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
735: {

740:   PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
741:   return(0);
742: }

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

747:    Logically collective on TSAdapt

749:    Input Arguments:
750: +  adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
751: .  name - name of the candidate scheme to add
752: .  order - order of the candidate scheme
753: .  stageorder - stage order of the candidate scheme
754: .  ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
755: .  cost - relative measure of the amount of work required for the candidate scheme
756: -  inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme

758:    Note:
759:    This routine is not available in Fortran.

761:    Level: developer

763: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
764: @*/
765: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
766: {
767:   PetscInt c;

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

779:   adapt->candidates.name[c]       = name;
780:   adapt->candidates.order[c]      = order;
781:   adapt->candidates.stageorder[c] = stageorder;
782:   adapt->candidates.ccfl[c]       = ccfl;
783:   adapt->candidates.cost[c]       = cost;
784:   adapt->candidates.n++;
785:   return(0);
786: }

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

791:    Not Collective

793:    Input Arguments:
794: .  adapt - time step adaptivity context

796:    Output Arguments:
797: +  n - number of candidate schemes, always at least 1
798: .  order - the order of each candidate scheme
799: .  stageorder - the stage order of each candidate scheme
800: .  ccfl - the CFL coefficient of each scheme
801: -  cost - the relative cost of each scheme

803:    Level: developer

805:    Note:
806:    The current scheme is always returned in the first slot

808: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
809: @*/
810: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
811: {
814:   if (n) *n = adapt->candidates.n;
815:   if (order) *order = adapt->candidates.order;
816:   if (stageorder) *stageorder = adapt->candidates.stageorder;
817:   if (ccfl) *ccfl = adapt->candidates.ccfl;
818:   if (cost) *cost = adapt->candidates.cost;
819:   return(0);
820: }

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

825:    Collective on TSAdapt

827:    Input Arguments:
828: +  adapt - adaptive contoller
829: -  h - current step size

831:    Output Arguments:
832: +  next_sc - optional, scheme to use for the next step
833: .  next_h - step size to use for the next step
834: -  accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size

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

840:    Level: developer

842: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
843: @*/
844: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
845: {
847:   PetscInt       ncandidates = adapt->candidates.n;
848:   PetscInt       scheme = 0;
849:   PetscReal      wlte = -1.0;
850:   PetscReal      wltea = -1.0;
851:   PetscReal      wlter = -1.0;

859:   if (next_sc) *next_sc = 0;

861:   /* Do not mess with adaptivity while handling events*/
862:   if (ts->event && ts->event->status != TSEVENT_NONE) {
863:     *next_h = h;
864:     *accept = PETSC_TRUE;
865:     return(0);
866:   }

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

873:   if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
874:     /* Increase/reduce step size if end time of next step is close to or overshoots max time */
875:     PetscReal t = ts->ptime + ts->time_step, h = *next_h;
876:     PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
877:     PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
878:     PetscReal b = adapt->matchstepfac[1];
879:     if (t < tmax && tend > tmax) *next_h = hmax;
880:     if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
881:     if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
882:   }

884:   if (adapt->monitor) {
885:     const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
886:     PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
887:     if (wlte < 0) {
888:       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);
889:     } else {
890:       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);
891:     }
892:     PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
893:   }
894:   return(0);
895: }

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

901:    Logicially Collective on TSAdapt

903:    Input Arguments:
904: +  adapt - adaptive controller context
905: -  cnt - the number of timesteps

907:    Options Database Key:
908: .  -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase

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

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

915:    Level: advanced

917: .seealso:
918: @*/
919: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
920: {
922:   adapt->timestepjustdecreased_delay = cnt;
923:   return(0);
924: }

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

930:    Collective on TSAdapt

932:    Input Arguments:
933: +  adapt - adaptive controller context
934: .  ts - time stepper
935: .  t - Current simulation time
936: -  Y - Current solution vector

938:    Output Arguments:
939: .  accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject

941:    Level: developer

943: .seealso:
944: @*/
945: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
946: {
947:   PetscErrorCode      ierr;
948:   SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;


955:   if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
956:   if (snesreason < 0) {
957:     *accept = PETSC_FALSE;
958:     if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
959:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
960:       PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
961:       if (adapt->monitor) {
962:         PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
963:         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);
964:         PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
965:       }
966:     }
967:   } else {
968:     *accept = PETSC_TRUE;
969:     TSFunctionDomainError(ts,t,Y,accept);
970:     if(*accept && adapt->checkstage) {
971:       (*adapt->checkstage)(adapt,ts,t,Y,accept);
972:     }
973:   }

975:   if(!(*accept) && !ts->reason) {
976:     PetscReal dt,new_dt;
977:     TSGetTimeStep(ts,&dt);
978:     new_dt = dt * adapt->scale_solve_failed;
979:     TSSetTimeStep(ts,new_dt);
980:     adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
981:     if (adapt->monitor) {
982:       PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
983:       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);
984:       PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
985:     }
986:   }
987:   return(0);
988: }

990: /*@
991:   TSAdaptCreate - create an adaptive controller context for time stepping

993:   Collective

995:   Input Parameter:
996: . comm - The communicator

998:   Output Parameter:
999: . adapt - new TSAdapt object

1001:   Level: developer

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

1006: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1007: @*/
1008: PetscErrorCode  TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1009: {
1011:   TSAdapt        adapt;

1015:   *inadapt = NULL;
1016:   TSAdaptInitializePackage();

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

1020:   adapt->always_accept      = PETSC_FALSE;
1021:   adapt->safety             = 0.9;
1022:   adapt->reject_safety      = 0.5;
1023:   adapt->clip[0]            = 0.1;
1024:   adapt->clip[1]            = 10.;
1025:   adapt->dt_min             = 1e-20;
1026:   adapt->dt_max             = 1e+20;
1027:   adapt->ignore_max         = -1.0;
1028:   adapt->glee_use_local     = PETSC_TRUE;
1029:   adapt->scale_solve_failed = 0.25;
1030:   /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1031:      to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1032:   adapt->matchstepfac[0]    = 0.01; /* allow 1% step size increase in the last step */
1033:   adapt->matchstepfac[1]    = 2.0;  /* halve last step if it is greater than what remains divided this factor */
1034:   adapt->wnormtype          = NORM_2;
1035:   adapt->timestepjustdecreased_delay = 0;

1037:   *inadapt = adapt;
1038:   return(0);
1039: }