Actual source code: adaptdsp.c
petsc-3.10.5 2019-03-28
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 controlers.
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>
334: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetFilter()
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>
362: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID()
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>
385: . -ts_adapt_dsp_pid <kkI,kkP,kkD>
386: . -ts_adapt_dsp_kbeta <b1,b2,b2>
387: - -ts_adapt_dsp_alpha <a2,a3>
389: .seealso: TS, TSAdapt, TSGetAdapt()
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: }