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[] = "@book{Brenan1995,\n"
9: " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n"
10: " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n"
11: " publisher = {Society for Industrial and Applied Mathematics},\n"
12: " year = {1995},\n"
13: " doi = {10.1137/1.9781611971224},\n}\n";
15: typedef struct {
16: PetscInt k, n;
17: PetscReal time[6 + 2];
18: Vec work[6 + 2];
19: Vec tvwork[6 + 2];
20: PetscReal shift;
21: Vec vec_dot; /* Xdot when !transientvar, else Cdot where C(X) is the transient variable. */
22: Vec vec_wrk;
23: Vec vec_lte;
25: PetscBool transientvar;
26: PetscInt order;
27: TSStepStatus status;
28: } TS_BDF;
30: /* Compute Lagrange polynomials on T[:n] evaluated at t.
31: * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i].
32: */
33: static inline void LagrangeBasisVals(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar L[])
34: {
35: PetscInt k, j;
36: for (k = 0; k < n; k++)
37: for (L[k] = 1, j = 0; j < n; j++)
38: if (j != k) L[k] *= (t - T[j]) / (T[k] - T[j]);
39: }
41: static inline void LagrangeBasisDers(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar dL[])
42: {
43: PetscInt k, j, i;
44: for (k = 0; k < n; k++)
45: for (dL[k] = 0, j = 0; j < n; j++)
46: if (j != k) {
47: PetscReal L = 1 / (T[k] - T[j]);
48: for (i = 0; i < n; i++)
49: if (i != j && i != k) L *= (t - T[i]) / (T[k] - T[i]);
50: dL[k] += L;
51: }
52: }
54: static PetscErrorCode TSBDF_GetVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
55: {
56: TS_BDF *bdf = (TS_BDF *)ts->data;
58: if (dm && dm != ts->dm) {
59: DMGetNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot);
60: DMGetNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot);
61: } else {
62: *Xdot = bdf->vec_dot;
63: *Ydot = bdf->vec_wrk;
64: }
65: return 0;
66: }
68: static PetscErrorCode TSBDF_RestoreVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
69: {
70: TS_BDF *bdf = (TS_BDF *)ts->data;
72: if (dm && dm != ts->dm) {
73: DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot);
74: DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot);
75: } else {
78: *Xdot = NULL;
79: *Ydot = NULL;
80: }
81: return 0;
82: }
84: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine, DM coarse, void *ctx)
85: {
86: return 0;
87: }
89: static PetscErrorCode DMRestrictHook_TSBDF(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
90: {
91: TS ts = (TS)ctx;
92: Vec Ydot, Ydot_c;
93: Vec Xdot, Xdot_c;
95: TSBDF_GetVecs(ts, fine, &Xdot, &Ydot);
96: TSBDF_GetVecs(ts, coarse, &Xdot_c, &Ydot_c);
98: MatRestrict(restrct, Ydot, Ydot_c);
99: VecPointwiseMult(Ydot_c, rscale, Ydot_c);
101: TSBDF_RestoreVecs(ts, fine, &Xdot, &Ydot);
102: TSBDF_RestoreVecs(ts, coarse, &Xdot_c, &Ydot_c);
103: return 0;
104: }
106: static PetscErrorCode TSBDF_Advance(TS ts, PetscReal t, Vec X)
107: {
108: TS_BDF *bdf = (TS_BDF *)ts->data;
109: PetscInt i, n = (PetscInt)(sizeof(bdf->work) / sizeof(Vec));
110: Vec tail = bdf->work[n - 1], tvtail = bdf->tvwork[n - 1];
112: for (i = n - 1; i >= 2; i--) {
113: bdf->time[i] = bdf->time[i - 1];
114: bdf->work[i] = bdf->work[i - 1];
115: bdf->tvwork[i] = bdf->tvwork[i - 1];
116: }
117: bdf->n = PetscMin(bdf->n + 1, n - 1);
118: bdf->time[1] = t;
119: bdf->work[1] = tail;
120: bdf->tvwork[1] = tvtail;
121: VecCopy(X, tail);
122: TSComputeTransientVariable(ts, tail, tvtail);
123: return 0;
124: }
126: static PetscErrorCode TSBDF_VecLTE(TS ts, PetscInt order, Vec lte)
127: {
128: TS_BDF *bdf = (TS_BDF *)ts->data;
129: PetscInt i, n = order + 1;
130: PetscReal *time = bdf->time;
131: Vec *vecs = bdf->work;
132: PetscScalar a[8], b[8], alpha[8];
134: LagrangeBasisDers(n + 0, time[0], time, a);
135: a[n] = 0;
136: LagrangeBasisDers(n + 1, time[0], time, b);
137: for (i = 0; i < n + 1; i++) alpha[i] = (a[i] - b[i]) / a[0];
138: VecZeroEntries(lte);
139: VecMAXPY(lte, n + 1, alpha, vecs);
140: return 0;
141: }
143: static PetscErrorCode TSBDF_Extrapolate(TS ts, PetscInt order, PetscReal t, Vec X)
144: {
145: TS_BDF *bdf = (TS_BDF *)ts->data;
146: PetscInt n = order + 1;
147: PetscReal *time = bdf->time + 1;
148: Vec *vecs = bdf->work + 1;
149: PetscScalar alpha[7];
151: n = PetscMin(n, bdf->n);
152: LagrangeBasisVals(n, t, time, alpha);
153: VecZeroEntries(X);
154: VecMAXPY(X, n, alpha, vecs);
155: return 0;
156: }
158: static PetscErrorCode TSBDF_Interpolate(TS ts, PetscInt order, PetscReal t, Vec X)
159: {
160: TS_BDF *bdf = (TS_BDF *)ts->data;
161: PetscInt n = order + 1;
162: PetscReal *time = bdf->time;
163: Vec *vecs = bdf->work;
164: PetscScalar alpha[7];
166: LagrangeBasisVals(n, t, time, alpha);
167: VecZeroEntries(X);
168: VecMAXPY(X, n, alpha, vecs);
169: return 0;
170: }
172: /* Compute the affine term V0 such that Xdot = shift*X + V0.
173: *
174: * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
175: */
176: static PetscErrorCode TSBDF_PreSolve(TS ts)
177: {
178: TS_BDF *bdf = (TS_BDF *)ts->data;
179: PetscInt i, n = PetscMax(bdf->k, 1) + 1;
180: Vec V, V0;
181: Vec vecs[7];
182: PetscScalar alpha[7];
184: TSBDF_GetVecs(ts, NULL, &V, &V0);
185: LagrangeBasisDers(n, bdf->time[0], bdf->time, alpha);
186: for (i = 1; i < n; i++) vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
187: VecZeroEntries(V0);
188: VecMAXPY(V0, n - 1, alpha + 1, vecs + 1);
189: bdf->shift = PetscRealPart(alpha[0]);
190: TSBDF_RestoreVecs(ts, NULL, &V, &V0);
191: return 0;
192: }
194: static PetscErrorCode TSBDF_SNESSolve(TS ts, Vec b, Vec x)
195: {
196: PetscInt nits, lits;
198: TSBDF_PreSolve(ts);
199: SNESSolve(ts->snes, b, x);
200: SNESGetIterationNumber(ts->snes, &nits);
201: SNESGetLinearSolveIterations(ts->snes, &lits);
202: ts->snes_its += nits;
203: ts->ksp_its += lits;
204: return 0;
205: }
207: static PetscErrorCode TSBDF_Restart(TS ts, PetscBool *accept)
208: {
209: TS_BDF *bdf = (TS_BDF *)ts->data;
211: bdf->k = 1;
212: bdf->n = 0;
213: TSBDF_Advance(ts, ts->ptime, ts->vec_sol);
215: bdf->time[0] = ts->ptime + ts->time_step / 2;
216: VecCopy(bdf->work[1], bdf->work[0]);
217: TSPreStage(ts, bdf->time[0]);
218: TSBDF_SNESSolve(ts, NULL, bdf->work[0]);
219: TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]);
220: TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], accept);
221: if (!*accept) return 0;
223: bdf->k = PetscMin(2, bdf->order);
224: bdf->n++;
225: VecCopy(bdf->work[0], bdf->work[2]);
226: bdf->time[2] = bdf->time[0];
227: TSComputeTransientVariable(ts, bdf->work[2], bdf->tvwork[2]);
228: return 0;
229: }
231: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
233: static PetscErrorCode TSStep_BDF(TS ts)
234: {
235: TS_BDF *bdf = (TS_BDF *)ts->data;
236: PetscInt rejections = 0;
237: PetscBool stageok, accept = PETSC_TRUE;
238: PetscReal next_time_step = ts->time_step;
240: PetscCitationsRegister(citation, &cited);
242: if (!ts->steprollback && !ts->steprestart) {
243: bdf->k = PetscMin(bdf->k + 1, bdf->order);
244: TSBDF_Advance(ts, ts->ptime, ts->vec_sol);
245: }
247: bdf->status = TS_STEP_INCOMPLETE;
248: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
249: if (ts->steprestart) {
250: TSBDF_Restart(ts, &stageok);
251: if (!stageok) goto reject_step;
252: }
254: bdf->time[0] = ts->ptime + ts->time_step;
255: TSBDF_Extrapolate(ts, bdf->k - (accept ? 0 : 1), bdf->time[0], bdf->work[0]);
256: TSPreStage(ts, bdf->time[0]);
257: TSBDF_SNESSolve(ts, NULL, bdf->work[0]);
258: TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]);
259: TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], &stageok);
260: if (!stageok) goto reject_step;
262: bdf->status = TS_STEP_PENDING;
263: TSAdaptCandidatesClear(ts->adapt);
264: TSAdaptCandidateAdd(ts->adapt, BDF_SchemeName[bdf->k], bdf->k, 1, 1.0, 1.0, PETSC_TRUE);
265: TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
266: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
267: if (!accept) {
268: ts->time_step = next_time_step;
269: goto reject_step;
270: }
272: VecCopy(bdf->work[0], ts->vec_sol);
273: ts->ptime += ts->time_step;
274: ts->time_step = next_time_step;
275: break;
277: reject_step:
278: ts->reject++;
279: accept = PETSC_FALSE;
280: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
281: PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections);
282: ts->reason = TS_DIVERGED_STEP_REJECTED;
283: }
284: }
285: return 0;
286: }
288: static PetscErrorCode TSInterpolate_BDF(TS ts, PetscReal t, Vec X)
289: {
290: TS_BDF *bdf = (TS_BDF *)ts->data;
292: TSBDF_Interpolate(ts, bdf->k, t, X);
293: return 0;
294: }
296: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
297: {
298: TS_BDF *bdf = (TS_BDF *)ts->data;
299: PetscInt k = bdf->k;
300: PetscReal wltea, wlter;
301: Vec X = bdf->work[0], Y = bdf->vec_lte;
303: k = PetscMin(k, bdf->n - 1);
304: TSBDF_VecLTE(ts, k, Y);
305: VecAXPY(Y, 1, X);
306: TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter);
307: if (order) *order = k + 1;
308: return 0;
309: }
311: static PetscErrorCode TSRollBack_BDF(TS ts)
312: {
313: TS_BDF *bdf = (TS_BDF *)ts->data;
315: VecCopy(bdf->work[1], ts->vec_sol);
316: return 0;
317: }
319: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes, Vec X, Vec F, TS ts)
320: {
321: TS_BDF *bdf = (TS_BDF *)ts->data;
322: DM dm, dmsave = ts->dm;
323: PetscReal t = bdf->time[0];
324: PetscReal shift = bdf->shift;
325: Vec V, V0;
327: SNESGetDM(snes, &dm);
328: TSBDF_GetVecs(ts, dm, &V, &V0);
329: if (bdf->transientvar) { /* shift*C(X) + V0 */
330: TSComputeTransientVariable(ts, X, V);
331: VecAYPX(V, shift, V0);
332: } else { /* shift*X + V0 */
333: VecWAXPY(V, shift, X, V0);
334: }
336: /* F = Function(t,X,V) */
337: ts->dm = dm;
338: TSComputeIFunction(ts, t, X, V, F, PETSC_FALSE);
339: ts->dm = dmsave;
341: TSBDF_RestoreVecs(ts, dm, &V, &V0);
342: return 0;
343: }
345: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes, Vec X, Mat J, Mat P, TS ts)
346: {
347: TS_BDF *bdf = (TS_BDF *)ts->data;
348: DM dm, dmsave = ts->dm;
349: PetscReal t = bdf->time[0];
350: PetscReal shift = bdf->shift;
351: Vec V, V0;
353: SNESGetDM(snes, &dm);
354: TSBDF_GetVecs(ts, dm, &V, &V0);
356: /* J,P = Jacobian(t,X,V) */
357: ts->dm = dm;
358: TSComputeIJacobian(ts, t, X, V, shift, J, P, PETSC_FALSE);
359: ts->dm = dmsave;
361: TSBDF_RestoreVecs(ts, dm, &V, &V0);
362: return 0;
363: }
365: static PetscErrorCode TSReset_BDF(TS ts)
366: {
367: TS_BDF *bdf = (TS_BDF *)ts->data;
368: size_t i, n = sizeof(bdf->work) / sizeof(Vec);
370: bdf->k = bdf->n = 0;
371: for (i = 0; i < n; i++) {
372: VecDestroy(&bdf->work[i]);
373: VecDestroy(&bdf->tvwork[i]);
374: }
375: VecDestroy(&bdf->vec_dot);
376: VecDestroy(&bdf->vec_wrk);
377: VecDestroy(&bdf->vec_lte);
378: if (ts->dm) DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts);
379: return 0;
380: }
382: static PetscErrorCode TSDestroy_BDF(TS ts)
383: {
384: TSReset_BDF(ts);
385: PetscFree(ts->data);
386: PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", NULL);
387: PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", NULL);
388: return 0;
389: }
391: static PetscErrorCode TSSetUp_BDF(TS ts)
392: {
393: TS_BDF *bdf = (TS_BDF *)ts->data;
394: size_t i, n = sizeof(bdf->work) / sizeof(Vec);
395: PetscReal low, high, two = 2;
397: TSHasTransientVariable(ts, &bdf->transientvar);
398: bdf->k = bdf->n = 0;
399: for (i = 0; i < n; i++) {
400: VecDuplicate(ts->vec_sol, &bdf->work[i]);
401: if (i && bdf->transientvar) VecDuplicate(ts->vec_sol, &bdf->tvwork[i]);
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(TS ts, PetscOptionItems *PetscOptionsObject)
419: {
420: PetscOptionsHeadBegin(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: PetscOptionsHeadEnd();
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) PetscViewerASCIIPrintf(viewer, " Order=%" PetscInt_FMT "\n", bdf->order);
439: return 0;
440: }
442: /* ------------------------------------------------------------ */
444: static PetscErrorCode TSBDFSetOrder_BDF(TS ts, PetscInt order)
445: {
446: TS_BDF *bdf = (TS_BDF *)ts->data;
448: if (order == bdf->order) return 0;
450: bdf->order = order;
451: return 0;
452: }
454: static PetscErrorCode TSBDFGetOrder_BDF(TS ts, PetscInt *order)
455: {
456: TS_BDF *bdf = (TS_BDF *)ts->data;
458: *order = bdf->order;
459: return 0;
460: }
462: /* ------------------------------------------------------------ */
464: /*MC
465: TSBDF - DAE solver using BDF methods
467: Level: beginner
469: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSType`
470: M*/
471: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
472: {
473: TS_BDF *bdf;
475: ts->ops->reset = TSReset_BDF;
476: ts->ops->destroy = TSDestroy_BDF;
477: ts->ops->view = TSView_BDF;
478: ts->ops->setup = TSSetUp_BDF;
479: ts->ops->setfromoptions = TSSetFromOptions_BDF;
480: ts->ops->step = TSStep_BDF;
481: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
482: ts->ops->rollback = TSRollBack_BDF;
483: ts->ops->interpolate = TSInterpolate_BDF;
484: ts->ops->snesfunction = SNESTSFormFunction_BDF;
485: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
486: ts->default_adapt_type = TSADAPTBASIC;
488: ts->usessnes = PETSC_TRUE;
490: PetscNew(&bdf);
491: ts->data = (void *)bdf;
493: bdf->status = TS_STEP_COMPLETE;
495: PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", TSBDFSetOrder_BDF);
496: PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", TSBDFGetOrder_BDF);
497: TSBDFSetOrder(ts, 2);
498: return 0;
499: }
501: /* ------------------------------------------------------------ */
503: /*@
504: TSBDFSetOrder - Set the order of the `TSBDF` method
506: Logically Collective
508: Input Parameters:
509: + ts - timestepping context
510: - order - order of the method
512: Options Database Key:
513: . -ts_bdf_order <order> - select the order
515: Level: intermediate
517: .seealso: `TSBDFGetOrder()`, `TS`, `TSBDF`
518: @*/
519: PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order)
520: {
523: PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order));
524: return 0;
525: }
527: /*@
528: TSBDFGetOrder - Get the order of the `TSBDF` method
530: Not Collective
532: Input Parameter:
533: . ts - timestepping context
535: Output Parameter:
536: . order - order of the method
538: Level: intermediate
540: .seealso: `TSBDFSetOrder()`, `TS`, `TSBDF`
541: @*/
542: PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order)
543: {
546: PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order));
547: return 0;
548: }