Actual source code: adaptglee.c

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

      4: typedef struct {
      5:   Vec Y;
      6: } TSAdapt_GLEE;

      8: static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
      9: {
     10:   TSAdapt_GLEE *glee = (TSAdapt_GLEE *)adapt->data;
     11:   Vec           X, Y, E;
     12:   PetscReal     enorm, enorma, enormr, hfac_lte, hfac_ltea, hfac_lter, h_lte, safety;
     13:   PetscInt      order;
     14:   PetscBool     bGTEMethod;

     16:   PetscFunctionBegin;
     17:   *next_sc = 0; /* Reuse the same order scheme */
     18:   safety   = adapt->safety;
     19:   PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSGLEE, &bGTEMethod));
     20:   order = adapt->candidates.order[0];

     22:   if (bGTEMethod) { /* the method is of GLEE type */
     23:     DM dm;

     25:     PetscCall(TSGetSolution(ts, &X));
     26:     if (!glee->Y && adapt->glee_use_local) {
     27:       PetscCall(VecDuplicate(X, &glee->Y)); /*create vector to store previous step global error*/
     28:       PetscCall(VecZeroEntries(glee->Y));   /*set error to zero on the first step - may not work if error is not zero initially*/
     29:     }
     30:     PetscCall(TSGetDM(ts, &dm));
     31:     PetscCall(DMGetGlobalVector(dm, &E));
     32:     PetscCall(TSGetTimeError(ts, 0, &E));

     34:     if (adapt->glee_use_local) PetscCall(VecAXPY(E, -1.0, glee->Y)); /* local error = current error - previous step error */

     36:     /* this should be called with the solution at the beginning of the step too*/
     37:     PetscCall(TSErrorWeightedENorm(ts, E, X, X, adapt->wnormtype, &enorm, &enorma, &enormr));
     38:     PetscCall(DMRestoreGlobalVector(dm, &E));
     39:   } else {
     40:     /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */
     41:     PetscCall(TSGetSolution(ts, &X));
     42:     if (!glee->Y) PetscCall(VecDuplicate(X, &glee->Y));
     43:     Y = glee->Y;
     44:     PetscCall(TSEvaluateStep(ts, order - 1, Y, NULL));
     45:     PetscCall(TSErrorWeightedNorm(ts, X, Y, adapt->wnormtype, &enorm, &enorma, &enormr));
     46:   }

     48:   if (enorm < 0) {
     49:     *accept = PETSC_TRUE;
     50:     *next_h = h;  /* Reuse the old step */
     51:     *wlte   = -1; /* Weighted error was not evaluated */
     52:     *wltea  = -1; /* Weighted absolute error was not evaluated */
     53:     *wlter  = -1; /* Weighted relative error was not evaluated */
     54:     PetscFunctionReturn(PETSC_SUCCESS);
     55:   }

     57:   if (enorm > 1. || enorma > 1. || enormr > 1.) {
     58:     if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
     59:     if (h < (1 + PETSC_SQRT_MACHINE_EPSILON) * adapt->dt_min) {
     60:       PetscCall(PetscInfo(adapt, "Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], accepting because step size %g is at minimum\n", (double)enorm, (double)enorma, (double)enormr, (double)h));
     61:       *accept = PETSC_TRUE;
     62:     } else if (adapt->always_accept) {
     63:       PetscCall(PetscInfo(adapt, "Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], accepting step of size %g because always_accept is set\n", (double)enorm, (double)enorma, (double)enormr, (double)h));
     64:       *accept = PETSC_TRUE;
     65:     } else {
     66:       PetscCall(PetscInfo(adapt, "Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], rejecting step of size %g\n", (double)enorm, (double)enorma, (double)enormr, (double)h));
     67:       *accept = PETSC_FALSE;
     68:     }
     69:   } else {
     70:     PetscCall(PetscInfo(adapt, "Estimated scaled truncation error [combined, absolute, relative] [%g, %g, %g], accepting step of size %g\n", (double)enorm, (double)enorma, (double)enormr, (double)h));
     71:     *accept = PETSC_TRUE;
     72:   }

     74:   if (bGTEMethod) {
     75:     if (*accept == PETSC_TRUE && adapt->glee_use_local) {
     76:       /* If step is accepted, then overwrite previous step error with the current error to be used on the next step */
     77:       /* WARNING: if the adapters are composable, then the accept test will not be reliable*/
     78:       PetscCall(TSGetTimeError(ts, 0, &glee->Y));
     79:     }

     81:     /* The optimal new step based on the current global truncation error. */
     82:     if (enorm > 0) {
     83:       /* factor based on the absolute tolerance */
     84:       hfac_ltea = safety * PetscPowReal(1. / enorma, ((PetscReal)1) / (order + 1));
     85:       /* factor based on the relative tolerance */
     86:       hfac_lter = safety * PetscPowReal(1. / enormr, ((PetscReal)1) / (order + 1));
     87:       /* pick the minimum time step among the relative and absolute tolerances */
     88:       hfac_lte = PetscMin(hfac_ltea, hfac_lter);
     89:     } else {
     90:       hfac_lte = safety * PETSC_INFINITY;
     91:     }
     92:     h_lte   = h * PetscClipInterval(hfac_lte, adapt->clip[0], adapt->clip[1]);
     93:     *next_h = PetscClipInterval(h_lte, adapt->dt_min, adapt->dt_max);
     94:   } else {
     95:     /* The optimal new step based purely on local truncation error for this step. */
     96:     if (enorm > 0) {
     97:       /* factor based on the absolute tolerance */
     98:       hfac_ltea = safety * PetscPowReal(enorma, ((PetscReal)-1) / order);
     99:       /* factor based on the relative tolerance */
    100:       hfac_lter = safety * PetscPowReal(enormr, ((PetscReal)-1) / order);
    101:       /* pick the minimum time step among the relative and absolute tolerances */
    102:       hfac_lte = PetscMin(hfac_ltea, hfac_lter);
    103:     } else {
    104:       hfac_lte = safety * PETSC_INFINITY;
    105:     }
    106:     h_lte   = h * PetscClipInterval(hfac_lte, adapt->clip[0], adapt->clip[1]);
    107:     *next_h = PetscClipInterval(h_lte, adapt->dt_min, adapt->dt_max);
    108:   }
    109:   *wlte  = enorm;
    110:   *wltea = enorma;
    111:   *wlter = enormr;
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
    116: {
    117:   TSAdapt_GLEE *glee = (TSAdapt_GLEE *)adapt->data;

    119:   PetscFunctionBegin;
    120:   PetscCall(VecDestroy(&glee->Y));
    121:   PetscFunctionReturn(PETSC_SUCCESS);
    122: }

    124: static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt)
    125: {
    126:   PetscFunctionBegin;
    127:   PetscCall(TSAdaptReset_GLEE(adapt));
    128:   PetscCall(PetscFree(adapt->data));
    129:   PetscFunctionReturn(PETSC_SUCCESS);
    130: }

    132: /*MC
    133:    TSADAPTGLEE - GLEE adaptive controller for time stepping

    135:    Level: intermediate

    137: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`
    138: M*/
    139: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
    140: {
    141:   TSAdapt_GLEE *glee;

    143:   PetscFunctionBegin;
    144:   PetscCall(PetscNew(&glee));
    145:   adapt->data         = (void *)glee;
    146:   adapt->ops->choose  = TSAdaptChoose_GLEE;
    147:   adapt->ops->reset   = TSAdaptReset_GLEE;
    148:   adapt->ops->destroy = TSAdaptDestroy_GLEE;
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }