Actual source code: adaptdsp.c

petsc-3.14.0 2020-09-29
Report Typos and Errors
  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:   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:   return(0);
 51: }

 53: static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
 54: {
 55:   TSAdapt_DSP *dsp = (TSAdapt_DSP*)adapt->data;
 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:   return(0);
 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);

 78:   *next_sc = 0;   /* Reuse the same order scheme */
 79:   *wltea   = -1;  /* Weighted absolute local truncation error is not used */
 80:   *wlter   = -1;  /* Weighted relative local truncation error is not used */

 82:   if (ts->ops->evaluatewlte) {
 83:     TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);
 84:     if (enorm >= 0 && order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order);
 85:   } else if (ts->ops->evaluatestep) {
 86:     DM  dm;
 87:     Vec Y;

 89:     if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered");
 90:     if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n);
 91:     order = adapt->candidates.order[0];
 92:     TSGetDM(ts,&dm);
 93:     DMGetGlobalVector(dm,&Y);
 94:     TSEvaluateStep(ts,order-1,Y,NULL);
 95:     TSErrorWeightedNorm(ts,ts->vec_sol,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
 96:     DMRestoreGlobalVector(dm,&Y);
 97:   }
 98:   if (enorm < 0) {
 99:     TSAdaptRestart_DSP(adapt);
100:     *accept = PETSC_TRUE;  /* Accept the step */
101:     *next_h = h;           /* Reuse the old step size */
102:     *wlte   = -1;          /* Weighted local truncation error was not evaluated */
103:     return(0);
104:   }

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

109:   /* Update history after rollback */
110:   if (!ts->steprollback)
111:     dsp->rollback = PETSC_FALSE;
112:   else if (!dsp->rollback) {
113:     dsp->rollback = PETSC_TRUE;
114:     TSAdaptRollBack_DSP(adapt);
115:   }
116:   /* Reset history after restart */
117:   if (ts->steprestart) {
118:     TSAdaptRestart_DSP(adapt);
119:   }

121:   {
122:     PetscReal k = (PetscReal)order;
123:     PetscReal b1 = dsp->kBeta[0];
124:     PetscReal b2 = dsp->kBeta[1];
125:     PetscReal b3 = dsp->kBeta[2];
126:     PetscReal a2 = dsp->Alpha[0];
127:     PetscReal a3 = dsp->Alpha[1];

129:     PetscReal ctr0;
130:     PetscReal ctr1 = dsp->cerror[0];
131:     PetscReal ctr2 = dsp->cerror[1];
132:     PetscReal rho0;
133:     PetscReal rho1 = dsp->hratio[0];
134:     PetscReal rho2 = dsp->hratio[1];

136:     /* Compute the step size ratio */
137:     enorm = PetscMax(enorm,PETSC_SMALL);
138:     ctr0  = PetscPowReal(1/enorm,1/k);
139:     rho0  = PetscPowReal(ctr0,b1);
140:     rho0 *= PetscPowReal(ctr1,b2);
141:     rho0 *= PetscPowReal(ctr2,b3);
142:     rho0 *= PetscPowReal(rho1,-a2);
143:     rho0 *= PetscPowReal(rho2,-a3);
144:     rho0  = Limiter(rho0,1);

146:     /* Determine whether the step is accepted or rejected */
147:     if (rho0 >= safety)
148:       *accept = PETSC_TRUE;
149:     else if (adapt->always_accept)
150:       *accept = PETSC_TRUE;
151:     else if (h < hmin)
152:       *accept = PETSC_TRUE;
153:     else
154:       *accept = PETSC_FALSE;

156:     /* Update history after accept */
157:     if (*accept) {
158:       dsp->cerror[2] = dsp->cerror[1];
159:       dsp->cerror[1] = dsp->cerror[0];
160:       dsp->cerror[0] = ctr0;
161:       dsp->hratio[2] = dsp->hratio[1];
162:       dsp->hratio[1] = dsp->hratio[0];
163:       dsp->hratio[0] = rho0;
164:       dsp->rollback  = PETSC_FALSE;
165:     }

167:     hfac = rho0;
168:   }

170:   hnew    = h * PetscClipInterval(hfac,adapt->clip[0],adapt->clip[1]);
171:   *next_h = PetscClipInterval(hnew,adapt->dt_min,adapt->dt_max);
172:   *wlte   = enorm;
173:   return(0);
174: }

176: static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt)
177: {

181:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",NULL);
182:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",NULL);
183:   PetscFree(adapt->data);
184:   return(0);
185: }

187: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt,PetscViewer viewer)
188: {
189:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;
190:   PetscBool      iascii;

194:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
195:   if (iascii) {
196:     double a2 = (double)dsp->Alpha[0], a3 = (double)dsp->Alpha[1];
197:     double b1 = (double)dsp->kBeta[0], b2 = (double)dsp->kBeta[1], b3 = (double)dsp->kBeta[2];
198:     PetscViewerASCIIPrintf(viewer,"filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n",b1,b2,b3,a2,a3);
199:   }
200:   return(0);
201: }

203: struct FilterTab {
204:   const char *name;
205:   PetscReal scale;
206:   PetscReal kBeta[3];
207:   PetscReal Alpha[2];
208: };

210: static struct FilterTab filterlist[] = {
211:   {"basic",    1, {  1,  0,  0 }, {   0,  0 }},

213:   {"PI30",     3, {  1,  0,  0 }, {   0,  0 }},
214:   {"PI42",     5, {  3, -1,  0 }, {   0,  0 }},
215:   {"PI33",     3, {  2, -1,  0 }, {   0,  0 }},
216:   {"PI34",    10, {  7, -4,  0 }, {   0,  0 }},

218:   {"PC11",     1, {  2, -1,  0 }, {  -1,  0 }},
219:   {"PC47",    10, { 11, -7,  0 }, { -10,  0 }},
220:   {"PC36",    10, {  9, -6,  0 }, { -10,  0 }},

222:   {"H0211",    2, {  1,  1,  0 }, {   1,  0 }},
223:   {"H211b",    4, {  1,  1,  0 }, {   1,  0 }},
224:   {"H211PI",   6, {  1,  1,  0 }, {   0,  0 }},

226:   {"H0312",    4, {  1,  2,  1 }, {   3,  1 }},
227:   {"H312b",    8, {  1,  2,  1 }, {   3,  1 }},
228:   {"H312PID", 18, {  1,  2,  1 }, {   0,  0 }},

230:   {"H0321",    4, {  5,  2, -3 }, {  -1, -3 }},
231:   {"H321",    18, {  6,  1, -5 }, { -15, -3 }},
232: };

234: static PetscErrorCode TSAdaptDSPSetFilter_DSP(TSAdapt adapt,const char *name)
235: {
236:   TSAdapt_DSP       *dsp = (TSAdapt_DSP*)adapt->data;
237:   PetscInt          i,count = (PetscInt)(sizeof(filterlist)/sizeof(filterlist[0]));
238:   struct FilterTab* tab = NULL;
239:   PetscBool         match;
240:   PetscErrorCode    ierr;

243:   for (i=0; i<count; i++) {
244:     PetscStrcasecmp(name,filterlist[i].name,&match);
245:     if (match) { tab = &filterlist[i]; break; }
246:   }
247:   if (!tab) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_UNKNOWN_TYPE,"Filter name %s not found",name);
248:   dsp->kBeta[0] = tab->kBeta[0]/tab->scale;
249:   dsp->kBeta[1] = tab->kBeta[1]/tab->scale;
250:   dsp->kBeta[2] = tab->kBeta[2]/tab->scale;
251:   dsp->Alpha[0] = tab->Alpha[0]/tab->scale;
252:   dsp->Alpha[1] = tab->Alpha[1]/tab->scale;
253:   return(0);
254: }

256: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)
257: {
258:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;

261:   dsp->kBeta[0] = kkI + kkP + kkD;
262:   dsp->kBeta[1] = -(kkP + 2*kkD);
263:   dsp->kBeta[2] = kkD;
264:   dsp->Alpha[0] = 0;
265:   dsp->Alpha[1] = 0;
266:   return(0);
267: }

269: static PetscErrorCode TSAdaptSetFromOptions_DSP(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
270: {
271:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;
272:   const char     *names[sizeof(filterlist)/sizeof(filterlist[0])];
273:   PetscInt       count = (PetscInt)(sizeof(filterlist)/sizeof(filterlist[0]));
274:   PetscInt       index = 2; /* PI42 */
275:   PetscReal      pid[3] = {1,0,0};
276:   PetscInt       i,n;
277:   PetscBool      set;

281:   for (i=0; i<count; i++) names[i] = filterlist[i].name;
282:   PetscOptionsHead(PetscOptionsObject,"DSP adaptive controller options");

284:   PetscOptionsEList("-ts_adapt_dsp_filter","Filter name","TSAdaptDSPSetFilter",names,count,names[index],&index,&set);
285:   if (set) { TSAdaptDSPSetFilter(adapt,names[index]);}

287:   PetscOptionsRealArray("-ts_adapt_dsp_pid","PID parameters <kkI,kkP,kkD>","TSAdaptDSPSetPID",pid,(n=3,&n),&set);
288:   if (set && !n) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONG,"Must provide at least one value for PID parameters");
289:   if (set) {TSAdaptDSPSetPID(adapt,pid[0],pid[1],pid[2]);}

291:   PetscOptionsRealArray("-ts_adapt_dsp_kbeta","Filter parameters","",dsp->kBeta,(n=3,&n),&set);
292:   if (set && !n) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONG,"Must provide at least one value for parameter kbeta");
293:   if (set) for (i=n; i<3; i++) dsp->kBeta[i] = 0;

295:   PetscOptionsRealArray("-ts_adapt_dsp_alpha","Filter parameters","",dsp->Alpha,(n=2,&n),&set);
296:   if (set && !n) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONG,"Must provide at least one value for parameter alpha");
297:   if (set) for (i=n; i<2; i++) dsp->Alpha[i] = 0;

299:   PetscOptionsTail();
300:   return(0);
301: }

303: /*@C
304:    TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter

306:    Collective on TSAdapt

308:    Input Arguments:
309: +  adapt - adaptive controller context
310: -  name - filter name

312:    Level: intermediate

314:    References: http://dx.doi.org/10.1145/641876.641877

316:    Notes:
317:     Valid filter names are
318: +  "basic" - similar to TSADAPTBASIC but with different criteria for step rejections.
319: .  "PI30", "PI42", "PI33", "PI34" - PI controllers.
320: .  "PC11", "PC47", "PC36" - predictive controllers.
321: .  "H0211", "H211b", "H211PI" - digital filters with orders dynamics=2, adaptivity=1, filter=1.
322: .  "H0312", "H312b", "H312PID" - digital filters with orders dynamics=3, adaptivity=1, filter=2.
323: -  "H0321", "H321" - digital filters with orders dynamics=3, adaptivity=2, filter=1.

325:    Options Database:
326: .   -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers

328: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID()
329: @*/
330: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt,const char *name)
331: {
336:   PetscTryMethod(adapt,"TSAdaptDSPSetFilter_C",(TSAdapt,const char*),(adapt,name));
337:   return(0);
338: }

340: /*@
341:    TSAdaptDSPSetPID - Set the PID controller parameters

343:    Input Arguments:
344: +  adapt - adaptive controller context
345: .  kkI - Integral parameter
346: .  kkP - Proportional parameter
347: -  kkD - Derivative parameter

349:    Level: intermediate

351:    References: http://dx.doi.org/10.1016/j.cam.2005.03.008

353:    Options Database:
354: .   -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters

356: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetFilter()
357: @*/
358: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)
359: {
366:   PetscTryMethod(adapt,"TSAdaptDSPSetPID_C",(TSAdapt,PetscReal,PetscReal,PetscReal),(adapt,kkI,kkP,kkD));
367:   return(0);
368: }

370: /*MC
371:    TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)

373:    Level: intermediate

375:    References: http://dx.doi.org/10.1145/641876.641877 and http://dx.doi.org/10.1016/j.cam.2005.03.008

377:    Options Database:
378: +   -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
379: .   -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
380: .   -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters
381: -   -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters

383: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID(), TSAdaptDSPSetFilter()
384: M*/
385: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
386: {
387:   TSAdapt_DSP    *dsp;

391:   PetscNewLog(adapt,&dsp);
392:   adapt->reject_safety = 1.0; /* unused */

394:   adapt->data                = (void*)dsp;
395:   adapt->ops->choose         = TSAdaptChoose_DSP;
396:   adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
397:   adapt->ops->destroy        = TSAdaptDestroy_DSP;
398:   adapt->ops->view           = TSAdaptView_DSP;

400:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",TSAdaptDSPSetFilter_DSP);
401:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",TSAdaptDSPSetPID_DSP);

403:   TSAdaptDSPSetFilter_DSP(adapt,"PI42");
404:   TSAdaptRestart_DSP(adapt);
405:   return(0);
406: }