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 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);
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(), TSAdaptSetScaleSolveFailed()
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 greater 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(), TSAdaptSetScaleSolveFailed()
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:    TSAdaptSetScaleSolveFailed - Scale step by this factor if solve fails

580:    Logically collective on TSAdapt

582:    Input Arguments:
583: +  adapt - adaptive controller context
584: -  scale - scale

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

589:    Level: intermediate

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

604: /*@
605:    TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size

607:    Not Collective

609:    Input Arguments:
610: .  adapt - adaptive controller context

612:    Ouput Arguments:
613: .  scale - scale factor

615:    Level: intermediate

617: .seealso: TSAdaptChoose(), TSAdaptSetScaleSolveFailed(), TSAdaptSetClip()
618: @*/
619: PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt,PetscReal *scale)
620: {
624:   if (scale)  *scale  = adapt->scale_solve_failed;
625:   return(0);
626: }

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

631:    Logically collective on TSAdapt

633:    Input Arguments:
634: +  adapt - time step adaptivity context, usually gotten with TSGetAdapt()
635: .  hmin - minimum time step
636: -  hmax - maximum time step

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

642:    Level: intermediate

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

653:   if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin);
654:   if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax);
655:   if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin;
656:   if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax;
657:   hmin = adapt->dt_min;
658:   hmax = adapt->dt_max;
659:   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);
660:   return(0);
661: }

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

666:    Not Collective

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

671:    Output Arguments:
672: +  hmin - minimum time step
673: -  hmax - maximum time step

675:    Level: intermediate

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

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

691: /*
692:    TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.

694:    Collective on TSAdapt

696:    Input Parameter:
697: .  adapt - the TSAdapt context

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

711:    Level: advanced

713:    Notes:
714:    This function is automatically called by TSSetFromOptions()

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

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

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

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

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

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

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

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

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

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

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

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

774: /*@
775:    TSAdaptCandidatesClear - clear any previously set candidate schemes

777:    Logically collective on TSAdapt

779:    Input Argument:
780: .  adapt - adaptive controller

782:    Level: developer

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

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

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

799:    Logically collective on TSAdapt

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

810:    Note:
811:    This routine is not available in Fortran.

813:    Level: developer

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

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

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

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

843:    Not Collective

845:    Input Arguments:
846: .  adapt - time step adaptivity context

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

855:    Level: developer

857:    Note:
858:    The current scheme is always returned in the first slot

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

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

877:    Collective on TSAdapt

879:    Input Arguments:
880: +  adapt - adaptive contoller
881: -  h - current step size

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

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

892:    Level: developer

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

911:   if (next_sc) *next_sc = 0;

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

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

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

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

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

953:    Logicially Collective on TSAdapt

955:    Input Arguments:
956: +  adapt - adaptive controller context
957: -  cnt - the number of timesteps

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

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

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

967:    Level: advanced

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

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

981:    Collective on TSAdapt

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

989:    Output Arguments:
990: .  accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject

992:    Level: developer

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


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

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

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

1052:   Collective

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

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

1060:   Level: developer

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

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

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

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

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

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