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 { 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; 13 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 */ 18 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 } 32 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 } 39 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 } 57 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]); 66 67 *next_h = PetscClipInterval(h_lte, adapt->dt_min, adapt->dt_max); 68 *wlte = enorm; 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 72 /*MC 73 TSADAPTBASIC - Basic adaptive controller for time stepping 74 75 Level: intermediate 76 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 } 85