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: 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];
22: if (bGTEMethod) { /* the method is of GLEE type */
23: DM dm;
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));
34: if (adapt->glee_use_local) PetscCall(VecAXPY(E, -1.0, glee->Y)); /* local error = current error - previous step error */
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: }
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: }
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: }
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: }
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: }
115: static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
116: {
117: TSAdapt_GLEE *glee = (TSAdapt_GLEE *)adapt->data;
119: PetscFunctionBegin;
120: PetscCall(VecDestroy(&glee->Y));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
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: }
132: /*MC
133: TSADAPTGLEE - GLEE adaptive controller for time stepping
135: Level: intermediate
137: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`
138: M*/
139: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
140: {
141: TSAdapt_GLEE *glee;
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: }