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;

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

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

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

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

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

 67:   *next_h = PetscClipInterval(h_lte, adapt->dt_min, adapt->dt_max);
 68:   *wlte   = enorm;
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

 72: /*MC
 73:    TSADAPTBASIC - Basic adaptive controller for time stepping

 75:    Level: intermediate

 77: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`
 78: M*/
 79: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
 80: {
 81:   PetscFunctionBegin;
 82:   adapt->ops->choose = TSAdaptChoose_Basic;
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }