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: PetscBool extrapolate;
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) L[k] *= (t - T[j]) / (T[k] - T[j]);
40: }
42: static inline void LagrangeBasisDers(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar dL[])
43: {
44: PetscInt k, j, i;
45: for (k = 0; k < n; k++)
46: for (dL[k] = 0, j = 0; j < n; j++)
47: if (j != k) {
48: PetscReal L = 1 / (T[k] - T[j]);
49: for (i = 0; i < n; i++)
50: if (i != j && i != k) L *= (t - T[i]) / (T[k] - T[i]);
51: dL[k] += L;
52: }
53: }
55: static PetscErrorCode TSBDF_GetVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
56: {
57: TS_BDF *bdf = (TS_BDF *)ts->data;
59: PetscFunctionBegin;
60: if (dm && dm != ts->dm) {
61: PetscCall(DMGetNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot));
62: PetscCall(DMGetNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot));
63: } else {
64: *Xdot = bdf->vec_dot;
65: *Ydot = bdf->vec_wrk;
66: }
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode TSBDF_RestoreVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
71: {
72: TS_BDF *bdf = (TS_BDF *)ts->data;
74: PetscFunctionBegin;
75: if (dm && dm != ts->dm) {
76: PetscCall(DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot));
77: PetscCall(DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot));
78: } else {
79: PetscCheck(*Xdot == bdf->vec_dot, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Vec does not match the cache");
80: PetscCheck(*Ydot == bdf->vec_wrk, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Vec does not match the cache");
81: *Xdot = NULL;
82: *Ydot = NULL;
83: }
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine, DM coarse, void *ctx)
88: {
89: PetscFunctionBegin;
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode DMRestrictHook_TSBDF(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
94: {
95: TS ts = (TS)ctx;
96: Vec Ydot, Ydot_c;
97: Vec Xdot, Xdot_c;
99: PetscFunctionBegin;
100: PetscCall(TSBDF_GetVecs(ts, fine, &Xdot, &Ydot));
101: PetscCall(TSBDF_GetVecs(ts, coarse, &Xdot_c, &Ydot_c));
103: PetscCall(MatRestrict(restrct, Ydot, Ydot_c));
104: PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c));
106: PetscCall(TSBDF_RestoreVecs(ts, fine, &Xdot, &Ydot));
107: PetscCall(TSBDF_RestoreVecs(ts, coarse, &Xdot_c, &Ydot_c));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode TSBDF_Advance(TS ts, PetscReal t, Vec X)
112: {
113: TS_BDF *bdf = (TS_BDF *)ts->data;
114: PetscInt i, n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
115: Vec tail = bdf->work[n - 1], tvtail = bdf->tvwork[n - 1];
117: PetscFunctionBegin;
118: for (i = n - 1; i >= 2; i--) {
119: bdf->time[i] = bdf->time[i - 1];
120: bdf->work[i] = bdf->work[i - 1];
121: bdf->tvwork[i] = bdf->tvwork[i - 1];
122: }
123: bdf->n = PetscMin(bdf->n + 1, n - 1);
124: bdf->time[1] = t;
125: bdf->work[1] = tail;
126: bdf->tvwork[1] = tvtail;
127: PetscCall(VecCopy(X, tail));
128: PetscCall(TSComputeTransientVariable(ts, tail, tvtail));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: static PetscErrorCode TSBDF_VecLTE(TS ts, PetscInt order, Vec lte)
133: {
134: TS_BDF *bdf = (TS_BDF *)ts->data;
135: PetscInt i, n = order + 1;
136: PetscReal *time = bdf->time;
137: Vec *vecs = bdf->work;
138: PetscScalar a[8], b[8], alpha[8];
140: PetscFunctionBegin;
141: LagrangeBasisDers(n + 0, time[0], time, a);
142: a[n] = 0;
143: LagrangeBasisDers(n + 1, time[0], time, b);
144: for (i = 0; i < n + 1; i++) alpha[i] = (a[i] - b[i]) / a[0];
145: PetscCall(VecZeroEntries(lte));
146: PetscCall(VecMAXPY(lte, n + 1, alpha, vecs));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode TSBDF_Extrapolate(TS ts, PetscInt order, PetscReal t, Vec X)
151: {
152: TS_BDF *bdf = (TS_BDF *)ts->data;
153: PetscInt n = order + 1;
154: PetscReal *time = bdf->time + 1;
155: Vec *vecs = bdf->work + 1;
156: PetscScalar alpha[7];
158: PetscFunctionBegin;
159: n = PetscMin(n, bdf->n);
160: LagrangeBasisVals(n, t, time, alpha);
161: PetscCall(VecZeroEntries(X));
162: PetscCall(VecMAXPY(X, n, alpha, vecs));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode TSBDF_Interpolate(TS ts, PetscInt order, PetscReal t, Vec X)
167: {
168: TS_BDF *bdf = (TS_BDF *)ts->data;
169: PetscInt n = order + 1;
170: PetscReal *time = bdf->time;
171: Vec *vecs = bdf->work;
172: PetscScalar alpha[7];
174: PetscFunctionBegin;
175: LagrangeBasisVals(n, t, time, alpha);
176: PetscCall(VecZeroEntries(X));
177: PetscCall(VecMAXPY(X, n, alpha, vecs));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /* Compute the affine term V0 such that Xdot = shift*X + V0.
182: *
183: * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
184: */
185: static PetscErrorCode TSBDF_PreSolve(TS ts)
186: {
187: TS_BDF *bdf = (TS_BDF *)ts->data;
188: PetscInt i, n = PetscMax(bdf->k, 1) + 1;
189: Vec V, V0;
190: Vec vecs[7];
191: PetscScalar alpha[7];
193: PetscFunctionBegin;
194: PetscCall(TSBDF_GetVecs(ts, NULL, &V, &V0));
195: LagrangeBasisDers(n, bdf->time[0], bdf->time, alpha);
196: for (i = 1; i < n; i++) vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
197: PetscCall(VecZeroEntries(V0));
198: PetscCall(VecMAXPY(V0, n - 1, alpha + 1, vecs + 1));
199: bdf->shift = PetscRealPart(alpha[0]);
200: PetscCall(TSBDF_RestoreVecs(ts, NULL, &V, &V0));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode TSBDF_SNESSolve(TS ts, Vec b, Vec x)
205: {
206: PetscInt nits, lits;
208: PetscFunctionBegin;
209: PetscCall(TSBDF_PreSolve(ts));
210: PetscCall(SNESSolve(ts->snes, b, x));
211: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
212: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
213: ts->snes_its += nits;
214: ts->ksp_its += lits;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode TSBDF_Restart(TS ts, PetscBool *accept)
219: {
220: TS_BDF *bdf = (TS_BDF *)ts->data;
222: PetscFunctionBegin;
223: bdf->k = 1;
224: bdf->n = 0;
225: PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol));
226: if (bdf->order == 1) {
227: *accept = PETSC_TRUE;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
230: bdf->time[0] = ts->ptime + ts->time_step / 2;
231: PetscCall(VecCopy(bdf->work[1], bdf->work[0]));
232: PetscCall(TSPreStage(ts, bdf->time[0]));
233: PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
234: PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
235: PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], accept));
236: if (!*accept) PetscFunctionReturn(PETSC_SUCCESS);
238: bdf->k = PetscMin(2, bdf->order);
239: bdf->n++;
240: PetscCall(VecCopy(bdf->work[0], bdf->work[2]));
241: bdf->time[2] = bdf->time[0];
242: PetscCall(TSComputeTransientVariable(ts, bdf->work[2], bdf->tvwork[2]));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
248: static PetscErrorCode TSStep_BDF(TS ts)
249: {
250: TS_BDF *bdf = (TS_BDF *)ts->data;
251: PetscInt rejections = 0;
252: PetscBool stageok, accept = PETSC_TRUE;
253: PetscReal next_time_step = ts->time_step;
255: PetscFunctionBegin;
256: PetscCall(PetscCitationsRegister(citation, &cited));
258: if (!ts->steprollback && !ts->steprestart) {
259: bdf->k = PetscMin(bdf->k + 1, bdf->order);
260: PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol));
261: }
263: bdf->status = TS_STEP_INCOMPLETE;
264: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
265: if (ts->steprestart) {
266: PetscCall(TSBDF_Restart(ts, &stageok));
267: if (!stageok) goto reject_step;
268: }
270: bdf->time[0] = ts->ptime + ts->time_step;
271: if (bdf->extrapolate) PetscCall(TSBDF_Extrapolate(ts, bdf->k - (accept ? 0 : 1), bdf->time[0], bdf->work[0]));
272: PetscCall(TSPreStage(ts, bdf->time[0]));
273: PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
274: PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
275: PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], &stageok));
276: if (!stageok) goto reject_step;
278: bdf->status = TS_STEP_PENDING;
279: PetscCall(TSAdaptCandidatesClear(ts->adapt));
280: PetscCall(TSAdaptCandidateAdd(ts->adapt, BDF_SchemeName[bdf->k], bdf->k, 1, 1.0, 1.0, PETSC_TRUE));
281: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
282: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
283: if (!accept) {
284: ts->time_step = next_time_step;
285: goto reject_step;
286: }
288: PetscCall(VecCopy(bdf->work[0], ts->vec_sol));
289: ts->ptime += ts->time_step;
290: ts->time_step = next_time_step;
291: break;
293: reject_step:
294: ts->reject++;
295: accept = PETSC_FALSE;
296: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
297: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
298: ts->reason = TS_DIVERGED_STEP_REJECTED;
299: }
300: }
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode TSInterpolate_BDF(TS ts, PetscReal t, Vec X)
305: {
306: TS_BDF *bdf = (TS_BDF *)ts->data;
308: PetscFunctionBegin;
309: PetscCall(TSBDF_Interpolate(ts, bdf->k, t, X));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
314: {
315: TS_BDF *bdf = (TS_BDF *)ts->data;
316: PetscInt k = bdf->k;
317: PetscReal wltea, wlter;
318: Vec X = bdf->work[0], Y = bdf->vec_lte;
320: PetscFunctionBegin;
321: k = PetscMin(k, bdf->n - 1);
322: if (k == 0) {
323: *wlte = -1;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
326: PetscCall(TSBDF_VecLTE(ts, k, Y));
327: PetscCall(VecAXPY(Y, 1, X));
328: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
329: if (order) *order = k + 1;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode TSResizeRegister_BDF(TS ts, PetscBool reg)
334: {
335: TS_BDF *bdf = (TS_BDF *)ts->data;
336: const char *names[] = {"", "ts:bdf:1", "ts:bdf:2", "ts:bdf:3", "ts:bdf:4", "ts:bdf:5", "ts:bdf:6", "ts:bdf:7"};
337: PetscInt i, maxn = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
339: PetscFunctionBegin;
340: PetscAssert(maxn == 8, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "names need to be redefined");
341: if (reg) {
342: for (i = 1; i < PetscMin(bdf->n + 1, maxn); i++) PetscCall(TSResizeRegisterVec(ts, names[i], bdf->work[i]));
343: } else {
344: for (i = 1; i < maxn; i++) {
345: PetscCall(TSResizeRetrieveVec(ts, names[i], &bdf->work[i]));
346: if (!bdf->work[i]) break;
347: PetscCall(PetscObjectReference((PetscObject)bdf->work[i]));
348: if (bdf->transientvar) {
349: PetscCall(VecDuplicate(bdf->work[i], &bdf->tvwork[i]));
350: PetscCall(TSComputeTransientVariable(ts, bdf->work[i], bdf->tvwork[i]));
351: }
352: }
353: }
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes, Vec X, Vec F, TS ts)
358: {
359: TS_BDF *bdf = (TS_BDF *)ts->data;
360: DM dm, dmsave = ts->dm;
361: PetscReal t = bdf->time[0];
362: PetscReal shift = bdf->shift;
363: Vec V, V0;
365: PetscFunctionBegin;
366: PetscCall(SNESGetDM(snes, &dm));
367: PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
368: if (bdf->transientvar) { /* shift*C(X) + V0 */
369: PetscCall(TSComputeTransientVariable(ts, X, V));
370: PetscCall(VecAYPX(V, shift, V0));
371: } else { /* shift*X + V0 */
372: PetscCall(VecWAXPY(V, shift, X, V0));
373: }
375: /* F = Function(t,X,V) */
376: ts->dm = dm;
377: PetscCall(TSComputeIFunction(ts, t, X, V, F, PETSC_FALSE));
378: ts->dm = dmsave;
380: PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes, Vec X, Mat J, Mat P, TS ts)
385: {
386: TS_BDF *bdf = (TS_BDF *)ts->data;
387: DM dm, dmsave = ts->dm;
388: PetscReal t = bdf->time[0];
389: PetscReal shift = bdf->shift;
390: Vec V, V0;
392: PetscFunctionBegin;
393: PetscCall(SNESGetDM(snes, &dm));
394: PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
396: /* J,P = Jacobian(t,X,V) */
397: ts->dm = dm;
398: PetscCall(TSComputeIJacobian(ts, t, X, V, shift, J, P, PETSC_FALSE));
399: ts->dm = dmsave;
401: PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: static PetscErrorCode TSReset_BDF(TS ts)
406: {
407: TS_BDF *bdf = (TS_BDF *)ts->data;
408: size_t i, n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
410: PetscFunctionBegin;
411: for (i = 0; i < n; i++) {
412: PetscCall(VecDestroy(&bdf->work[i]));
413: PetscCall(VecDestroy(&bdf->tvwork[i]));
414: }
415: PetscCall(VecDestroy(&bdf->vec_dot));
416: PetscCall(VecDestroy(&bdf->vec_wrk));
417: PetscCall(VecDestroy(&bdf->vec_lte));
418: if (ts->dm) PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode TSDestroy_BDF(TS ts)
423: {
424: PetscFunctionBegin;
425: PetscCall(TSReset_BDF(ts));
426: PetscCall(PetscFree(ts->data));
427: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", NULL));
428: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", NULL));
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: static PetscErrorCode TSSetUp_BDF(TS ts)
433: {
434: TS_BDF *bdf = (TS_BDF *)ts->data;
435: size_t n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
436: PetscReal low, high, two = 2;
437: PetscInt cnt = 0;
439: PetscFunctionBegin;
440: PetscCall(TSHasTransientVariable(ts, &bdf->transientvar));
441: for (size_t i = 0; i < n; i++) {
442: if (!bdf->work[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->work[i]));
443: else cnt++;
444: if (i && bdf->transientvar && !bdf->tvwork[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->tvwork[i]));
445: }
446: if (!cnt) bdf->k = bdf->n = 0;
447: PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_dot));
448: PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_wrk));
449: PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_lte));
450: PetscCall(TSGetDM(ts, &ts->dm));
451: PetscCall(DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
453: PetscCall(TSGetAdapt(ts, &ts->adapt));
454: PetscCall(TSAdaptCandidatesClear(ts->adapt));
455: PetscCall(TSAdaptGetClip(ts->adapt, &low, &high));
456: PetscCall(TSAdaptSetClip(ts->adapt, low, PetscMin(high, two)));
458: PetscCall(TSGetSNES(ts, &ts->snes));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: static PetscErrorCode TSSetFromOptions_BDF(TS ts, PetscOptionItems PetscOptionsObject)
463: {
464: TS_BDF *bdf = (TS_BDF *)ts->data;
466: PetscFunctionBegin;
467: PetscOptionsHeadBegin(PetscOptionsObject, "BDF ODE solver options");
468: {
469: PetscBool flg;
470: PetscInt order;
471: PetscCall(TSBDFGetOrder(ts, &order));
472: PetscCall(PetscOptionsInt("-ts_bdf_order", "Order of the BDF method", "TSBDFSetOrder", order, &order, &flg));
473: if (flg) PetscCall(TSBDFSetOrder(ts, order));
474: PetscCall(PetscOptionsBool("-ts_bdf_initial_guess_extrapolate", "Extrapolate the initial guess of the nonlinear solve from previous time steps", "", bdf->extrapolate, &bdf->extrapolate, NULL));
475: }
476: PetscOptionsHeadEnd();
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode TSView_BDF(TS ts, PetscViewer viewer)
481: {
482: TS_BDF *bdf = (TS_BDF *)ts->data;
483: PetscBool isascii;
485: PetscFunctionBegin;
486: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
487: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Order=%" PetscInt_FMT "\n", bdf->order));
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: /* ------------------------------------------------------------ */
493: static PetscErrorCode TSBDFSetOrder_BDF(TS ts, PetscInt order)
494: {
495: TS_BDF *bdf = (TS_BDF *)ts->data;
497: PetscFunctionBegin;
498: if (order == bdf->order) PetscFunctionReturn(PETSC_SUCCESS);
499: PetscCheck(order >= 1 && order <= 6, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "BDF Order %" PetscInt_FMT " not implemented", order);
500: bdf->order = order;
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: static PetscErrorCode TSBDFGetOrder_BDF(TS ts, PetscInt *order)
505: {
506: TS_BDF *bdf = (TS_BDF *)ts->data;
508: PetscFunctionBegin;
509: *order = bdf->order;
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: /* ------------------------------------------------------------ */
515: /*MC
516: TSBDF - DAE solver using BDF methods
518: Level: beginner
520: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSType`
521: M*/
522: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
523: {
524: TS_BDF *bdf;
526: PetscFunctionBegin;
527: ts->ops->reset = TSReset_BDF;
528: ts->ops->destroy = TSDestroy_BDF;
529: ts->ops->view = TSView_BDF;
530: ts->ops->setup = TSSetUp_BDF;
531: ts->ops->setfromoptions = TSSetFromOptions_BDF;
532: ts->ops->step = TSStep_BDF;
533: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
534: ts->ops->interpolate = TSInterpolate_BDF;
535: ts->ops->resizeregister = TSResizeRegister_BDF;
536: ts->ops->snesfunction = SNESTSFormFunction_BDF;
537: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
538: ts->default_adapt_type = TSADAPTBASIC;
540: ts->usessnes = PETSC_TRUE;
542: PetscCall(PetscNew(&bdf));
543: ts->data = (void *)bdf;
545: bdf->extrapolate = PETSC_TRUE;
546: bdf->status = TS_STEP_COMPLETE;
547: for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(bdf->work); i++) bdf->work[i] = bdf->tvwork[i] = NULL;
549: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", TSBDFSetOrder_BDF));
550: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", TSBDFGetOrder_BDF));
551: PetscCall(TSBDFSetOrder(ts, 2));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /* ------------------------------------------------------------ */
557: /*@
558: TSBDFSetOrder - Set the order of the `TSBDF` method
560: Logically Collective
562: Input Parameters:
563: + ts - timestepping context
564: - order - order of the method
566: Options Database Key:
567: . -ts_bdf_order <order> - select the order
569: Level: intermediate
571: .seealso: `TSBDFGetOrder()`, `TS`, `TSBDF`
572: @*/
573: PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order)
574: {
575: PetscFunctionBegin;
578: PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*@
583: TSBDFGetOrder - Get the order of the `TSBDF` method
585: Not Collective
587: Input Parameter:
588: . ts - timestepping context
590: Output Parameter:
591: . order - order of the method
593: Level: intermediate
595: .seealso: `TSBDFSetOrder()`, `TS`, `TSBDF`
596: @*/
597: PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order)
598: {
599: PetscFunctionBegin;
601: PetscAssertPointer(order, 2);
602: PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order));
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }