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;

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

 54: static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
 55: {
 56:   TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;

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

 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);

 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 */

 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;

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

107:   PetscCall(PetscCitationsRegister(citation[0], &cited[0]));
108:   PetscCall(PetscCitationsRegister(citation[1], &cited[1]));

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));

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];

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];

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);

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;

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

161:     hfac = rho0;
162:   }

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

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

179: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt, PetscViewer viewer)
180: {
181:   TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;
182:   PetscBool    iascii;

184:   PetscFunctionBegin;
185:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
186:   if (iascii) {
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: }

194: struct FilterTab {
195:   const char *name;
196:   PetscReal   scale;
197:   PetscReal   kBeta[3];
198:   PetscReal   Alpha[2];
199: };

201: static struct FilterTab filterlist[] = {
202:   {"basic",   1,  {1, 0, 0},   {0, 0}   },

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

209:   {"PC11",    1,  {2, -1, 0},  {-1, 0}  },
210:   {"PC47",    10, {11, -7, 0}, {-10, 0} },
211:   {"PC36",    10, {9, -6, 0},  {-10, 0} },

213:   {"H0211",   2,  {1, 1, 0},   {1, 0}   },
214:   {"H211b",   4,  {1, 1, 0},   {1, 0}   },
215:   {"H211PI",  6,  {1, 1, 0},   {0, 0}   },

217:   {"H0312",   4,  {1, 2, 1},   {3, 1}   },
218:   {"H312b",   8,  {1, 2, 1},   {3, 1}   },
219:   {"H312PID", 18, {1, 2, 1},   {0, 0}   },

221:   {"H0321",   4,  {5, 2, -3},  {-1, -3} },
222:   {"H321",    18, {6, 1, -5},  {-15, -3}},
223: };

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;

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

249: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
250: {
251:   TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data;

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

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;

272:   PetscFunctionBegin;
273:   for (i = 0; i < count; i++) names[i] = filterlist[i].name;
274:   PetscOptionsHeadBegin(PetscOptionsObject, "DSP adaptive controller options");

276:   PetscCall(PetscOptionsEList("-ts_adapt_dsp_filter", "Filter name", "TSAdaptDSPSetFilter", names, count, names[index], &index, &set));
277:   if (set) PetscCall(TSAdaptDSPSetFilter(adapt, names[index]));

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]));

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;

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;

293:   PetscOptionsHeadEnd();
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: /*@
298:   TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter {cite}`soderlind2006adaptive` {cite}`soderlind2003digital`

300:   Collective

302:   Input Parameters:
303: + adapt - adaptive controller context
304: - name  - filter name

306:   Options Database Key:
307: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers

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.

317:   Level: intermediate

319: .seealso: [](ch_ts), `TSADAPTDSP`, `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`
320: @*/
321: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt, const char name[])
322: {
323:   PetscFunctionBegin;
325:   PetscAssertPointer(name, 2);
326:   PetscTryMethod(adapt, "TSAdaptDSPSetFilter_C", (TSAdapt, const char *), (adapt, name));
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

330: /*@
331:   TSAdaptDSPSetPID - Set the PID controller parameters {cite}`soderlind2006adaptive`  {cite}`soderlind2003digital`

333:   Input Parameters:
334: + adapt - adaptive controller context
335: . kkI   - Integral parameter
336: . kkP   - Proportional parameter
337: - kkD   - Derivative parameter

339:   Options Database Key:
340: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters

342:   Level: intermediate

344: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetFilter()`
345: @*/
346: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD)
347: {
348:   PetscFunctionBegin;
353:   PetscTryMethod(adapt, "TSAdaptDSPSetPID_C", (TSAdapt, PetscReal, PetscReal, PetscReal), (adapt, kkI, kkP, kkD));
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }

357: /*MC
358:    TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)
359:    {cite}`soderlind2006adaptive`  {cite}`soderlind2003digital`

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

367:    Level: intermediate

369: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`, `TSAdaptDSPSetFilter()`, `TSAdaptType`
370: M*/
371: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
372: {
373:   TSAdapt_DSP *dsp;

375:   PetscFunctionBegin;
376:   PetscCall(PetscNew(&dsp));
377:   adapt->reject_safety = 1.0; /* unused */

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;

385:   PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", TSAdaptDSPSetFilter_DSP));
386:   PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", TSAdaptDSPSetPID_DSP));

388:   PetscCall(TSAdaptDSPSetFilter_DSP(adapt, "PI42"));
389:   PetscCall(TSAdaptRestart_DSP(adapt));
390:   PetscFunctionReturn(PETSC_SUCCESS);
391: }