Actual source code: adaptdsp.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdm.h>
4: static const char *citation[] = {
5: "@article{Soderlind2003,\n"
6: " author = {S\"{o}derlind, Gustaf},\n"
7: " title = {Digital Filters in Adaptive Time-stepping},\n"
8: " journal = {ACM Transactions on Mathematical Software},\n"
9: " volume = {29},\n"
10: " number = {1},\n"
11: " pages = {1--26},\n"
12: " year = {2003},\n"
13: " issn = {0098-3500},\n"
14: " doi = {http://dx.doi.org/10.1145/641876.641877},\n"
15: "}\n",
16: "@article{Soderlind2006,\n"
17: " author = {Gustaf S\"{o}derlind and Lina Wang},\n"
18: " title = {Adaptive time-stepping and computational stability},\n"
19: " journal = {Journal of Computational and Applied Mathematics},\n"
20: " volume = {185},\n"
21: " number = {2},\n"
22: " pages = {225--243},\n"
23: " year = {2006},\n"
24: " issn = {0377-0427},\n"
25: " doi = {http://dx.doi.org/10.1016/j.cam.2005.03.008},\n"
26: "}\n",
27: };
28: static PetscBool cited[] = {PETSC_FALSE, PETSC_FALSE};
30: typedef struct {
31: PetscReal kBeta[3]; /* filter parameters */
32: PetscReal Alpha[2]; /* filter parameters */
33: PetscReal cerror[3]; /* control error (controller input) history */
34: PetscReal hratio[3]; /* stepsize ratio (controller output) history */
35: PetscBool rollback;
36: } TSAdapt_DSP;
38: static PetscReal Limiter(PetscReal value, PetscReal kappa)
39: {
40: return 1 + kappa * PetscAtanReal((value - 1) / kappa);
41: }
43: static PetscErrorCode TSAdaptRestart_DSP(TSAdapt adapt)
44: {
45: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
47: PetscFunctionBegin;
48: dsp->cerror[0] = dsp->hratio[0] = 1.0;
49: dsp->cerror[1] = dsp->hratio[1] = 1.0;
50: dsp->cerror[2] = dsp->hratio[2] = 1.0;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
55: {
56: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
58: PetscFunctionBegin;
59: dsp->cerror[0] = dsp->cerror[1];
60: dsp->cerror[1] = dsp->cerror[2];
61: dsp->cerror[2] = 1.0;
62: dsp->hratio[0] = dsp->hratio[1];
63: dsp->hratio[1] = dsp->hratio[2];
64: dsp->hratio[2] = 1.0;
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: static PetscErrorCode TSAdaptChoose_DSP(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
69: {
70: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
71: PetscInt order = PETSC_DECIDE;
72: PetscReal enorm = -1;
73: PetscReal enorma, enormr;
74: PetscReal safety = adapt->safety * (PetscReal)0.9;
75: PetscReal hnew, hfac = PETSC_INFINITY;
76: PetscReal hmin = adapt->dt_min * (1 + PETSC_SQRT_MACHINE_EPSILON);
78: PetscFunctionBegin;
79: *next_sc = 0; /* Reuse the same order scheme */
80: *wltea = -1; /* Weighted absolute local truncation error is not used */
81: *wlter = -1; /* Weighted relative local truncation error is not used */
83: if (ts->ops->evaluatewlte) {
84: PetscCall(TSEvaluateWLTE(ts, adapt->wnormtype, &order, &enorm));
85: PetscCheck(enorm < 0 || order >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Computed error order %" PetscInt_FMT " must be positive", order);
86: } else if (ts->ops->evaluatestep) {
87: DM dm;
88: Vec Y;
90: PetscCheck(adapt->candidates.n >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONGSTATE, "No candidate has been registered");
91: 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);
92: order = adapt->candidates.order[0];
93: PetscCall(TSGetDM(ts, &dm));
94: PetscCall(DMGetGlobalVector(dm, &Y));
95: PetscCall(TSEvaluateStep(ts, order - 1, Y, NULL));
96: PetscCall(TSErrorWeightedNorm(ts, ts->vec_sol, Y, adapt->wnormtype, &enorm, &enorma, &enormr));
97: PetscCall(DMRestoreGlobalVector(dm, &Y));
98: }
99: if (enorm < 0) {
100: PetscCall(TSAdaptRestart_DSP(adapt));
101: *accept = PETSC_TRUE; /* Accept the step */
102: *next_h = h; /* Reuse the old step size */
103: *wlte = -1; /* Weighted local truncation error was not evaluated */
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: PetscCall(PetscCitationsRegister(citation[0], &cited[0]));
108: PetscCall(PetscCitationsRegister(citation[1], &cited[1]));
110: /* Update history after rollback */
111: if (!ts->steprollback) dsp->rollback = PETSC_FALSE;
112: else if (!dsp->rollback) {
113: dsp->rollback = PETSC_TRUE;
114: PetscCall(TSAdaptRollBack_DSP(adapt));
115: }
116: /* Reset history after restart */
117: if (ts->steprestart) PetscCall(TSAdaptRestart_DSP(adapt));
119: {
120: PetscReal k = (PetscReal)order;
121: PetscReal b1 = dsp->kBeta[0];
122: PetscReal b2 = dsp->kBeta[1];
123: PetscReal b3 = dsp->kBeta[2];
124: PetscReal a2 = dsp->Alpha[0];
125: PetscReal a3 = dsp->Alpha[1];
127: PetscReal ctr0;
128: PetscReal ctr1 = dsp->cerror[0];
129: PetscReal ctr2 = dsp->cerror[1];
130: PetscReal rho0;
131: PetscReal rho1 = dsp->hratio[0];
132: PetscReal rho2 = dsp->hratio[1];
134: /* Compute the step size ratio */
135: enorm = PetscMax(enorm, PETSC_SMALL);
136: ctr0 = PetscPowReal(1 / enorm, 1 / k);
137: rho0 = PetscPowReal(ctr0, b1);
138: rho0 *= PetscPowReal(ctr1, b2);
139: rho0 *= PetscPowReal(ctr2, b3);
140: rho0 *= PetscPowReal(rho1, -a2);
141: rho0 *= PetscPowReal(rho2, -a3);
142: rho0 = Limiter(rho0, 1);
144: /* Determine whether the step is accepted or rejected */
145: if (rho0 >= safety) *accept = PETSC_TRUE;
146: else if (adapt->always_accept) *accept = PETSC_TRUE;
147: else if (h < hmin) *accept = PETSC_TRUE;
148: else *accept = PETSC_FALSE;
150: /* Update history after accept */
151: if (*accept) {
152: dsp->cerror[2] = dsp->cerror[1];
153: dsp->cerror[1] = dsp->cerror[0];
154: dsp->cerror[0] = ctr0;
155: dsp->hratio[2] = dsp->hratio[1];
156: dsp->hratio[1] = dsp->hratio[0];
157: dsp->hratio[0] = rho0;
158: dsp->rollback = PETSC_FALSE;
159: }
161: hfac = rho0;
162: }
164: hnew = h * PetscClipInterval(hfac, adapt->clip[0], adapt->clip[1]);
165: *next_h = PetscClipInterval(hnew, adapt->dt_min, adapt->dt_max);
166: *wlte = enorm;
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt)
171: {
172: PetscFunctionBegin;
173: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", NULL));
174: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", NULL));
175: PetscCall(PetscFree(adapt->data));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt, PetscViewer viewer)
180: {
181: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
182: PetscBool iascii;
184: PetscFunctionBegin;
185: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
186: if (iascii) {
187: double a2 = (double)dsp->Alpha[0], a3 = (double)dsp->Alpha[1];
188: double b1 = (double)dsp->kBeta[0], b2 = (double)dsp->kBeta[1], b3 = (double)dsp->kBeta[2];
189: PetscCall(PetscViewerASCIIPrintf(viewer, "filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n", b1, b2, b3, a2, a3));
190: }
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: struct FilterTab {
195: const char *name;
196: PetscReal scale;
197: PetscReal kBeta[3];
198: PetscReal Alpha[2];
199: };
201: static struct FilterTab filterlist[] = {
202: {"basic", 1, {1, 0, 0}, {0, 0} },
204: {"PI30", 3, {1, 0, 0}, {0, 0} },
205: {"PI42", 5, {3, -1, 0}, {0, 0} },
206: {"PI33", 3, {2, -1, 0}, {0, 0} },
207: {"PI34", 10, {7, -4, 0}, {0, 0} },
209: {"PC11", 1, {2, -1, 0}, {-1, 0} },
210: {"PC47", 10, {11, -7, 0}, {-10, 0} },
211: {"PC36", 10, {9, -6, 0}, {-10, 0} },
213: {"H0211", 2, {1, 1, 0}, {1, 0} },
214: {"H211b", 4, {1, 1, 0}, {1, 0} },
215: {"H211PI", 6, {1, 1, 0}, {0, 0} },
217: {"H0312", 4, {1, 2, 1}, {3, 1} },
218: {"H312b", 8, {1, 2, 1}, {3, 1} },
219: {"H312PID", 18, {1, 2, 1}, {0, 0} },
221: {"H0321", 4, {5, 2, -3}, {-1, -3} },
222: {"H321", 18, {6, 1, -5}, {-15, -3}},
223: };
225: static PetscErrorCode TSAdaptDSPSetFilter_DSP(TSAdapt adapt, const char *name)
226: {
227: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
228: PetscInt i, count = PETSC_STATIC_ARRAY_LENGTH(filterlist);
229: struct FilterTab *tab = NULL;
230: PetscBool match;
232: PetscFunctionBegin;
233: for (i = 0; i < count; i++) {
234: PetscCall(PetscStrcasecmp(name, filterlist[i].name, &match));
235: if (match) {
236: tab = &filterlist[i];
237: break;
238: }
239: }
240: PetscCheck(tab, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_UNKNOWN_TYPE, "Filter name %s not found", name);
241: dsp->kBeta[0] = tab->kBeta[0] / tab->scale;
242: dsp->kBeta[1] = tab->kBeta[1] / tab->scale;
243: dsp->kBeta[2] = tab->kBeta[2] / tab->scale;
244: dsp->Alpha[0] = tab->Alpha[0] / tab->scale;
245: dsp->Alpha[1] = tab->Alpha[1] / tab->scale;
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
250: {
251: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
253: PetscFunctionBegin;
254: dsp->kBeta[0] = kkI + kkP + kkD;
255: dsp->kBeta[1] = -(kkP + 2 * kkD);
256: dsp->kBeta[2] = kkD;
257: dsp->Alpha[0] = 0;
258: dsp->Alpha[1] = 0;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: static PetscErrorCode TSAdaptSetFromOptions_DSP(TSAdapt adapt, PetscOptionItems *PetscOptionsObject)
263: {
264: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
265: const char *names[PETSC_STATIC_ARRAY_LENGTH(filterlist)];
266: PetscInt count = PETSC_STATIC_ARRAY_LENGTH(filterlist);
267: PetscInt index = 2; /* PI42 */
268: PetscReal pid[3] = {1, 0, 0};
269: PetscInt i, n;
270: PetscBool set;
272: PetscFunctionBegin;
273: for (i = 0; i < count; i++) names[i] = filterlist[i].name;
274: PetscOptionsHeadBegin(PetscOptionsObject, "DSP adaptive controller options");
276: PetscCall(PetscOptionsEList("-ts_adapt_dsp_filter", "Filter name", "TSAdaptDSPSetFilter", names, count, names[index], &index, &set));
277: if (set) PetscCall(TSAdaptDSPSetFilter(adapt, names[index]));
279: PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_pid", "PID parameters <kkI,kkP,kkD>", "TSAdaptDSPSetPID", pid, (n = 3, &n), &set));
280: PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for PID parameters");
281: if (set) PetscCall(TSAdaptDSPSetPID(adapt, pid[0], pid[1], pid[2]));
283: PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_kbeta", "Filter parameters", "", dsp->kBeta, (n = 3, &n), &set));
284: PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for parameter kbeta");
285: if (set)
286: for (i = n; i < 3; i++) dsp->kBeta[i] = 0;
288: PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_alpha", "Filter parameters", "", dsp->Alpha, (n = 2, &n), &set));
289: PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for parameter alpha");
290: if (set)
291: for (i = n; i < 2; i++) dsp->Alpha[i] = 0;
293: PetscOptionsHeadEnd();
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@
298: TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
300: Collective
302: Input Parameters:
303: + adapt - adaptive controller context
304: - name - filter name
306: Options Database Key:
307: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
309: Filter names\:
310: . basic - similar to `TSADAPTBASIC` but with different criteria for step rejections.
311: . PI30, PI42, PI33, PI34 - PI controllers.
312: . PC11, PC47, PC36 - predictive controllers.
313: . H0211, H211b, H211PI - digital filters with orders dynamics=2, adaptivity=1, filter=1.
314: . H0312, H312b, H312PID - digital filters with orders dynamics=3, adaptivity=1, filter=2.
315: - H0321, H321 - digital filters with orders dynamics=3, adaptivity=2, filter=1.
317: Level: intermediate
319: .seealso: [](ch_ts), `TSADAPTDSP`, `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`
320: @*/
321: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt, const char name[])
322: {
323: PetscFunctionBegin;
325: PetscAssertPointer(name, 2);
326: PetscTryMethod(adapt, "TSAdaptDSPSetFilter_C", (TSAdapt, const char *), (adapt, name));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*@
331: TSAdaptDSPSetPID - Set the PID controller parameters {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
333: Input Parameters:
334: + adapt - adaptive controller context
335: . kkI - Integral parameter
336: . kkP - Proportional parameter
337: - kkD - Derivative parameter
339: Options Database Key:
340: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
342: Level: intermediate
344: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetFilter()`
345: @*/
346: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
347: {
348: PetscFunctionBegin;
353: PetscTryMethod(adapt, "TSAdaptDSPSetPID_C", (TSAdapt, PetscReal, PetscReal, PetscReal), (adapt, kkI, kkP, kkD));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*MC
358: TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)
359: {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
361: Options Database Keys:
362: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
363: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
364: . -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters
365: - -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters
367: Level: intermediate
369: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`, `TSAdaptDSPSetFilter()`, `TSAdaptType`
370: M*/
371: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
372: {
373: TSAdapt_DSP *dsp;
375: PetscFunctionBegin;
376: PetscCall(PetscNew(&dsp));
377: adapt->reject_safety = 1.0; /* unused */
379: adapt->data = (void *)dsp;
380: adapt->ops->choose = TSAdaptChoose_DSP;
381: adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
382: adapt->ops->destroy = TSAdaptDestroy_DSP;
383: adapt->ops->view = TSAdaptView_DSP;
385: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", TSAdaptDSPSetFilter_DSP));
386: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", TSAdaptDSPSetPID_DSP));
388: PetscCall(TSAdaptDSPSetFilter_DSP(adapt, "PI42"));
389: PetscCall(TSAdaptRestart_DSP(adapt));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }