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