1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2 #include <petscdm.h>
3
TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept,PetscReal * wlte,PetscReal * wltea,PetscReal * wlter)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*/
TSAdaptCreate_Basic(TSAdapt adapt)79 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
80 {
81 PetscFunctionBegin;
82 adapt->ops->choose = TSAdaptChoose_Basic;
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85