"Fossies" - the Fresh Open Source Software Archive

Member "petsc-3.12.3/src/ts/adapt/impls/dsp/adaptdsp.c.html" (3 Jan 2020, 43724 Bytes) of package /linux/misc/petsc-3.12.3.tar.gz:


Caution: In this restricted "Fossies" environment the current HTML page may not be correctly presentated and may have some non-functional links. You can here alternatively try to browse the pure source code or just view or download the uninterpreted raw source code. If the rendering is insufficient you may try to find and view the page on the petsc-3.12.3.tar.gz project site itself.

petsc-3.12.3 2020-01-03
Report Typos and Errors
  1:  #include <petsc/private/tsimpl.h>

  3: static const char *citation[] = {
  4:   "@article{Soderlind2003,\n"
  5:   " author = {S\"{o}derlind, Gustaf},\n"
  6:   " title = {Digital Filters in Adaptive Time-stepping},\n"
  7:   " journal = {ACM Transactions on Mathematical Software},\n"
  8:   " volume = {29},\n"
  9:   " number = {1},\n"
 10:   " pages = {1--26},\n"
 11:   " year = {2003},\n"
 12:   " issn = {0098-3500},\n"
 13:   " doi = {http://dx.doi.org/10.1145/641876.641877},\n"
 14:   "}\n",
 15:   "@article{Soderlind2006,\n"
 16:   " author = {Gustaf S\"{o}derlind and Lina Wang},\n"
 17:   " title = {Adaptive time-stepping and computational stability},\n"
 18:   " journal = {Journal of Computational and Applied Mathematics},\n"
 19:   " volume = {185},\n"
 20:   " number = {2},\n"
 21:   " pages = {225--243},\n"
 22:   " year = {2006},\n"
 23:   " issn = {0377-0427},\n"
 24:   " doi = {http://dx.doi.org/10.1016/j.cam.2005.03.008},\n"
 25:   "}\n",
 26: };
 27: static PetscBool cited[] = {PETSC_FALSE,PETSC_FALSE};

 29: typedef struct {
 30:   Vec       Y;
 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:     if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered");
 87:     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);
 88:     if (!dsp->Y) {VecDuplicate(ts->vec_sol,&dsp->Y);}
 89:     order = adapt->candidates.order[0];
 90:     TSEvaluateStep(ts,order-1,dsp->Y,NULL);
 91:     TSErrorWeightedNorm(ts,ts->vec_sol,dsp->Y,adapt->wnormtype,&enorm,&enorma,&enormr);
 92:   }
 93:   if (enorm < 0) {
 94:     TSAdaptRestart_DSP(adapt);
 95:     *accept = PETSC_TRUE;  /* Accept the step */
 96:     *next_h = h;           /* Reuse the old step size */
 97:     *wlte   = -1;          /* Weighted local truncation error was not evaluated */
 98:     return(0);
 99:   }

101:   PetscCitationsRegister(citation[0],&cited[0]);
102:   PetscCitationsRegister(citation[1],&cited[1]);

104:   /* Update history after rollback */
105:   if (!ts->steprollback)
106:     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) {
113:     TSAdaptRestart_DSP(adapt);
114:   }

116:   {
117:     PetscReal k = (PetscReal)order;
118:     PetscReal b1 = dsp->kBeta[0];
119:     PetscReal b2 = dsp->kBeta[1];
120:     PetscReal b3 = dsp->kBeta[2];
121:     PetscReal a2 = dsp->Alpha[0];
122:     PetscReal a3 = dsp->Alpha[0];

124:     PetscReal ctr0;
125:     PetscReal ctr1 = dsp->cerror[0];
126:     PetscReal ctr2 = dsp->cerror[1];
127:     PetscReal rho0;
128:     PetscReal rho1 = dsp->hratio[0];
129:     PetscReal rho2 = dsp->hratio[1];

131:     /* Compute the step size ratio */
132:     enorm = PetscMax(enorm,PETSC_SMALL);
133:     ctr0  = PetscPowReal(1/enorm,1/k);
134:     rho0  = PetscPowReal(ctr0,b1);
135:     rho0 *= PetscPowReal(ctr1,b2);
136:     rho0 *= PetscPowReal(ctr2,b3);
137:     rho0 *= PetscPowReal(rho1,-a2);
138:     rho0 *= PetscPowReal(rho2,-a3);
139:     rho0  = Limiter(rho0,1);

141:     /* Determine whether the step is accepted or rejected */
142:     if (rho0 >= safety)
143:       *accept = PETSC_TRUE;
144:     else if (adapt->always_accept)
145:       *accept = PETSC_TRUE;
146:     else if (h < hmin)
147:       *accept = PETSC_TRUE;
148:     else
149:       *accept = PETSC_FALSE;

151:     /* Update history after accept */
152:     if (*accept) {
153:       dsp->cerror[2] = dsp->cerror[1];
154:       dsp->cerror[1] = dsp->cerror[0];
155:       dsp->cerror[0] = ctr0;
156:       dsp->hratio[2] = dsp->hratio[1];
157:       dsp->hratio[1] = dsp->hratio[0];
158:       dsp->hratio[0] = rho0;
159:       dsp->rollback  = PETSC_FALSE;
160:     }

162:     hfac = rho0;
163:   }

165:   hnew    = h * PetscClipInterval(hfac,adapt->clip[0],adapt->clip[1]);
166:   *next_h = PetscClipInterval(hnew,adapt->dt_min,adapt->dt_max);
167:   *wlte   = enorm;
168:   return(0);
169: }

171: static PetscErrorCode TSAdaptReset_DSP(TSAdapt adapt)
172: {
173:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;

177:   VecDestroy(&dsp->Y);
178:   return(0);
179: }

181: static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt)
182: {

186:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",NULL);
187:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",NULL);
188:   TSAdaptReset_DSP(adapt);
189:   PetscFree(adapt->data);
190:   return(0);
191: }

193: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt,PetscViewer viewer)
194: {
195:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;
196:   PetscBool      iascii;

200:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
201:   if (iascii) {
202:     double a2 = (double)dsp->Alpha[0], a3 = (double)dsp->Alpha[1];
203:     double b1 = (double)dsp->kBeta[0], b2 = (double)dsp->kBeta[1], b3 = (double)dsp->kBeta[2];
204:     PetscViewerASCIIPrintf(viewer,"filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n",b1,b2,b3,a2,a3);
205:   }
206:   return(0);
207: }

209: struct FilterTab {
210:   const char *name;
211:   PetscReal scale;
212:   PetscReal kBeta[3];
213:   PetscReal Alpha[2];
214: };

216: static struct FilterTab filterlist[] = {
217:   {"basic",    1, {  1,  0,  0 }, {   0,  0 }},

219:   {"PI30",     3, {  1,  0,  0 }, {   0,  0 }},
220:   {"PI42",     5, {  3, -1,  0 }, {   0,  0 }},
221:   {"PI33",     3, {  2, -1,  0 }, {   0,  0 }},
222:   {"PI34",    10, {  7, -4,  0 }, {   0,  0 }},

224:   {"PC11",     1, {  2, -1,  0 }, {  -1,  0 }},
225:   {"PC47",    10, { 11, -7,  0 }, { -10,  0 }},
226:   {"PC36",    10, {  9, -6,  0 }, { -10,  0 }},

228:   {"H0211",    2, {  1,  1,  0 }, {   1,  0 }},
229:   {"H211b",    4, {  1,  1,  0 }, {   1,  0 }},
230:   {"H211PI",   6, {  1,  1,  0 }, {   0,  0 }},

232:   {"H0312",    4, {  1,  2,  1 }, {   3,  1 }},
233:   {"H312b",    8, {  1,  2,  1 }, {   3,  1 }},
234:   {"H312PID", 18, {  1,  2,  1 }, {   0,  0 }},

236:   {"H0321",    4, {  5,  2,- 3 }, {  -1, -3 }},
237:   {"H321",    18, {  6,  1,- 5 }, { -15, -3 }},
238: };

240: static PetscErrorCode TSAdaptDSPSetFilter_DSP(TSAdapt adapt,const char *name)
241: {
242:   TSAdapt_DSP       *dsp = (TSAdapt_DSP*)adapt->data;
243:   PetscInt          i,count = (PetscInt)(sizeof(filterlist)/sizeof(filterlist[0]));
244:   struct FilterTab* tab = NULL;
245:   PetscBool         match;
246:   PetscErrorCode    ierr;

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

262: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)
263: {
264:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;

267:   dsp->kBeta[0] = kkI + kkP + kkD;
268:   dsp->kBeta[1] = -(kkP + 2*kkD);
269:   dsp->kBeta[2] = kkD;
270:   dsp->Alpha[0] = 0;
271:   dsp->Alpha[1] = 0;
272:   return(0);
273: }

275: static PetscErrorCode TSAdaptSetFromOptions_DSP(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
276: {
277:   TSAdapt_DSP    *dsp = (TSAdapt_DSP*)adapt->data;
278:   const char     *names[sizeof(filterlist)/sizeof(filterlist[0])];
279:   PetscInt       count = (PetscInt)(sizeof(filterlist)/sizeof(filterlist[0]));
280:   PetscInt       index = 2; /* PI42 */
281:   PetscReal      pid[3] = {1,0,0};
282:   PetscInt       i,n;
283:   PetscBool      set;

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

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

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

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

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

305:   PetscOptionsTail();
306:   return(0);
307: }

309: /*@C
310:    TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter

312:    Collective on TSAdapt

314:    Input Arguments:
315: +  adapt - adaptive controller context
316: -  name - filter name

318:    Level: intermediate

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

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

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

334: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID()
335: @*/
336: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt,const char *name)
337: {
342:   PetscTryMethod(adapt,"TSAdaptDSPSetFilter_C",(TSAdapt,const char*),(adapt,name));
343:   return(0);
344: }

346: /*@
347:    TSAdaptDSPSetPID - Set the PID controller parameters

349:    Input Arguments:
350: +  adapt - adaptive controller context
351: .  kkI - Integral parameter
352: .  kkP - Proportional parameter
353: -  kkD - Derivative parameter

355:    Level: intermediate

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

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

362: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetFilter()
363: @*/
364: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)
365: {
372:   PetscTryMethod(adapt,"TSAdaptDSPSetPID_C",(TSAdapt,PetscReal,PetscReal,PetscReal),(adapt,kkI,kkP,kkD));
373:   return(0);
374: }

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

379:    Level: intermediate

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

383:    Options Database:
384: +   -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
385: .   -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
386: .   -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters
387: -   -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters

389: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID(), TSAdaptDSPSetFilter()
390: M*/
391: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
392: {
393:   TSAdapt_DSP    *dsp;

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

400:   adapt->data                = (void*)dsp;
401:   adapt->ops->choose         = TSAdaptChoose_DSP;
402:   adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
403:   adapt->ops->destroy        = TSAdaptDestroy_DSP;
404:   adapt->ops->reset          = TSAdaptReset_DSP;
405:   adapt->ops->view           = TSAdaptView_DSP;

407:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",TSAdaptDSPSetFilter_DSP);
408:   PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",TSAdaptDSPSetPID_DSP);

410:   TSAdaptDSPSetFilter_DSP(adapt,"PI42");
411:   TSAdaptRestart_DSP(adapt);
412:   return(0);
413: }