1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2 #include <petscdm.h>
3
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};
29
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;
37
Limiter(PetscReal value,PetscReal kappa)38 static PetscReal Limiter(PetscReal value, PetscReal kappa)
39 {
40 return 1 + kappa * PetscAtanReal((value - 1) / kappa);
41 }
42
TSAdaptRestart_DSP(TSAdapt adapt)43 static PetscErrorCode TSAdaptRestart_DSP(TSAdapt adapt)
44 {
45 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
46
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 }
53
TSAdaptRollBack_DSP(TSAdapt adapt)54 static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
55 {
56 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
57
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 }
67
TSAdaptChoose_DSP(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept,PetscReal * wlte,PetscReal * wltea,PetscReal * wlter)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);
77
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 */
82
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;
89
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 }
106
107 PetscCall(PetscCitationsRegister(citation[0], &cited[0]));
108 PetscCall(PetscCitationsRegister(citation[1], &cited[1]));
109
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));
118
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];
126
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];
133
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);
143
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;
149
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 }
160
161 hfac = rho0;
162 }
163
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 }
169
TSAdaptDestroy_DSP(TSAdapt adapt)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 }
178
TSAdaptView_DSP(TSAdapt adapt,PetscViewer viewer)179 static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt, PetscViewer viewer)
180 {
181 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
182 PetscBool isascii;
183
184 PetscFunctionBegin;
185 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
186 if (isascii) {
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 }
193
194 struct FilterTab {
195 const char *name;
196 PetscReal scale;
197 PetscReal kBeta[3];
198 PetscReal Alpha[2];
199 };
200
201 static struct FilterTab filterlist[] = {
202 {"basic", 1, {1, 0, 0}, {0, 0} },
203
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} },
208
209 {"PC11", 1, {2, -1, 0}, {-1, 0} },
210 {"PC47", 10, {11, -7, 0}, {-10, 0} },
211 {"PC36", 10, {9, -6, 0}, {-10, 0} },
212
213 {"H0211", 2, {1, 1, 0}, {1, 0} },
214 {"H211b", 4, {1, 1, 0}, {1, 0} },
215 {"H211PI", 6, {1, 1, 0}, {0, 0} },
216
217 {"H0312", 4, {1, 2, 1}, {3, 1} },
218 {"H312b", 8, {1, 2, 1}, {3, 1} },
219 {"H312PID", 18, {1, 2, 1}, {0, 0} },
220
221 {"H0321", 4, {5, 2, -3}, {-1, -3} },
222 {"H321", 18, {6, 1, -5}, {-15, -3}},
223 };
224
TSAdaptDSPSetFilter_DSP(TSAdapt adapt,const char * name)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;
231
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 }
248
TSAdaptDSPSetPID_DSP(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)249 static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
250 {
251 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
252
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 }
261
TSAdaptSetFromOptions_DSP(TSAdapt adapt,PetscOptionItems PetscOptionsObject)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;
271
272 PetscFunctionBegin;
273 for (i = 0; i < count; i++) names[i] = filterlist[i].name;
274 PetscOptionsHeadBegin(PetscOptionsObject, "DSP adaptive controller options");
275
276 PetscCall(PetscOptionsEList("-ts_adapt_dsp_filter", "Filter name", "TSAdaptDSPSetFilter", names, count, names[index], &index, &set));
277 if (set) PetscCall(TSAdaptDSPSetFilter(adapt, names[index]));
278
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]));
282
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;
287
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;
292
293 PetscOptionsHeadEnd();
294 PetscFunctionReturn(PETSC_SUCCESS);
295 }
296
297 /*@
298 TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
299
300 Collective
301
302 Input Parameters:
303 + adapt - adaptive controller context
304 - name - filter name
305
306 Options Database Key:
307 + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
308
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.
316
317 Level: intermediate
318
319 .seealso: [](ch_ts), `TSADAPTDSP`, `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`
320 @*/
TSAdaptDSPSetFilter(TSAdapt adapt,const char name[])321 PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt, const char name[])
322 {
323 PetscFunctionBegin;
324 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1);
325 PetscAssertPointer(name, 2);
326 PetscTryMethod(adapt, "TSAdaptDSPSetFilter_C", (TSAdapt, const char *), (adapt, name));
327 PetscFunctionReturn(PETSC_SUCCESS);
328 }
329
330 /*@
331 TSAdaptDSPSetPID - Set the PID controller parameters {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
332
333 Input Parameters:
334 + adapt - adaptive controller context
335 . kkI - Integral parameter
336 . kkP - Proportional parameter
337 - kkD - Derivative parameter
338
339 Options Database Key:
340 . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
341
342 Level: intermediate
343
344 .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetFilter()`
345 @*/
TSAdaptDSPSetPID(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)346 PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
347 {
348 PetscFunctionBegin;
349 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1);
350 PetscValidLogicalCollectiveReal(adapt, kkI, 2);
351 PetscValidLogicalCollectiveReal(adapt, kkP, 3);
352 PetscValidLogicalCollectiveReal(adapt, kkD, 4);
353 PetscTryMethod(adapt, "TSAdaptDSPSetPID_C", (TSAdapt, PetscReal, PetscReal, PetscReal), (adapt, kkI, kkP, kkD));
354 PetscFunctionReturn(PETSC_SUCCESS);
355 }
356
357 /*MC
358 TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)
359 {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`
360
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
366
367 Level: intermediate
368
369 .seealso: [](ch_ts), [](sec_ts_error_control), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`, `TSAdaptDSPSetFilter()`, `TSAdaptType`
370 M*/
TSAdaptCreate_DSP(TSAdapt adapt)371 PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
372 {
373 TSAdapt_DSP *dsp;
374
375 PetscFunctionBegin;
376 PetscCall(PetscNew(&dsp));
377 adapt->reject_safety = 1.0; /* unused */
378
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;
384
385 PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", TSAdaptDSPSetFilter_DSP));
386 PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", TSAdaptDSPSetPID_DSP));
387
388 PetscCall(TSAdaptDSPSetFilter_DSP(adapt, "PI42"));
389 PetscCall(TSAdaptRestart_DSP(adapt));
390 PetscFunctionReturn(PETSC_SUCCESS);
391 }
392