Actual source code: bdf.c
1: /*
2: Code for timestepping with BDF methods
3: */
4: #include <petsc/private/tsimpl.h>
5: #include <petscdm.h>
7: static PetscBool cited = PETSC_FALSE;
8: static const char citation[] =
9: "@book{Brenan1995,\n"
10: " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n"
11: " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n"
12: " publisher = {Society for Industrial and Applied Mathematics},\n"
13: " year = {1995},\n"
14: " doi = {10.1137/1.9781611971224},\n}\n";
16: typedef struct {
17: PetscInt k,n;
18: PetscReal time[6+2];
19: Vec work[6+2];
20: Vec tvwork[6+2];
21: PetscReal shift;
22: Vec vec_dot; /* Xdot when !transientvar, else Cdot where C(X) is the transient variable. */
23: Vec vec_wrk;
24: Vec vec_lte;
26: PetscBool transientvar;
27: PetscInt order;
28: TSStepStatus status;
29: } TS_BDF;
31: /* Compute Lagrange polynomials on T[:n] evaluated at t.
32: * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i].
33: */
34: static inline void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
35: {
36: PetscInt k,j;
37: for (k=0; k<n; k++)
38: for (L[k]=1, j=0; j<n; j++)
39: if (j != k)
40: L[k] *= (t - T[j])/(T[k] - T[j]);
41: }
43: static inline void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
44: {
45: PetscInt k,j,i;
46: for (k=0; k<n; k++)
47: for (dL[k]=0, j=0; j<n; j++)
48: if (j != k) {
49: PetscReal L = 1/(T[k] - T[j]);
50: for (i=0; i<n; i++)
51: if (i != j && i != k)
52: L *= (t - T[i])/(T[k] - T[i]);
53: dL[k] += L;
54: }
55: }
57: static PetscErrorCode TSBDF_GetVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
58: {
59: TS_BDF *bdf = (TS_BDF*)ts->data;
61: if (dm && dm != ts->dm) {
62: DMGetNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
63: DMGetNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
64: } else {
65: *Xdot = bdf->vec_dot;
66: *Ydot = bdf->vec_wrk;
67: }
68: return 0;
69: }
71: static PetscErrorCode TSBDF_RestoreVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
72: {
73: TS_BDF *bdf = (TS_BDF*)ts->data;
75: if (dm && dm != ts->dm) {
76: DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
77: DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
78: } else {
81: *Xdot = NULL;
82: *Ydot = NULL;
83: }
84: return 0;
85: }
87: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine,DM coarse,void *ctx)
88: {
89: return 0;
90: }
92: static PetscErrorCode DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
93: {
94: TS ts = (TS)ctx;
95: Vec Ydot,Ydot_c;
96: Vec Xdot,Xdot_c;
98: TSBDF_GetVecs(ts,fine,&Xdot,&Ydot);
99: TSBDF_GetVecs(ts,coarse,&Xdot_c,&Ydot_c);
101: MatRestrict(restrct,Ydot,Ydot_c);
102: VecPointwiseMult(Ydot_c,rscale,Ydot_c);
104: TSBDF_RestoreVecs(ts,fine,&Xdot,&Ydot);
105: TSBDF_RestoreVecs(ts,coarse,&Xdot_c,&Ydot_c);
106: return 0;
107: }
109: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
110: {
111: TS_BDF *bdf = (TS_BDF*)ts->data;
112: PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
113: Vec tail = bdf->work[n-1],tvtail = bdf->tvwork[n-1];
115: for (i=n-1; i>=2; i--) {
116: bdf->time[i] = bdf->time[i-1];
117: bdf->work[i] = bdf->work[i-1];
118: bdf->tvwork[i] = bdf->tvwork[i-1];
119: }
120: bdf->n = PetscMin(bdf->n+1,n-1);
121: bdf->time[1] = t;
122: bdf->work[1] = tail;
123: bdf->tvwork[1] = tvtail;
124: VecCopy(X,tail);
125: TSComputeTransientVariable(ts,tail,tvtail);
126: return 0;
127: }
129: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
130: {
131: TS_BDF *bdf = (TS_BDF*)ts->data;
132: PetscInt i,n = order+1;
133: PetscReal *time = bdf->time;
134: Vec *vecs = bdf->work;
135: PetscScalar a[8],b[8],alpha[8];
137: LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
138: LagrangeBasisDers(n+1,time[0],time,b);
139: for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
140: VecZeroEntries(lte);
141: VecMAXPY(lte,n+1,alpha,vecs);
142: return 0;
143: }
145: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
146: {
147: TS_BDF *bdf = (TS_BDF*)ts->data;
148: PetscInt n = order+1;
149: PetscReal *time = bdf->time+1;
150: Vec *vecs = bdf->work+1;
151: PetscScalar alpha[7];
153: n = PetscMin(n,bdf->n);
154: LagrangeBasisVals(n,t,time,alpha);
155: VecZeroEntries(X);
156: VecMAXPY(X,n,alpha,vecs);
157: return 0;
158: }
160: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
161: {
162: TS_BDF *bdf = (TS_BDF*)ts->data;
163: PetscInt n = order+1;
164: PetscReal *time = bdf->time;
165: Vec *vecs = bdf->work;
166: PetscScalar alpha[7];
168: LagrangeBasisVals(n,t,time,alpha);
169: VecZeroEntries(X);
170: VecMAXPY(X,n,alpha,vecs);
171: return 0;
172: }
174: /* Compute the affine term V0 such that Xdot = shift*X + V0.
175: *
176: * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
177: */
178: static PetscErrorCode TSBDF_PreSolve(TS ts)
179: {
180: TS_BDF *bdf = (TS_BDF*)ts->data;
181: PetscInt i,n = PetscMax(bdf->k,1) + 1;
182: Vec V,V0;
183: Vec vecs[7];
184: PetscScalar alpha[7];
186: TSBDF_GetVecs(ts,NULL,&V,&V0);
187: LagrangeBasisDers(n,bdf->time[0],bdf->time,alpha);
188: for (i=1; i<n; i++) {
189: vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
190: }
191: VecZeroEntries(V0);
192: VecMAXPY(V0,n-1,alpha+1,vecs+1);
193: bdf->shift = PetscRealPart(alpha[0]);
194: TSBDF_RestoreVecs(ts,NULL,&V,&V0);
195: return 0;
196: }
198: static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x)
199: {
200: PetscInt nits,lits;
202: TSBDF_PreSolve(ts);
203: SNESSolve(ts->snes,b,x);
204: SNESGetIterationNumber(ts->snes,&nits);
205: SNESGetLinearSolveIterations(ts->snes,&lits);
206: ts->snes_its += nits; ts->ksp_its += lits;
207: return 0;
208: }
210: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
211: {
212: TS_BDF *bdf = (TS_BDF*)ts->data;
214: bdf->k = 1; bdf->n = 0;
215: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
217: bdf->time[0] = ts->ptime + ts->time_step/2;
218: VecCopy(bdf->work[1],bdf->work[0]);
219: TSPreStage(ts,bdf->time[0]);
220: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
221: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
222: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
223: if (!*accept) return 0;
225: bdf->k = PetscMin(2,bdf->order); bdf->n++;
226: VecCopy(bdf->work[0],bdf->work[2]);
227: bdf->time[2] = bdf->time[0];
228: TSComputeTransientVariable(ts,bdf->work[2],bdf->tvwork[2]);
229: return 0;
230: }
232: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
234: static PetscErrorCode TSStep_BDF(TS ts)
235: {
236: TS_BDF *bdf = (TS_BDF*)ts->data;
237: PetscInt rejections = 0;
238: PetscBool stageok,accept = PETSC_TRUE;
239: PetscReal next_time_step = ts->time_step;
241: PetscCitationsRegister(citation,&cited);
243: if (!ts->steprollback && !ts->steprestart) {
244: bdf->k = PetscMin(bdf->k+1,bdf->order);
245: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
246: }
248: bdf->status = TS_STEP_INCOMPLETE;
249: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
251: if (ts->steprestart) {
252: TSBDF_Restart(ts,&stageok);
253: if (!stageok) goto reject_step;
254: }
256: bdf->time[0] = ts->ptime + ts->time_step;
257: TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
258: TSPreStage(ts,bdf->time[0]);
259: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
260: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
261: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
262: if (!stageok) goto reject_step;
264: bdf->status = TS_STEP_PENDING;
265: TSAdaptCandidatesClear(ts->adapt);
266: TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
267: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
268: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
269: if (!accept) { ts->time_step = next_time_step; goto reject_step; }
271: VecCopy(bdf->work[0],ts->vec_sol);
272: ts->ptime += ts->time_step;
273: ts->time_step = next_time_step;
274: break;
276: reject_step:
277: ts->reject++; accept = PETSC_FALSE;
278: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
279: PetscInfo(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
280: ts->reason = TS_DIVERGED_STEP_REJECTED;
281: }
282: }
283: return 0;
284: }
286: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
287: {
288: TS_BDF *bdf = (TS_BDF*)ts->data;
290: TSBDF_Interpolate(ts,bdf->k,t,X);
291: return 0;
292: }
294: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
295: {
296: TS_BDF *bdf = (TS_BDF*)ts->data;
297: PetscInt k = bdf->k;
298: PetscReal wltea,wlter;
299: Vec X = bdf->work[0], Y = bdf->vec_lte;
301: k = PetscMin(k,bdf->n-1);
302: TSBDF_VecLTE(ts,k,Y);
303: VecAXPY(Y,1,X);
304: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
305: if (order) *order = k + 1;
306: return 0;
307: }
309: static PetscErrorCode TSRollBack_BDF(TS ts)
310: {
311: TS_BDF *bdf = (TS_BDF*)ts->data;
313: VecCopy(bdf->work[1],ts->vec_sol);
314: return 0;
315: }
317: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts)
318: {
319: TS_BDF *bdf = (TS_BDF*)ts->data;
320: DM dm, dmsave = ts->dm;
321: PetscReal t = bdf->time[0];
322: PetscReal shift = bdf->shift;
323: Vec V,V0;
325: SNESGetDM(snes,&dm);
326: TSBDF_GetVecs(ts,dm,&V,&V0);
327: if (bdf->transientvar) { /* shift*C(X) + V0 */
328: TSComputeTransientVariable(ts,X,V);
329: VecAYPX(V,shift,V0);
330: } else { /* shift*X + V0 */
331: VecWAXPY(V,shift,X,V0);
332: }
334: /* F = Function(t,X,V) */
335: ts->dm = dm;
336: TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
337: ts->dm = dmsave;
339: TSBDF_RestoreVecs(ts,dm,&V,&V0);
340: return 0;
341: }
343: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts)
344: {
345: TS_BDF *bdf = (TS_BDF*)ts->data;
346: DM dm, dmsave = ts->dm;
347: PetscReal t = bdf->time[0];
348: PetscReal shift = bdf->shift;
349: Vec V,V0;
351: SNESGetDM(snes,&dm);
352: TSBDF_GetVecs(ts,dm,&V,&V0);
354: /* J,P = Jacobian(t,X,V) */
355: ts->dm = dm;
356: TSComputeIJacobian(ts,t,X,V,shift,J,P,PETSC_FALSE);
357: ts->dm = dmsave;
359: TSBDF_RestoreVecs(ts,dm,&V,&V0);
360: return 0;
361: }
363: static PetscErrorCode TSReset_BDF(TS ts)
364: {
365: TS_BDF *bdf = (TS_BDF*)ts->data;
366: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
368: bdf->k = bdf->n = 0;
369: for (i=0; i<n; i++) {
370: VecDestroy(&bdf->work[i]);
371: VecDestroy(&bdf->tvwork[i]);
372: }
373: VecDestroy(&bdf->vec_dot);
374: VecDestroy(&bdf->vec_wrk);
375: VecDestroy(&bdf->vec_lte);
376: if (ts->dm) DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);
377: return 0;
378: }
380: static PetscErrorCode TSDestroy_BDF(TS ts)
381: {
382: TSReset_BDF(ts);
383: PetscFree(ts->data);
384: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
385: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
386: return 0;
387: }
389: static PetscErrorCode TSSetUp_BDF(TS ts)
390: {
391: TS_BDF *bdf = (TS_BDF*)ts->data;
392: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
393: PetscReal low,high,two = 2;
395: TSHasTransientVariable(ts,&bdf->transientvar);
396: bdf->k = bdf->n = 0;
397: for (i=0; i<n; i++) {
398: VecDuplicate(ts->vec_sol,&bdf->work[i]);
399: if (i && bdf->transientvar) {
400: VecDuplicate(ts->vec_sol,&bdf->tvwork[i]);
401: }
402: }
403: VecDuplicate(ts->vec_sol,&bdf->vec_dot);
404: VecDuplicate(ts->vec_sol,&bdf->vec_wrk);
405: VecDuplicate(ts->vec_sol,&bdf->vec_lte);
406: TSGetDM(ts,&ts->dm);
407: DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);
409: TSGetAdapt(ts,&ts->adapt);
410: TSAdaptCandidatesClear(ts->adapt);
411: TSAdaptGetClip(ts->adapt,&low,&high);
412: TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));
414: TSGetSNES(ts,&ts->snes);
415: return 0;
416: }
418: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
419: {
420: PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
421: {
422: PetscBool flg;
423: PetscInt order;
424: TSBDFGetOrder(ts,&order);
425: PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
426: if (flg) TSBDFSetOrder(ts,order);
427: }
428: PetscOptionsTail();
429: return 0;
430: }
432: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
433: {
434: TS_BDF *bdf = (TS_BDF*)ts->data;
435: PetscBool iascii;
437: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
438: if (iascii) {
439: PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);
440: }
441: return 0;
442: }
444: /* ------------------------------------------------------------ */
446: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
447: {
448: TS_BDF *bdf = (TS_BDF*)ts->data;
450: if (order == bdf->order) return 0;
452: bdf->order = order;
453: return 0;
454: }
456: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
457: {
458: TS_BDF *bdf = (TS_BDF*)ts->data;
460: *order = bdf->order;
461: return 0;
462: }
464: /* ------------------------------------------------------------ */
466: /*MC
467: TSBDF - DAE solver using BDF methods
469: Level: beginner
471: .seealso: TS, TSCreate(), TSSetType()
472: M*/
473: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
474: {
475: TS_BDF *bdf;
477: ts->ops->reset = TSReset_BDF;
478: ts->ops->destroy = TSDestroy_BDF;
479: ts->ops->view = TSView_BDF;
480: ts->ops->setup = TSSetUp_BDF;
481: ts->ops->setfromoptions = TSSetFromOptions_BDF;
482: ts->ops->step = TSStep_BDF;
483: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
484: ts->ops->rollback = TSRollBack_BDF;
485: ts->ops->interpolate = TSInterpolate_BDF;
486: ts->ops->snesfunction = SNESTSFormFunction_BDF;
487: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
488: ts->default_adapt_type = TSADAPTBASIC;
490: ts->usessnes = PETSC_TRUE;
492: PetscNewLog(ts,&bdf);
493: ts->data = (void*)bdf;
495: bdf->status = TS_STEP_COMPLETE;
497: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
498: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
499: TSBDFSetOrder(ts,2);
500: return 0;
501: }
503: /* ------------------------------------------------------------ */
505: /*@
506: TSBDFSetOrder - Set the order of the BDF method
508: Logically Collective on TS
510: Input Parameters:
511: + ts - timestepping context
512: - order - order of the method
514: Options Database:
515: . -ts_bdf_order <order> - select the order
517: Level: intermediate
519: @*/
520: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
521: {
524: PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
525: return 0;
526: }
528: /*@
529: TSBDFGetOrder - Get the order of the BDF method
531: Not Collective
533: Input Parameter:
534: . ts - timestepping context
536: Output Parameter:
537: . order - order of the method
539: Level: intermediate
541: @*/
542: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
543: {
546: PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
547: return 0;
548: }