Actual source code: adaptbasic.c

      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), [](sec_ts_error_control), `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: }