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: dsp->cerror[0] = dsp->hratio[0] = 1.0;
47: dsp->cerror[1] = dsp->hratio[1] = 1.0;
48: dsp->cerror[2] = dsp->hratio[2] = 1.0;
49: return 0;
50: }
52: static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt)
53: {
54: TSAdapt_DSP *dsp = (TSAdapt_DSP*)adapt->data;
55: dsp->cerror[0] = dsp->cerror[1];
56: dsp->cerror[1] = dsp->cerror[2];
57: dsp->cerror[2] = 1.0;
58: dsp->hratio[0] = dsp->hratio[1];
59: dsp->hratio[1] = dsp->hratio[2];
60: dsp->hratio[2] = 1.0;
61: return 0;
62: }
64: static PetscErrorCode TSAdaptChoose_DSP(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter)
65: {
66: TSAdapt_DSP *dsp = (TSAdapt_DSP*)adapt->data;
67: PetscInt order = PETSC_DECIDE;
68: PetscReal enorm = -1;
69: PetscReal enorma,enormr;
70: PetscReal safety = adapt->safety * (PetscReal)0.9;
71: PetscReal hnew,hfac = PETSC_INFINITY;
72: PetscReal hmin = adapt->dt_min*(1 + PETSC_SQRT_MACHINE_EPSILON);
74: *next_sc = 0; /* Reuse the same order scheme */
75: *wltea = -1; /* Weighted absolute local truncation error is not used */
76: *wlter = -1; /* Weighted relative local truncation error is not used */
78: if (ts->ops->evaluatewlte) {
79: TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);
81: } else if (ts->ops->evaluatestep) {
82: DM dm;
83: Vec Y;
87: order = adapt->candidates.order[0];
88: TSGetDM(ts,&dm);
89: DMGetGlobalVector(dm,&Y);
90: TSEvaluateStep(ts,order-1,Y,NULL);
91: TSErrorWeightedNorm(ts,ts->vec_sol,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
92: DMRestoreGlobalVector(dm,&Y);
93: }
94: if (enorm < 0) {
95: TSAdaptRestart_DSP(adapt);
96: *accept = PETSC_TRUE; /* Accept the step */
97: *next_h = h; /* Reuse the old step size */
98: *wlte = -1; /* Weighted local truncation error was not evaluated */
99: return 0;
100: }
102: PetscCitationsRegister(citation[0],&cited[0]);
103: PetscCitationsRegister(citation[1],&cited[1]);
105: /* Update history after rollback */
106: if (!ts->steprollback)
107: dsp->rollback = PETSC_FALSE;
108: else if (!dsp->rollback) {
109: dsp->rollback = PETSC_TRUE;
110: TSAdaptRollBack_DSP(adapt);
111: }
112: /* Reset history after restart */
113: if (ts->steprestart) {
114: TSAdaptRestart_DSP(adapt);
115: }
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)
144: *accept = PETSC_TRUE;
145: else if (adapt->always_accept)
146: *accept = PETSC_TRUE;
147: else if (h < hmin)
148: *accept = PETSC_TRUE;
149: else
150: *accept = PETSC_FALSE;
152: /* Update history after accept */
153: if (*accept) {
154: dsp->cerror[2] = dsp->cerror[1];
155: dsp->cerror[1] = dsp->cerror[0];
156: dsp->cerror[0] = ctr0;
157: dsp->hratio[2] = dsp->hratio[1];
158: dsp->hratio[1] = dsp->hratio[0];
159: dsp->hratio[0] = rho0;
160: dsp->rollback = PETSC_FALSE;
161: }
163: hfac = rho0;
164: }
166: hnew = h * PetscClipInterval(hfac,adapt->clip[0],adapt->clip[1]);
167: *next_h = PetscClipInterval(hnew,adapt->dt_min,adapt->dt_max);
168: *wlte = enorm;
169: return 0;
170: }
172: static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt)
173: {
174: PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",NULL);
175: PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",NULL);
176: PetscFree(adapt->data);
177: return 0;
178: }
180: static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt,PetscViewer viewer)
181: {
182: TSAdapt_DSP *dsp = (TSAdapt_DSP*)adapt->data;
183: PetscBool iascii;
185: 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: PetscViewerASCIIPrintf(viewer,"filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n",b1,b2,b3,a2,a3);
190: }
191: return 0;
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 = (PetscInt)(sizeof(filterlist)/sizeof(filterlist[0]));
229: struct FilterTab* tab = NULL;
230: PetscBool match;
232: for (i=0; i<count; i++) {
233: PetscStrcasecmp(name,filterlist[i].name,&match);
234: if (match) { tab = &filterlist[i]; break; }
235: }
237: dsp->kBeta[0] = tab->kBeta[0]/tab->scale;
238: dsp->kBeta[1] = tab->kBeta[1]/tab->scale;
239: dsp->kBeta[2] = tab->kBeta[2]/tab->scale;
240: dsp->Alpha[0] = tab->Alpha[0]/tab->scale;
241: dsp->Alpha[1] = tab->Alpha[1]/tab->scale;
242: return 0;
243: }
245: static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)
246: {
247: TSAdapt_DSP *dsp = (TSAdapt_DSP*)adapt->data;
249: dsp->kBeta[0] = kkI + kkP + kkD;
250: dsp->kBeta[1] = -(kkP + 2*kkD);
251: dsp->kBeta[2] = kkD;
252: dsp->Alpha[0] = 0;
253: dsp->Alpha[1] = 0;
254: return 0;
255: }
257: static PetscErrorCode TSAdaptSetFromOptions_DSP(PetscOptionItems *PetscOptionsObject,TSAdapt adapt)
258: {
259: TSAdapt_DSP *dsp = (TSAdapt_DSP*)adapt->data;
260: const char *names[sizeof(filterlist)/sizeof(filterlist[0])];
261: PetscInt count = (PetscInt)(sizeof(filterlist)/sizeof(filterlist[0]));
262: PetscInt index = 2; /* PI42 */
263: PetscReal pid[3] = {1,0,0};
264: PetscInt i,n;
265: PetscBool set;
267: for (i=0; i<count; i++) names[i] = filterlist[i].name;
268: PetscOptionsHead(PetscOptionsObject,"DSP adaptive controller options");
270: PetscOptionsEList("-ts_adapt_dsp_filter","Filter name","TSAdaptDSPSetFilter",names,count,names[index],&index,&set);
271: if (set) TSAdaptDSPSetFilter(adapt,names[index]);
273: PetscOptionsRealArray("-ts_adapt_dsp_pid","PID parameters <kkI,kkP,kkD>","TSAdaptDSPSetPID",pid,(n=3,&n),&set);
275: if (set) TSAdaptDSPSetPID(adapt,pid[0],pid[1],pid[2]);
277: PetscOptionsRealArray("-ts_adapt_dsp_kbeta","Filter parameters","",dsp->kBeta,(n=3,&n),&set);
279: if (set) for (i=n; i<3; i++) dsp->kBeta[i] = 0;
281: PetscOptionsRealArray("-ts_adapt_dsp_alpha","Filter parameters","",dsp->Alpha,(n=2,&n),&set);
283: if (set) for (i=n; i<2; i++) dsp->Alpha[i] = 0;
285: PetscOptionsTail();
286: return 0;
287: }
289: /*@C
290: TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter
292: Collective on TSAdapt
294: Input Parameters:
295: + adapt - adaptive controller context
296: - name - filter name
298: Level: intermediate
300: References:
301: . * - http://dx.doi.org/10.1145/641876.641877
303: Notes:
304: Valid filter names are
305: + "basic" - similar to TSADAPTBASIC but with different criteria for step rejections.
306: . "PI30", "PI42", "PI33", "PI34" - PI controllers.
307: . "PC11", "PC47", "PC36" - predictive controllers.
308: . "H0211", "H211b", "H211PI" - digital filters with orders dynamics=2, adaptivity=1, filter=1.
309: . "H0312", "H312b", "H312PID" - digital filters with orders dynamics=3, adaptivity=1, filter=2.
310: - "H0321", "H321" - digital filters with orders dynamics=3, adaptivity=2, filter=1.
312: Options Database:
313: . -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
315: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID()
316: @*/
317: PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt,const char *name)
318: {
321: PetscTryMethod(adapt,"TSAdaptDSPSetFilter_C",(TSAdapt,const char*),(adapt,name));
322: return 0;
323: }
325: /*@
326: TSAdaptDSPSetPID - Set the PID controller parameters
328: Input Parameters:
329: + adapt - adaptive controller context
330: . kkI - Integral parameter
331: . kkP - Proportional parameter
332: - kkD - Derivative parameter
334: Level: intermediate
336: References:
337: . * - http://dx.doi.org/10.1016/j.cam.2005.03.008
339: Options Database:
340: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
342: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetFilter()
343: @*/
344: PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD)
345: {
350: PetscTryMethod(adapt,"TSAdaptDSPSetPID_C",(TSAdapt,PetscReal,PetscReal,PetscReal),(adapt,kkI,kkP,kkD));
351: return 0;
352: }
354: /*MC
355: TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP)
357: Level: intermediate
359: References:
360: + * - http://dx.doi.org/10.1145/641876.641877
361: - * - http://dx.doi.org/10.1016/j.cam.2005.03.008
363: Options Database:
364: + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers
365: . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters
366: . -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters
367: - -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters
369: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID(), TSAdaptDSPSetFilter()
370: M*/
371: PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
372: {
373: TSAdapt_DSP *dsp;
375: PetscNewLog(adapt,&dsp);
376: adapt->reject_safety = 1.0; /* unused */
378: adapt->data = (void*)dsp;
379: adapt->ops->choose = TSAdaptChoose_DSP;
380: adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
381: adapt->ops->destroy = TSAdaptDestroy_DSP;
382: adapt->ops->view = TSAdaptView_DSP;
384: PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",TSAdaptDSPSetFilter_DSP);
385: PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",TSAdaptDSPSetPID_DSP);
387: TSAdaptDSPSetFilter_DSP(adapt,"PI42");
388: TSAdaptRestart_DSP(adapt);
389: return 0;
390: }