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: dsp->cerror[0] = dsp->hratio[0] = 1.0;
47: dsp->cerror[1] = dsp->hratio[1] = 1.0;
48: dsp->cerror[2] = dsp->hratio[2] = 1.0;
49: return 0;
50: }
52: static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
53: {
54: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
55: dsp->cerror[0] = dsp->cerror[1];
56: dsp->cerror[1] = dsp->cerror[2];
57: dsp->cerror[2] = 1.0;
58: dsp->hratio[0] = dsp->hratio[1];
59: dsp->hratio[1] = dsp->hratio[2];
60: dsp->hratio[2] = 1.0;
61: return 0;
62: }
64: static PetscErrorCode TSAdaptChoose_DSP(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
65: {
66: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
67: PetscInt order = PETSC_DECIDE;
68: PetscReal enorm = -1;
69: PetscReal enorma, enormr;
70: PetscReal safety = adapt->safety * (PetscReal)0.9;
71: PetscReal hnew, hfac = PETSC_INFINITY;
72: PetscReal hmin = adapt->dt_min * (1 + PETSC_SQRT_MACHINE_EPSILON);
74: *next_sc = 0; /* Reuse the same order scheme */
75: *wltea = -1; /* Weighted absolute local truncation error is not used */
76: *wlter = -1; /* Weighted relative local truncation error is not used */
78: if (ts->ops->evaluatewlte) {
79: TSEvaluateWLTE(ts, adapt->wnormtype, &order, &enorm);
81: } else if (ts->ops->evaluatestep) {
82: DM dm;
83: Vec Y;
87: order = adapt->candidates.order[0];
88: TSGetDM(ts, &dm);
89: DMGetGlobalVector(dm, &Y);
90: TSEvaluateStep(ts, order - 1, Y, NULL);
91: TSErrorWeightedNorm(ts, ts->vec_sol, Y, adapt->wnormtype, &enorm, &enorma, &enormr);
92: DMRestoreGlobalVector(dm, &Y);
93: }
94: if (enorm < 0) {
95: TSAdaptRestart_DSP(adapt);
96: *accept = PETSC_TRUE; /* Accept the step */
97: *next_h = h; /* Reuse the old step size */
98: *wlte = -1; /* Weighted local truncation error was not evaluated */
99: return 0;
100: }
102: PetscCitationsRegister(citation[0], &cited[0]);
103: PetscCitationsRegister(citation[1], &cited[1]);
105: /* Update history after rollback */
106: if (!ts->steprollback) dsp->rollback = PETSC_FALSE;
107: else if (!dsp->rollback) {
108: dsp->rollback = PETSC_TRUE;
109: TSAdaptRollBack_DSP(adapt);
110: }
111: /* Reset history after restart */
112: if (ts->steprestart) TSAdaptRestart_DSP(adapt);
114: {
115: PetscReal k = (PetscReal)order;
116: PetscReal b1 = dsp->kBeta[0];
117: PetscReal b2 = dsp->kBeta[1];
118: PetscReal b3 = dsp->kBeta[2];
119: PetscReal a2 = dsp->Alpha[0];
120: PetscReal a3 = dsp->Alpha[1];
122: PetscReal ctr0;
123: PetscReal ctr1 = dsp->cerror[0];
124: PetscReal ctr2 = dsp->cerror[1];
125: PetscReal rho0;
126: PetscReal rho1 = dsp->hratio[0];
127: PetscReal rho2 = dsp->hratio[1];
129: /* Compute the step size ratio */
130: enorm = PetscMax(enorm, PETSC_SMALL);
131: ctr0 = PetscPowReal(1 / enorm, 1 / k);
132: rho0 = PetscPowReal(ctr0, b1);
133: rho0 *= PetscPowReal(ctr1, b2);
134: rho0 *= PetscPowReal(ctr2, b3);
135: rho0 *= PetscPowReal(rho1, -a2);
136: rho0 *= PetscPowReal(rho2, -a3);
137: rho0 = Limiter(rho0, 1);
139: /* Determine whether the step is accepted or rejected */
140: if (rho0 >= safety) *accept = PETSC_TRUE;
141: else if (adapt->always_accept) *accept = PETSC_TRUE;
142: else if (h < hmin) *accept = PETSC_TRUE;
143: else *accept = PETSC_FALSE;
145: /* Update history after accept */
146: if (*accept) {
147: dsp->cerror[2] = dsp->cerror[1];
148: dsp->cerror[1] = dsp->cerror[0];
149: dsp->cerror[0] = ctr0;
150: dsp->hratio[2] = dsp->hratio[1];
151: dsp->hratio[1] = dsp->hratio[0];
152: dsp->hratio[0] = rho0;
153: dsp->rollback = PETSC_FALSE;
154: }
156: hfac = rho0;
157: }
159: hnew = h * PetscClipInterval(hfac, adapt->clip[0], adapt->clip[1]);
160: *next_h = PetscClipInterval(hnew, adapt->dt_min, adapt->dt_max);
161: *wlte = enorm;
162: return 0;
163: }
165: static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt)
166: {
167: PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", NULL);
168: PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", NULL);
169: PetscFree(adapt->data);
170: return 0;
171: }
173: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt, PetscViewer viewer)
174: {
175: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
176: PetscBool iascii;
178: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
179: if (iascii) {
180: double a2 = (double)dsp->Alpha[0], a3 = (double)dsp->Alpha[1];
181: double b1 = (double)dsp->kBeta[0], b2 = (double)dsp->kBeta[1], b3 = (double)dsp->kBeta[2];
182: PetscViewerASCIIPrintf(viewer, "filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n", b1, b2, b3, a2, a3);
183: }
184: return 0;
185: }
187: struct FilterTab {
188: const char *name;
189: PetscReal scale;
190: PetscReal kBeta[3];
191: PetscReal Alpha[2];
192: };
194: static struct FilterTab filterlist[] = {
195: {"basic", 1, {1, 0, 0}, {0, 0} },
197: {"PI30", 3, {1, 0, 0}, {0, 0} },
198: {"PI42", 5, {3, -1, 0}, {0, 0} },
199: {"PI33", 3, {2, -1, 0}, {0, 0} },
200: {"PI34", 10, {7, -4, 0}, {0, 0} },
202: {"PC11", 1, {2, -1, 0}, {-1, 0} },
203: {"PC47", 10, {11, -7, 0}, {-10, 0} },
204: {"PC36", 10, {9, -6, 0}, {-10, 0} },
206: {"H0211", 2, {1, 1, 0}, {1, 0} },
207: {"H211b", 4, {1, 1, 0}, {1, 0} },
208: {"H211PI", 6, {1, 1, 0}, {0, 0} },
210: {"H0312", 4, {1, 2, 1}, {3, 1} },
211: {"H312b", 8, {1, 2, 1}, {3, 1} },
212: {"H312PID", 18, {1, 2, 1}, {0, 0} },
214: {"H0321", 4, {5, 2, -3}, {-1, -3} },
215: {"H321", 18, {6, 1, -5}, {-15, -3}},
216: };
218: static PetscErrorCode TSAdaptDSPSetFilter_DSP(TSAdapt adapt, const char *name)
219: {
220: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
221: PetscInt i, count = (PetscInt)(sizeof(filterlist) / sizeof(filterlist[0]));
222: struct FilterTab *tab = NULL;
223: PetscBool match;
225: for (i = 0; i < count; i++) {
226: PetscStrcasecmp(name, filterlist[i].name, &match);
227: if (match) {
228: tab = &filterlist[i];
229: break;
230: }
231: }
233: dsp->kBeta[0] = tab->kBeta[0] / tab->scale;
234: dsp->kBeta[1] = tab->kBeta[1] / tab->scale;
235: dsp->kBeta[2] = tab->kBeta[2] / tab->scale;
236: dsp->Alpha[0] = tab->Alpha[0] / tab->scale;
237: dsp->Alpha[1] = tab->Alpha[1] / tab->scale;
238: return 0;
239: }
241: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
242: {
243: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
245: dsp->kBeta[0] = kkI + kkP + kkD;
246: dsp->kBeta[1] = -(kkP + 2 * kkD);
247: dsp->kBeta[2] = kkD;
248: dsp->Alpha[0] = 0;
249: dsp->Alpha[1] = 0;
250: return 0;
251: }
253: static PetscErrorCode TSAdaptSetFromOptions_DSP(TSAdapt adapt, PetscOptionItems *PetscOptionsObject)
254: {
255: TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
256: const char *names[sizeof(filterlist) / sizeof(filterlist[0])];
257: PetscInt count = (PetscInt)(sizeof(filterlist) / sizeof(filterlist[0]));
258: PetscInt index = 2; /* PI42 */
259: PetscReal pid[3] = {1, 0, 0};
260: PetscInt i, n;
261: PetscBool set;
263: for (i = 0; i < count; i++) names[i] = filterlist[i].name;
264: PetscOptionsHeadBegin(PetscOptionsObject, "DSP adaptive controller options");
266: PetscOptionsEList("-ts_adapt_dsp_filter", "Filter name", "TSAdaptDSPSetFilter", names, count, names[index], &index, &set);
267: if (set) TSAdaptDSPSetFilter(adapt, names[index]);
269: PetscOptionsRealArray("-ts_adapt_dsp_pid", "PID parameters <kkI,kkP,kkD>", "TSAdaptDSPSetPID", pid, (n = 3, &n), &set);
271: if (set) TSAdaptDSPSetPID(adapt, pid[0], pid[1], pid[2]);
273: PetscOptionsRealArray("-ts_adapt_dsp_kbeta", "Filter parameters", "", dsp->kBeta, (n = 3, &n), &set);
275: if (set)
276: for (i = n; i < 3; i++) dsp->kBeta[i] = 0;
278: PetscOptionsRealArray("-ts_adapt_dsp_alpha", "Filter parameters", "", dsp->Alpha, (n = 2, &n), &set);
280: if (set)
281: for (i = n; i < 2; i++) dsp->Alpha[i] = 0;
283: PetscOptionsHeadEnd();
284: return 0;
285: }
287: /*@C
288: TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter
290: Collective
292: Input Parameters:
293: + adapt - adaptive controller context
294: - name - filter name
296: Options Database Key:
297: . -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
299: Filter names:
300: + basic - similar to `TSADAPTBASIC` but with different criteria for step rejections.
301: . PI30, PI42, PI33, PI34 - PI controllers.
302: . PC11, PC47, PC36 - predictive controllers.
303: . H0211, H211b, H211PI - digital filters with orders dynamics=2, adaptivity=1, filter=1.
304: . H0312, H312b, H312PID - digital filters with orders dynamics=3, adaptivity=1, filter=2.
305: - H0321, H321 - digital filters with orders dynamics=3, adaptivity=2, filter=1.
307: Level: intermediate
309: References:
310: . * - http://dx.doi.org/10.1145/641876.641877
312: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`
313: @*/
314: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt, const char *name)
315: {
318: PetscTryMethod(adapt, "TSAdaptDSPSetFilter_C", (TSAdapt, const char *), (adapt, name));
319: return 0;
320: }
322: /*@
323: TSAdaptDSPSetPID - Set the PID controller parameters
325: Input Parameters:
326: + adapt - adaptive controller context
327: . kkI - Integral parameter
328: . kkP - Proportional parameter
329: - kkD - Derivative parameter
331: Options Database Key:
332: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
334: Level: intermediate
336: References:
337: . * - http://dx.doi.org/10.1016/j.cam.2005.03.008
339: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetFilter()`
340: @*/
341: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
342: {
347: PetscTryMethod(adapt, "TSAdaptDSPSetPID_C", (TSAdapt, PetscReal, PetscReal, PetscReal), (adapt, kkI, kkP, kkD));
348: return 0;
349: }
351: /*MC
352: TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)
354: Options Database Key:
355: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
356: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
357: . -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters
358: - -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters
360: Level: intermediate
362: References:
363: + * - http://dx.doi.org/10.1145/641876.641877
364: - * - http://dx.doi.org/10.1016/j.cam.2005.03.008
366: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`, `TSAdaptDSPSetFilter()`, `TSAdaptType`
367: M*/
368: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
369: {
370: TSAdapt_DSP *dsp;
372: PetscNew(&dsp);
373: adapt->reject_safety = 1.0; /* unused */
375: adapt->data = (void *)dsp;
376: adapt->ops->choose = TSAdaptChoose_DSP;
377: adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
378: adapt->ops->destroy = TSAdaptDestroy_DSP;
379: adapt->ops->view = TSAdaptView_DSP;
381: PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", TSAdaptDSPSetFilter_DSP);
382: PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", TSAdaptDSPSetPID_DSP);
384: TSAdaptDSPSetFilter_DSP(adapt, "PI42");
385: TSAdaptRestart_DSP(adapt);
386: return 0;
387: }