Actual source code: tsadapt.c


  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 Parameters:
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);
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 Parameters:
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 Parameters:
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: {
378:   adapt->checkstage = func;
379:   return(0);
380: }

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

386:    Logically collective on TSAdapt

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

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

395:    Level: intermediate

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

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

411:    Logically collective on TSAdapt

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

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

422:    Level: intermediate

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

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

444:    Not Collective

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

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

453:    Level: intermediate

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

468: /*@
469:    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.

471:    Logically collective on TSAdapt

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

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

480:    Level: intermediate

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

493: /*@
494:    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).

496:    Not Collective

498:    Input Parameter:
499: .  adapt - adaptive controller context

501:    Output Parameter:
502: .  max_ignore - threshold for solution components that are ignored during error estimation

504:    Level: intermediate

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

517: /*@
518:    TSAdaptSetClip - Sets the admissible decrease/increase factor in step size

520:    Logically collective on TSAdapt

522:    Input Parameters:
523: +  adapt - adaptive controller context
524: .  low - admissible decrease factor
525: -  high - admissible increase factor

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

530:    Level: intermediate

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

548: /*@
549:    TSAdaptGetClip - Gets the admissible decrease/increase factor in step size

551:    Not Collective

553:    Input Parameter:
554: .  adapt - adaptive controller context

556:    Output Parameters:
557: +  low - optional, admissible decrease factor
558: -  high - optional, admissible increase factor

560:    Level: intermediate

562: .seealso: TSAdaptChoose(), TSAdaptSetClip(), TSAdaptSetScaleSolveFailed()
563: @*/
564: PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high)
565: {
570:   if (low)  *low  = adapt->clip[0];
571:   if (high) *high = adapt->clip[1];
572:   return(0);
573: }

575: /*@
576:    TSAdaptSetScaleSolveFailed - Scale step by this factor if solve fails

578:    Logically collective on TSAdapt

580:    Input Parameters:
581: +  adapt - adaptive controller context
582: -  scale - scale

584:    Options Database Keys:
585: .  -ts_adapt_scale_solve_failed <scale> - to set scale step by this factor if solve fails

587:    Level: intermediate

589: .seealso: TSAdaptChoose(), TSAdaptGetScaleSolveFailed(), TSAdaptGetClip()
590: @*/
591: PetscErrorCode TSAdaptSetScaleSolveFailed(TSAdapt adapt,PetscReal scale)
592: {
596:   if (scale != PETSC_DEFAULT && scale <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be positive",(double)scale);
597:   if (scale != PETSC_DEFAULT && scale  > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale factor %g must be less than one",(double)scale);
598:   if (scale != PETSC_DEFAULT) adapt->scale_solve_failed = scale;
599:   return(0);
600: }

602: /*@
603:    TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size

605:    Not Collective

607:    Input Parameter:
608: .  adapt - adaptive controller context

610:    Output Parameter:
611: .  scale - scale factor

613:    Level: intermediate

615: .seealso: TSAdaptChoose(), TSAdaptSetScaleSolveFailed(), TSAdaptSetClip()
616: @*/
617: PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal *scale)
618: {
622:   if (scale)  *scale  = adapt->scale_solve_failed;
623:   return(0);
624: }

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

629:    Logically collective on TSAdapt

631:    Input Parameters:
632: +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
633: .  hmin - minimum time step
634: -  hmax - maximum time step

636:    Options Database Keys:
637: +  -ts_adapt_dt_min <min> - to set minimum time step
638: -  -ts_adapt_dt_max <max> - to set maximum time step

640:    Level: intermediate

642: .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose()
643: @*/
644: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
645: {

651:   if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
652:   if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
653:   if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
654:   if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
655:   hmin = adapt->dt_min;
656:   hmax = adapt->dt_max;
657:   if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must greater than minimum time step %g",(double)hmax,(double)hmin);
658:   return(0);
659: }

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

664:    Not Collective

666:    Input Parameter:
667: .  adapt - time step adaptivity context, usually gotten with TSGetAdapt()

669:    Output Parameters:
670: +  hmin - minimum time step
671: -  hmax - maximum time step

673:    Level: intermediate

675: .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose()
676: @*/
677: PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax)
678: {

684:   if (hmin) *hmin = adapt->dt_min;
685:   if (hmax) *hmax = adapt->dt_max;
686:   return(0);
687: }

689: /*
690:    TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.

692:    Collective on TSAdapt

694:    Input Parameter:
695: .  adapt - the TSAdapt context

697:    Options Database Keys:
698: +  -ts_adapt_type <type> - algorithm to use for adaptivity
699: .  -ts_adapt_always_accept - always accept steps regardless of error/stability goals
700: .  -ts_adapt_safety <safety> - safety factor relative to target error/stability goal
701: .  -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected
702: .  -ts_adapt_clip <low,high> - admissible time step decrease and increase factors
703: .  -ts_adapt_dt_min <min> - minimum timestep to use
704: .  -ts_adapt_dt_max <max> - maximum timestep to use
705: .  -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails
706: .  -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates
707: -  -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver

709:    Level: advanced

711:    Notes:
712:    This function is automatically called by TSSetFromOptions()

714: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(),
715:           TSAdaptSetClip(), TSAdaptSetScaleSolveFailed(), TSAdaptSetStepLimits(), TSAdaptSetMonitor()
716: */
717: PetscErrorCode  TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
718: {
720:   char           type[256] = TSADAPTBASIC;
721:   PetscReal      safety,reject_safety,clip[2],scale,hmin,hmax;
722:   PetscBool      set,flg;
723:   PetscInt       two;

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

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

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

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

748:   hmin = adapt->dt_min; hmax = adapt->dt_max;
749:   PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);
750:   PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);
751:   if (set || flg) {TSAdaptSetStepLimits(adapt,hmin,hmax);}

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

756:   PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","TSAdaptSetScaleSolveFailed",adapt->scale_solve_failed,&scale,&set);
757:   if (set) {TSAdaptSetScaleSolveFailed(adapt,scale);}

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

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

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

767:   if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
768:   PetscOptionsTail();
769:   return(0);
770: }

772: /*@
773:    TSAdaptCandidatesClear - clear any previously set candidate schemes

775:    Logically collective on TSAdapt

777:    Input Parameter:
778: .  adapt - adaptive controller

780:    Level: developer

782: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
783: @*/
784: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
785: {

790:   PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
791:   return(0);
792: }

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

797:    Logically collective on TSAdapt

799:    Input Parameters:
800: +  adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
801: .  name - name of the candidate scheme to add
802: .  order - order of the candidate scheme
803: .  stageorder - stage order of the candidate scheme
804: .  ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
805: .  cost - relative measure of the amount of work required for the candidate scheme
806: -  inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme

808:    Note:
809:    This routine is not available in Fortran.

811:    Level: developer

813: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
814: @*/
815: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
816: {
817:   PetscInt c;

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

829:   adapt->candidates.name[c]       = name;
830:   adapt->candidates.order[c]      = order;
831:   adapt->candidates.stageorder[c] = stageorder;
832:   adapt->candidates.ccfl[c]       = ccfl;
833:   adapt->candidates.cost[c]       = cost;
834:   adapt->candidates.n++;
835:   return(0);
836: }

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

841:    Not Collective

843:    Input Parameter:
844: .  adapt - time step adaptivity context

846:    Output Parameters:
847: +  n - number of candidate schemes, always at least 1
848: .  order - the order of each candidate scheme
849: .  stageorder - the stage order of each candidate scheme
850: .  ccfl - the CFL coefficient of each scheme
851: -  cost - the relative cost of each scheme

853:    Level: developer

855:    Note:
856:    The current scheme is always returned in the first slot

858: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
859: @*/
860: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
861: {
864:   if (n) *n = adapt->candidates.n;
865:   if (order) *order = adapt->candidates.order;
866:   if (stageorder) *stageorder = adapt->candidates.stageorder;
867:   if (ccfl) *ccfl = adapt->candidates.ccfl;
868:   if (cost) *cost = adapt->candidates.cost;
869:   return(0);
870: }

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

875:    Collective on TSAdapt

877:    Input Parameters:
878: +  adapt - adaptive contoller
879: .  ts - time stepper
880: -  h - current step size

882:    Output Parameters:
883: +  next_sc - optional, scheme to use for the next step
884: .  next_h - step size to use for the next step
885: -  accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size

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

891:    Level: developer

893: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
894: @*/
895: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
896: {
898:   PetscInt       ncandidates = adapt->candidates.n;
899:   PetscInt       scheme = 0;
900:   PetscReal      wlte = -1.0;
901:   PetscReal      wltea = -1.0;
902:   PetscReal      wlter = -1.0;

910:   if (next_sc) *next_sc = 0;

912:   /* Do not mess with adaptivity while handling events*/
913:   if (ts->event && ts->event->status != TSEVENT_NONE) {
914:     *next_h = h;
915:     *accept = PETSC_TRUE;
916:     return(0);
917:   }

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

924:   if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
925:     /* Increase/reduce step size if end time of next step is close to or overshoots max time */
926:     PetscReal t = ts->ptime + ts->time_step, h = *next_h;
927:     PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t;
928:     PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]);
929:     PetscReal b = adapt->matchstepfac[1];
930:     if (t < tmax && tend > tmax) *next_h = hmax;
931:     if (t < tmax && tend < tmax && h*b > hmax) *next_h = hmax/2;
932:     if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax;
933:   }

935:   if (adapt->monitor) {
936:     const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : "";
937:     PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
938:     if (wlte < 0) {
939:       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);
940:     } else {
941:       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);
942:     }
943:     PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
944:   }
945:   return(0);
946: }

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

952:    Logicially Collective on TSAdapt

954:    Input Parameters:
955: +  adapt - adaptive controller context
956: -  cnt - the number of timesteps

958:    Options Database Key:
959: .  -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase

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

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

966:    Level: advanced

968: .seealso:
969: @*/
970: PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt,PetscInt cnt)
971: {
973:   adapt->timestepjustdecreased_delay = cnt;
974:   return(0);
975: }

977: /*@
978:    TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible)

980:    Collective on TSAdapt

982:    Input Parameters:
983: +  adapt - adaptive controller context
984: .  ts - time stepper
985: .  t - Current simulation time
986: -  Y - Current solution vector

988:    Output Parameter:
989: .  accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject

991:    Level: developer

993: .seealso:
994: @*/
995: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept)
996: {
997:   PetscErrorCode      ierr;
998:   SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;


1005:   if (ts->snes) {SNESGetConvergedReason(ts->snes,&snesreason);}
1006:   if (snesreason < 0) {
1007:     *accept = PETSC_FALSE;
1008:     if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
1009:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
1010:       PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
1011:       if (adapt->monitor) {
1012:         PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1013:         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);
1014:         PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1015:       }
1016:     }
1017:   } else {
1018:     *accept = PETSC_TRUE;
1019:     TSFunctionDomainError(ts,t,Y,accept);
1020:     if (*accept && adapt->checkstage) {
1021:       (*adapt->checkstage)(adapt,ts,t,Y,accept);
1022:       if (!*accept) {
1023:         PetscInfo1(ts,"Step=%D, solution rejected by user function provided by TSSetFunctionDomainError()\n",ts->steps);
1024:         if (adapt->monitor) {
1025:           PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1026:           PetscViewerASCIIPrintf(adapt->monitor,"    TSAdapt %s step %3D stage rejected by user function provided by TSSetFunctionDomainError()\n",((PetscObject)adapt)->type_name,ts->steps);
1027:           PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1028:         }
1029:       }
1030:     }
1031:   }

1033:   if (!(*accept) && !ts->reason) {
1034:     PetscReal dt,new_dt;
1035:     TSGetTimeStep(ts,&dt);
1036:     new_dt = dt * adapt->scale_solve_failed;
1037:     TSSetTimeStep(ts,new_dt);
1038:     adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay;
1039:     if (adapt->monitor) {
1040:       PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1041:       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);
1042:       PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
1043:     }
1044:   }
1045:   return(0);
1046: }

1048: /*@
1049:   TSAdaptCreate - create an adaptive controller context for time stepping

1051:   Collective

1053:   Input Parameter:
1054: . comm - The communicator

1056:   Output Parameter:
1057: . adapt - new TSAdapt object

1059:   Level: developer

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

1064: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
1065: @*/
1066: PetscErrorCode  TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
1067: {
1069:   TSAdapt        adapt;

1073:   *inadapt = NULL;
1074:   TSAdaptInitializePackage();

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

1078:   adapt->always_accept      = PETSC_FALSE;
1079:   adapt->safety             = 0.9;
1080:   adapt->reject_safety      = 0.5;
1081:   adapt->clip[0]            = 0.1;
1082:   adapt->clip[1]            = 10.;
1083:   adapt->dt_min             = 1e-20;
1084:   adapt->dt_max             = 1e+20;
1085:   adapt->ignore_max         = -1.0;
1086:   adapt->glee_use_local     = PETSC_TRUE;
1087:   adapt->scale_solve_failed = 0.25;
1088:   /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case
1089:      to prevent from situations were unreasonably small time steps are taken in order to match the final time */
1090:   adapt->matchstepfac[0]    = 0.01; /* allow 1% step size increase in the last step */
1091:   adapt->matchstepfac[1]    = 2.0;  /* halve last step if it is greater than what remains divided this factor */
1092:   adapt->wnormtype          = NORM_2;
1093:   adapt->timestepjustdecreased_delay = 0;

1095:   *inadapt = adapt;
1096:   return(0);
1097: }