1: #include <petsc/private/tsimpl.h>
  2: #include <petscdm.h>

  4: static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter)
  5: {
  6:   Vec            Y;
  7:   DM             dm;
  8:   PetscInt       order  = PETSC_DECIDE;
  9:   PetscReal      enorm  = -1;
 10:   PetscReal      enorma,enormr;
 11:   PetscReal      safety = adapt->safety;
 12:   PetscReal      hfac_lte,h_lte;

 16:   *next_sc = 0;   /* Reuse the same order scheme */
 17:   *wltea   = -1;  /* Weighted absolute local truncation error is not used */
 18:   *wlter   = -1;  /* Weighted relative local truncation error is not used */

 20:   if (ts->ops->evaluatewlte) {
 21:     TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);
 22:     if (enorm >= 0 && order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order);
 23:   } else if (ts->ops->evaluatestep) {
 24:     if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered");
 25:     if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n);
 26:     order = adapt->candidates.order[0];
 27:     TSGetDM(ts,&dm);
 28:     DMGetGlobalVector(dm,&Y);
 29:     TSEvaluateStep(ts,order-1,Y,NULL);
 30:     TSErrorWeightedNorm(ts,ts->vec_sol,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
 31:     DMRestoreGlobalVector(dm,&Y);
 32:   }

 34:   if (enorm < 0) {
 35:     *accept  = PETSC_TRUE;
 36:     *next_h  = h;            /* Reuse the old step */
 37:     *wlte    = -1;           /* Weighted local truncation error was not evaluated */
 38:     return(0);
 39:   }

 41:   /* Determine whether the step is accepted of rejected */
 42:   if (enorm > 1) {
 43:     if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
 44:     if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
 45:       PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);
 46:       *accept = PETSC_TRUE;
 47:     } else if (adapt->always_accept) {
 48:       PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);
 49:       *accept = PETSC_TRUE;
 50:     } else {
 51:       PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);
 52:       *accept = PETSC_FALSE;
 53:     }
 54:   } else {
 55:     PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);
 56:     *accept = PETSC_TRUE;
 57:   }

 59:   /* The optimal new step based purely on local truncation error for this step. */
 60:   if (enorm > 0)
 61:     hfac_lte = safety * PetscPowReal(enorm,((PetscReal)-1)/order);
 62:   else
 63:     hfac_lte = safety * PETSC_INFINITY;
 64:   if (adapt->timestepjustdecreased){
 65:     hfac_lte = PetscMin(hfac_lte,1.0);
 66:     adapt->timestepjustdecreased--;
 67:   }
 68:   h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);

 70:   *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
 71:   *wlte   = enorm;
 72:   return(0);
 73: }

 75: /*MC
 76:    TSADAPTBASIC - Basic adaptive controller for time stepping

 78:    Level: intermediate

 80: .seealso: TS, TSAdapt, TSGetAdapt()
 81: M*/
 82: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
 83: {
 85:   adapt->ops->choose = TSAdaptChoose_Basic;
 86:   return(0);
 87: }