1: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/

  3: typedef struct {
  4:   PetscBool always_accept;
  5:   PetscReal clip[2];            /* admissible decrease/increase factors */
  6:   PetscReal safety;             /* safety factor relative to target error */
  7:   PetscReal reject_safety;      /* extra safety factor if the last step was rejected */
  8:   Vec       Y;
  9: } TSAdapt_Basic;

 13: static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte)
 14: {
 15:   TSAdapt_Basic  *basic = (TSAdapt_Basic*)adapt->data;
 17:   Vec            X,Y;
 18:   PetscReal      enorm,hfac_lte,h_lte,safety;
 19:   PetscInt       order,stepno;

 22:   TSGetTimeStepNumber(ts,&stepno);
 23:   TSGetSolution(ts,&X);
 24:   if (!basic->Y) {VecDuplicate(X,&basic->Y);}
 25:   Y     = basic->Y;
 26:   order = adapt->candidates.order[0];
 27:   TSEvaluateStep(ts,order-1,Y,NULL);

 29:   safety = basic->safety;
 30:   TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm);
 31:   if (enorm > 1.) {
 32:     if (!*accept) safety *= basic->reject_safety; /* The last attempt also failed, shorten more aggressively */
 33:     if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
 34:       PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);
 35:       *accept = PETSC_TRUE;
 36:     } else if (basic->always_accept) {
 37:       PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);
 38:       *accept = PETSC_TRUE;
 39:     } else {
 40:       PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);
 41:       *accept = PETSC_FALSE;
 42:     }
 43:   } else {
 44:     PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);
 45:     *accept = PETSC_TRUE;
 46:   }

 48:   /* The optimal new step based purely on local truncation error for this step. */
 49:   hfac_lte = safety * PetscPowReal(enorm,-1./order);
 50:   h_lte    = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]);

 52:   *next_sc = 0;
 53:   *next_h  = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
 54:   *wlte    = enorm;
 55:   return(0);
 56: }

 60: static PetscErrorCode TSAdaptReset_Basic(TSAdapt adapt)
 61: {
 62:   TSAdapt_Basic  *basic = (TSAdapt_Basic*)adapt->data;

 66:   VecDestroy(&basic->Y);
 67:   return(0);
 68: }

 72: static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt)
 73: {

 77:   TSAdaptReset_Basic(adapt);
 78:   PetscFree(adapt->data);
 79:   return(0);
 80: }

 84: static PetscErrorCode TSAdaptSetFromOptions_Basic(PetscOptions *PetscOptionsObject,TSAdapt adapt)
 85: {
 86:   TSAdapt_Basic  *basic = (TSAdapt_Basic*)adapt->data;
 88:   PetscInt       two;
 89:   PetscBool      set;

 92:   PetscOptionsHead(PetscOptionsObject,"Basic adaptive controller options");
 93:   two  = 2;
 94:   PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase in step size","",basic->clip,&two,&set);
 95:   if (set && (two != 2 || basic->clip[0] > basic->clip[1])) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip");
 96:   PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,NULL);
 97:   PetscOptionsReal("-ts_adapt_basic_reject_safety","Extra safety factor to apply if the last step was rejected","",basic->reject_safety,&basic->reject_safety,NULL);
 98:   PetscOptionsBool("-ts_adapt_basic_always_accept","Always accept the step regardless of whether local truncation error meets goal","",basic->always_accept,&basic->always_accept,NULL);
 99:   PetscOptionsTail();
100:   return(0);
101: }

105: static PetscErrorCode TSAdaptView_Basic(TSAdapt adapt,PetscViewer viewer)
106: {
107:   TSAdapt_Basic  *basic = (TSAdapt_Basic*)adapt->data;
109:   PetscBool      iascii;

112:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
113:   if (iascii) {
114:     if (basic->always_accept) {PetscViewerASCIIPrintf(viewer,"  Basic: always accepting steps\n");}
115:     PetscViewerASCIIPrintf(viewer,"  Basic: clip fastest decrease %g, fastest increase %g\n",(double)basic->clip[0],(double)basic->clip[1]);
116:     PetscViewerASCIIPrintf(viewer,"  Basic: safety factor %g, extra factor after step rejection %g\n",(double)basic->safety,(double)basic->reject_safety);
117:   }
118:   return(0);
119: }

123: /*MC
124:    TSADAPTBASIC - Basic adaptive controller for time stepping

126:    Level: intermediate

128: .seealso: TS, TSAdapt, TSSetAdapt()
129: M*/
130: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
131: {
133:   TSAdapt_Basic  *a;

136:   PetscNewLog(adapt,&a);
137:   adapt->data                = (void*)a;
138:   adapt->ops->choose         = TSAdaptChoose_Basic;
139:   adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic;
140:   adapt->ops->destroy        = TSAdaptDestroy_Basic;
141:   adapt->ops->view           = TSAdaptView_Basic;

143:   a->clip[0]       = 0.1;
144:   a->clip[1]       = 10.;
145:   a->safety        = 0.9;
146:   a->reject_safety = 0.5;
147:   a->always_accept = PETSC_FALSE;
148:   return(0);
149: }