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;
46: PetscFunctionBegin;
47: dsp->cerror[0] = dsp->hratio[0] = 1.0;
48: dsp->cerror[1] = dsp->hratio[1] = 1.0;
49: dsp->cerror[2] = dsp->hratio[2] = 1.0;
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
54: {
55: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
56: PetscFunctionBegin;
57: dsp->cerror[0] = dsp->cerror[1];
58: dsp->cerror[1] = dsp->cerror[2];
59: dsp->cerror[2] = 1.0;
60: dsp->hratio[0] = dsp->hratio[1];
61: dsp->hratio[1] = dsp->hratio[2];
62: dsp->hratio[2] = 1.0;
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode TSAdaptChoose_DSP(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
67: {
68: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
69: PetscInt order = PETSC_DECIDE;
70: PetscReal enorm = -1;
71: PetscReal enorma, enormr;
72: PetscReal safety = adapt->safety * (PetscReal)0.9;
73: PetscReal hnew, hfac = PETSC_INFINITY;
74: PetscReal hmin = adapt->dt_min * (1 + PETSC_SQRT_MACHINE_EPSILON);
76: PetscFunctionBegin;
77: *next_sc = 0; /* Reuse the same order scheme */
78: *wltea = -1; /* Weighted absolute local truncation error is not used */
79: *wlter = -1; /* Weighted relative local truncation error is not used */
81: if (ts->ops->evaluatewlte) {
82: PetscCall(TSEvaluateWLTE(ts, adapt->wnormtype, &order, &enorm));
83: PetscCheck(enorm < 0 || order >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Computed error order %" PetscInt_FMT " must be positive", order);
84: } else if (ts->ops->evaluatestep) {
85: DM dm;
86: Vec Y;
88: PetscCheck(adapt->candidates.n >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONGSTATE, "No candidate has been registered");
89: 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);
90: order = adapt->candidates.order[0];
91: PetscCall(TSGetDM(ts, &dm));
92: PetscCall(DMGetGlobalVector(dm, &Y));
93: PetscCall(TSEvaluateStep(ts, order - 1, Y, NULL));
94: PetscCall(TSErrorWeightedNorm(ts, ts->vec_sol, Y, adapt->wnormtype, &enorm, &enorma, &enormr));
95: PetscCall(DMRestoreGlobalVector(dm, &Y));
96: }
97: if (enorm < 0) {
98: PetscCall(TSAdaptRestart_DSP(adapt));
99: *accept = PETSC_TRUE; /* Accept the step */
100: *next_h = h; /* Reuse the old step size */
101: *wlte = -1; /* Weighted local truncation error was not evaluated */
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: PetscCall(PetscCitationsRegister(citation[0], &cited[0]));
106: PetscCall(PetscCitationsRegister(citation[1], &cited[1]));
108: /* Update history after rollback */
109: if (!ts->steprollback) dsp->rollback = PETSC_FALSE;
110: else if (!dsp->rollback) {
111: dsp->rollback = PETSC_TRUE;
112: PetscCall(TSAdaptRollBack_DSP(adapt));
113: }
114: /* Reset history after restart */
115: if (ts->steprestart) PetscCall(TSAdaptRestart_DSP(adapt));
117: {
118: PetscReal k = (PetscReal)order;
119: PetscReal b1 = dsp->kBeta[0];
120: PetscReal b2 = dsp->kBeta[1];
121: PetscReal b3 = dsp->kBeta[2];
122: PetscReal a2 = dsp->Alpha[0];
123: PetscReal a3 = dsp->Alpha[1];
125: PetscReal ctr0;
126: PetscReal ctr1 = dsp->cerror[0];
127: PetscReal ctr2 = dsp->cerror[1];
128: PetscReal rho0;
129: PetscReal rho1 = dsp->hratio[0];
130: PetscReal rho2 = dsp->hratio[1];
132: /* Compute the step size ratio */
133: enorm = PetscMax(enorm, PETSC_SMALL);
134: ctr0 = PetscPowReal(1 / enorm, 1 / k);
135: rho0 = PetscPowReal(ctr0, b1);
136: rho0 *= PetscPowReal(ctr1, b2);
137: rho0 *= PetscPowReal(ctr2, b3);
138: rho0 *= PetscPowReal(rho1, -a2);
139: rho0 *= PetscPowReal(rho2, -a3);
140: rho0 = Limiter(rho0, 1);
142: /* Determine whether the step is accepted or rejected */
143: if (rho0 >= safety) *accept = PETSC_TRUE;
144: else if (adapt->always_accept) *accept = PETSC_TRUE;
145: else if (h < hmin) *accept = PETSC_TRUE;
146: else *accept = PETSC_FALSE;
148: /* Update history after accept */
149: if (*accept) {
150: dsp->cerror[2] = dsp->cerror[1];
151: dsp->cerror[1] = dsp->cerror[0];
152: dsp->cerror[0] = ctr0;
153: dsp->hratio[2] = dsp->hratio[1];
154: dsp->hratio[1] = dsp->hratio[0];
155: dsp->hratio[0] = rho0;
156: dsp->rollback = PETSC_FALSE;
157: }
159: hfac = rho0;
160: }
162: hnew = h * PetscClipInterval(hfac, adapt->clip[0], adapt->clip[1]);
163: *next_h = PetscClipInterval(hnew, adapt->dt_min, adapt->dt_max);
164: *wlte = enorm;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt)
169: {
170: PetscFunctionBegin;
171: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", NULL));
172: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", NULL));
173: PetscCall(PetscFree(adapt->data));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt, PetscViewer viewer)
178: {
179: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
180: PetscBool iascii;
182: PetscFunctionBegin;
183: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
184: if (iascii) {
185: double a2 = (double)dsp->Alpha[0], a3 = (double)dsp->Alpha[1];
186: double b1 = (double)dsp->kBeta[0], b2 = (double)dsp->kBeta[1], b3 = (double)dsp->kBeta[2];
187: PetscCall(PetscViewerASCIIPrintf(viewer, "filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n", b1, b2, b3, a2, a3));
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: struct FilterTab {
193: const char *name;
194: PetscReal scale;
195: PetscReal kBeta[3];
196: PetscReal Alpha[2];
197: };
199: static struct FilterTab filterlist[] = {
200: {"basic", 1, {1, 0, 0}, {0, 0} },
202: {"PI30", 3, {1, 0, 0}, {0, 0} },
203: {"PI42", 5, {3, -1, 0}, {0, 0} },
204: {"PI33", 3, {2, -1, 0}, {0, 0} },
205: {"PI34", 10, {7, -4, 0}, {0, 0} },
207: {"PC11", 1, {2, -1, 0}, {-1, 0} },
208: {"PC47", 10, {11, -7, 0}, {-10, 0} },
209: {"PC36", 10, {9, -6, 0}, {-10, 0} },
211: {"H0211", 2, {1, 1, 0}, {1, 0} },
212: {"H211b", 4, {1, 1, 0}, {1, 0} },
213: {"H211PI", 6, {1, 1, 0}, {0, 0} },
215: {"H0312", 4, {1, 2, 1}, {3, 1} },
216: {"H312b", 8, {1, 2, 1}, {3, 1} },
217: {"H312PID", 18, {1, 2, 1}, {0, 0} },
219: {"H0321", 4, {5, 2, -3}, {-1, -3} },
220: {"H321", 18, {6, 1, -5}, {-15, -3}},
221: };
223: static PetscErrorCode TSAdaptDSPSetFilter_DSP(TSAdapt adapt, const char *name)
224: {
225: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
226: PetscInt i, count = PETSC_STATIC_ARRAY_LENGTH(filterlist);
227: struct FilterTab *tab = NULL;
228: PetscBool match;
230: PetscFunctionBegin;
231: for (i = 0; i < count; i++) {
232: PetscCall(PetscStrcasecmp(name, filterlist[i].name, &match));
233: if (match) {
234: tab = &filterlist[i];
235: break;
236: }
237: }
238: PetscCheck(tab, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_UNKNOWN_TYPE, "Filter name %s not found", name);
239: dsp->kBeta[0] = tab->kBeta[0] / tab->scale;
240: dsp->kBeta[1] = tab->kBeta[1] / tab->scale;
241: dsp->kBeta[2] = tab->kBeta[2] / tab->scale;
242: dsp->Alpha[0] = tab->Alpha[0] / tab->scale;
243: dsp->Alpha[1] = tab->Alpha[1] / tab->scale;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
248: {
249: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
251: PetscFunctionBegin;
252: dsp->kBeta[0] = kkI + kkP + kkD;
253: dsp->kBeta[1] = -(kkP + 2 * kkD);
254: dsp->kBeta[2] = kkD;
255: dsp->Alpha[0] = 0;
256: dsp->Alpha[1] = 0;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode TSAdaptSetFromOptions_DSP(TSAdapt adapt, PetscOptionItems *PetscOptionsObject)
261: {
262: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
263: const char *names[PETSC_STATIC_ARRAY_LENGTH(filterlist)];
264: PetscInt count = PETSC_STATIC_ARRAY_LENGTH(filterlist);
265: PetscInt index = 2; /* PI42 */
266: PetscReal pid[3] = {1, 0, 0};
267: PetscInt i, n;
268: PetscBool set;
270: PetscFunctionBegin;
271: for (i = 0; i < count; i++) names[i] = filterlist[i].name;
272: PetscOptionsHeadBegin(PetscOptionsObject, "DSP adaptive controller options");
274: PetscCall(PetscOptionsEList("-ts_adapt_dsp_filter", "Filter name", "TSAdaptDSPSetFilter", names, count, names[index], &index, &set));
275: if (set) PetscCall(TSAdaptDSPSetFilter(adapt, names[index]));
277: PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_pid", "PID parameters <kkI,kkP,kkD>", "TSAdaptDSPSetPID", pid, (n = 3, &n), &set));
278: PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for PID parameters");
279: if (set) PetscCall(TSAdaptDSPSetPID(adapt, pid[0], pid[1], pid[2]));
281: PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_kbeta", "Filter parameters", "", dsp->kBeta, (n = 3, &n), &set));
282: PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for parameter kbeta");
283: if (set)
284: for (i = n; i < 3; i++) dsp->kBeta[i] = 0;
286: PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_alpha", "Filter parameters", "", dsp->Alpha, (n = 2, &n), &set));
287: PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for parameter alpha");
288: if (set)
289: for (i = n; i < 2; i++) dsp->Alpha[i] = 0;
291: PetscOptionsHeadEnd();
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@C
296: TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
298: Collective
300: Input Parameters:
301: + adapt - adaptive controller context
302: - name - filter name
304: Options Database Key:
305: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
307: Filter names\:
308: . basic - similar to `TSADAPTBASIC` but with different criteria for step rejections.
309: . PI30, PI42, PI33, PI34 - PI controllers.
310: . PC11, PC47, PC36 - predictive controllers.
311: . H0211, H211b, H211PI - digital filters with orders dynamics=2, adaptivity=1, filter=1.
312: . H0312, H312b, H312PID - digital filters with orders dynamics=3, adaptivity=1, filter=2.
313: - H0321, H321 - digital filters with orders dynamics=3, adaptivity=2, filter=1.
315: Level: intermediate
317: .seealso: [](ch_ts), `TSADAPTDSP`, `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`
318: @*/
319: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt, const char *name)
320: {
321: PetscFunctionBegin;
323: PetscAssertPointer(name, 2);
324: PetscTryMethod(adapt, "TSAdaptDSPSetFilter_C", (TSAdapt, const char *), (adapt, name));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: TSAdaptDSPSetPID - Set the PID controller parameters {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
331: Input Parameters:
332: + adapt - adaptive controller context
333: . kkI - Integral parameter
334: . kkP - Proportional parameter
335: - kkD - Derivative parameter
337: Options Database Key:
338: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
340: Level: intermediate
342: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetFilter()`
343: @*/
344: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
345: {
346: PetscFunctionBegin;
351: PetscTryMethod(adapt, "TSAdaptDSPSetPID_C", (TSAdapt, PetscReal, PetscReal, PetscReal), (adapt, kkI, kkP, kkD));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*MC
356: TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)
357: {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
359: Options Database Keys:
360: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
361: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
362: . -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters
363: - -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters
365: Level: intermediate
367: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`, `TSAdaptDSPSetFilter()`, `TSAdaptType`
368: M*/
369: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
370: {
371: TSAdapt_DSP *dsp;
373: PetscFunctionBegin;
374: PetscCall(PetscNew(&dsp));
375: adapt->reject_safety = 1.0; /* unused */
377: adapt->data = (void *)dsp;
378: adapt->ops->choose = TSAdaptChoose_DSP;
379: adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
380: adapt->ops->destroy = TSAdaptDestroy_DSP;
381: adapt->ops->view = TSAdaptView_DSP;
383: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", TSAdaptDSPSetFilter_DSP));
384: PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", TSAdaptDSPSetPID_DSP));
386: PetscCall(TSAdaptDSPSetFilter_DSP(adapt, "PI42"));
387: PetscCall(TSAdaptRestart_DSP(adapt));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }