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:   *next_sc = 0; /* Reuse the same order scheme */
 17:   safety   = adapt->safety;
 18:   PetscObjectTypeCompare((PetscObject)ts, TSGLEE, &bGTEMethod);
 19:   order = adapt->candidates.order[0];

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

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

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

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

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

 56:   if (enorm > 1. || enorma > 1. || enormr > 1.) {
 57:     if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
 58:     if (h < (1 + PETSC_SQRT_MACHINE_EPSILON) * adapt->dt_min) {
 59:       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);
 60:       *accept = PETSC_TRUE;
 61:     } else if (adapt->always_accept) {
 62:       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);
 63:       *accept = PETSC_TRUE;
 64:     } else {
 65:       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);
 66:       *accept = PETSC_FALSE;
 67:     }
 68:   } else {
 69:     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);
 70:     *accept = PETSC_TRUE;
 71:   }

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

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

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

118:   VecDestroy(&glee->Y);
119:   return 0;
120: }

122: static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt)
123: {
124:   TSAdaptReset_GLEE(adapt);
125:   PetscFree(adapt->data);
126:   return 0;
127: }

129: /*MC
130:    TSADAPTGLEE - GLEE adaptive controller for time stepping

132:    Level: intermediate

134: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`
135: M*/
136: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
137: {
138:   TSAdapt_GLEE *glee;

140:   PetscNew(&glee);
141:   adapt->data         = (void *)glee;
142:   adapt->ops->choose  = TSAdaptChoose_GLEE;
143:   adapt->ops->reset   = TSAdaptReset_GLEE;
144:   adapt->ops->destroy = TSAdaptDestroy_GLEE;
145:   return 0;
146: }