1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdm.h> 3 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 Vec Y; 6 DM dm; 7 PetscInt order = PETSC_DECIDE; 8 PetscReal enorm = -1; 9 PetscReal enorma, enormr; 10 PetscReal safety = adapt->safety; 11 PetscReal hfac_lte, h_lte; 12 13 PetscFunctionBegin; 14 *next_sc = 0; /* Reuse the same order scheme */ 15 *wltea = -1; /* Weighted absolute local truncation error is not used */ 16 *wlter = -1; /* Weighted relative local truncation error is not used */ 17 18 if (ts->ops->evaluatewlte) { 19 PetscCall(TSEvaluateWLTE(ts, adapt->wnormtype, &order, &enorm)); 20 PetscCheck(enorm < 0 || order >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Computed error order %" PetscInt_FMT " must be positive", order); 21 } else if (ts->ops->evaluatestep) { 22 PetscCheck(adapt->candidates.n >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONGSTATE, "No candidate has been registered"); 23 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); 24 order = adapt->candidates.order[0]; 25 PetscCall(TSGetDM(ts, &dm)); 26 PetscCall(DMGetGlobalVector(dm, &Y)); 27 PetscCall(TSEvaluateStep(ts, order - 1, Y, NULL)); 28 PetscCall(TSErrorWeightedNorm(ts, ts->vec_sol, Y, adapt->wnormtype, &enorm, &enorma, &enormr)); 29 PetscCall(DMRestoreGlobalVector(dm, &Y)); 30 } 31 32 if (enorm < 0) { 33 *accept = PETSC_TRUE; 34 *next_h = h; /* Reuse the old step */ 35 *wlte = -1; /* Weighted local truncation error was not evaluated */ 36 PetscFunctionReturn(0); 37 } 38 39 /* Determine whether the step is accepted of rejected */ 40 if (enorm > 1) { 41 if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */ 42 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON) * adapt->dt_min) { 43 PetscCall(PetscInfo(adapt, "Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n", (double)enorm, (double)h)); 44 *accept = PETSC_TRUE; 45 } else if (adapt->always_accept) { 46 PetscCall(PetscInfo(adapt, "Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n", (double)enorm, (double)h)); 47 *accept = PETSC_TRUE; 48 } else { 49 PetscCall(PetscInfo(adapt, "Estimated scaled local truncation error %g, rejecting step of size %g\n", (double)enorm, (double)h)); 50 *accept = PETSC_FALSE; 51 } 52 } else { 53 PetscCall(PetscInfo(adapt, "Estimated scaled local truncation error %g, accepting step of size %g\n", (double)enorm, (double)h)); 54 *accept = PETSC_TRUE; 55 } 56 57 /* The optimal new step based purely on local truncation error for this step. */ 58 if (enorm > 0) hfac_lte = safety * PetscPowReal(enorm, ((PetscReal)-1) / order); 59 else hfac_lte = safety * PETSC_INFINITY; 60 if (adapt->timestepjustdecreased) { 61 hfac_lte = PetscMin(hfac_lte, 1.0); 62 adapt->timestepjustdecreased--; 63 } 64 h_lte = h * PetscClipInterval(hfac_lte, adapt->clip[0], adapt->clip[1]); 65 66 *next_h = PetscClipInterval(h_lte, adapt->dt_min, adapt->dt_max); 67 *wlte = enorm; 68 PetscFunctionReturn(0); 69 } 70 71 /*MC 72 TSADAPTBASIC - Basic adaptive controller for time stepping 73 74 Level: intermediate 75 76 .seealso: `TS`, `TSAdapt`, `TSGetAdapt()` 77 M*/ 78 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) { 79 PetscFunctionBegin; 80 adapt->ops->choose = TSAdaptChoose_Basic; 81 PetscFunctionReturn(0); 82 } 83