1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2 #include <petscdm.h>
3
4 typedef struct {
5 Vec Y;
6 } TSAdapt_GLEE;
7
TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept,PetscReal * wlte,PetscReal * wltea,PetscReal * wlter)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;
15
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];
21
22 if (bGTEMethod) { /* the method is of GLEE type */
23 DM dm;
24
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));
33
34 if (adapt->glee_use_local) PetscCall(VecAXPY(E, -1.0, glee->Y)); /* local error = current error - previous step error */
35
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 }
47
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 }
56
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 }
73
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 }
80
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 }
114
TSAdaptReset_GLEE(TSAdapt adapt)115 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
116 {
117 TSAdapt_GLEE *glee = (TSAdapt_GLEE *)adapt->data;
118
119 PetscFunctionBegin;
120 PetscCall(VecDestroy(&glee->Y));
121 PetscFunctionReturn(PETSC_SUCCESS);
122 }
123
TSAdaptDestroy_GLEE(TSAdapt adapt)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 }
131
132 /*MC
133 TSADAPTGLEE - GLEE adaptive controller for time stepping
134
135 Level: intermediate
136
137 .seealso: [](ch_ts), [](sec_ts_error_control), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`
138 M*/
TSAdaptCreate_GLEE(TSAdapt adapt)139 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
140 {
141 TSAdapt_GLEE *glee;
142
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 }
151