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: }