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: PetscFunctionBegin;
59: if (dm && dm != ts->dm) {
60: PetscCall(DMGetNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot));
61: PetscCall(DMGetNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot));
62: } else {
63: *Xdot = bdf->vec_dot;
64: *Ydot = bdf->vec_wrk;
65: }
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode TSBDF_RestoreVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
70: {
71: TS_BDF *bdf = (TS_BDF *)ts->data;
73: PetscFunctionBegin;
74: if (dm && dm != ts->dm) {
75: PetscCall(DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot));
76: PetscCall(DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot));
77: } else {
78: PetscCheck(*Xdot == bdf->vec_dot, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Vec does not match the cache");
79: PetscCheck(*Ydot == bdf->vec_wrk, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Vec does not match the cache");
80: *Xdot = NULL;
81: *Ydot = NULL;
82: }
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine, DM coarse, void *ctx)
87: {
88: PetscFunctionBegin;
89: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
99: PetscCall(TSBDF_GetVecs(ts, fine, &Xdot, &Ydot));
100: PetscCall(TSBDF_GetVecs(ts, coarse, &Xdot_c, &Ydot_c));
102: PetscCall(MatRestrict(restrct, Ydot, Ydot_c));
103: PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c));
105: PetscCall(TSBDF_RestoreVecs(ts, fine, &Xdot, &Ydot));
106: PetscCall(TSBDF_RestoreVecs(ts, coarse, &Xdot_c, &Ydot_c));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode TSBDF_Advance(TS ts, PetscReal t, Vec X)
111: {
112: TS_BDF *bdf = (TS_BDF *)ts->data;
113: PetscInt i, n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
114: Vec tail = bdf->work[n - 1], tvtail = bdf->tvwork[n - 1];
116: PetscFunctionBegin;
117: for (i = n - 1; i >= 2; i--) {
118: bdf->time[i] = bdf->time[i - 1];
119: bdf->work[i] = bdf->work[i - 1];
120: bdf->tvwork[i] = bdf->tvwork[i - 1];
121: }
122: bdf->n = PetscMin(bdf->n + 1, n - 1);
123: bdf->time[1] = t;
124: bdf->work[1] = tail;
125: bdf->tvwork[1] = tvtail;
126: PetscCall(VecCopy(X, tail));
127: PetscCall(TSComputeTransientVariable(ts, tail, tvtail));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode TSBDF_VecLTE(TS ts, PetscInt order, Vec lte)
132: {
133: TS_BDF *bdf = (TS_BDF *)ts->data;
134: PetscInt i, n = order + 1;
135: PetscReal *time = bdf->time;
136: Vec *vecs = bdf->work;
137: PetscScalar a[8], b[8], alpha[8];
139: PetscFunctionBegin;
140: LagrangeBasisDers(n + 0, time[0], time, a);
141: a[n] = 0;
142: LagrangeBasisDers(n + 1, time[0], time, b);
143: for (i = 0; i < n + 1; i++) alpha[i] = (a[i] - b[i]) / a[0];
144: PetscCall(VecZeroEntries(lte));
145: PetscCall(VecMAXPY(lte, n + 1, alpha, vecs));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode TSBDF_Extrapolate(TS ts, PetscInt order, PetscReal t, Vec X)
150: {
151: TS_BDF *bdf = (TS_BDF *)ts->data;
152: PetscInt n = order + 1;
153: PetscReal *time = bdf->time + 1;
154: Vec *vecs = bdf->work + 1;
155: PetscScalar alpha[7];
157: PetscFunctionBegin;
158: n = PetscMin(n, bdf->n);
159: LagrangeBasisVals(n, t, time, alpha);
160: PetscCall(VecZeroEntries(X));
161: PetscCall(VecMAXPY(X, n, alpha, vecs));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode TSBDF_Interpolate(TS ts, PetscInt order, PetscReal t, Vec X)
166: {
167: TS_BDF *bdf = (TS_BDF *)ts->data;
168: PetscInt n = order + 1;
169: PetscReal *time = bdf->time;
170: Vec *vecs = bdf->work;
171: PetscScalar alpha[7];
173: PetscFunctionBegin;
174: LagrangeBasisVals(n, t, time, alpha);
175: PetscCall(VecZeroEntries(X));
176: PetscCall(VecMAXPY(X, n, alpha, vecs));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /* Compute the affine term V0 such that Xdot = shift*X + V0.
181: *
182: * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
183: */
184: static PetscErrorCode TSBDF_PreSolve(TS ts)
185: {
186: TS_BDF *bdf = (TS_BDF *)ts->data;
187: PetscInt i, n = PetscMax(bdf->k, 1) + 1;
188: Vec V, V0;
189: Vec vecs[7];
190: PetscScalar alpha[7];
192: PetscFunctionBegin;
193: PetscCall(TSBDF_GetVecs(ts, NULL, &V, &V0));
194: LagrangeBasisDers(n, bdf->time[0], bdf->time, alpha);
195: for (i = 1; i < n; i++) vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
196: PetscCall(VecZeroEntries(V0));
197: PetscCall(VecMAXPY(V0, n - 1, alpha + 1, vecs + 1));
198: bdf->shift = PetscRealPart(alpha[0]);
199: PetscCall(TSBDF_RestoreVecs(ts, NULL, &V, &V0));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode TSBDF_SNESSolve(TS ts, Vec b, Vec x)
204: {
205: PetscInt nits, lits;
207: PetscFunctionBegin;
208: PetscCall(TSBDF_PreSolve(ts));
209: PetscCall(SNESSolve(ts->snes, b, x));
210: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
211: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
212: ts->snes_its += nits;
213: ts->ksp_its += lits;
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode TSBDF_Restart(TS ts, PetscBool *accept)
218: {
219: TS_BDF *bdf = (TS_BDF *)ts->data;
221: PetscFunctionBegin;
222: bdf->k = 1;
223: bdf->n = 0;
224: PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol));
225: if (bdf->order == 1) {
226: *accept = PETSC_TRUE;
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
229: bdf->time[0] = ts->ptime + ts->time_step / 2;
230: PetscCall(VecCopy(bdf->work[1], bdf->work[0]));
231: PetscCall(TSPreStage(ts, bdf->time[0]));
232: PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
233: PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
234: PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], accept));
235: if (!*accept) PetscFunctionReturn(PETSC_SUCCESS);
237: bdf->k = PetscMin(2, bdf->order);
238: bdf->n++;
239: PetscCall(VecCopy(bdf->work[0], bdf->work[2]));
240: bdf->time[2] = bdf->time[0];
241: PetscCall(TSComputeTransientVariable(ts, bdf->work[2], bdf->tvwork[2]));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
247: static PetscErrorCode TSStep_BDF(TS ts)
248: {
249: TS_BDF *bdf = (TS_BDF *)ts->data;
250: PetscInt rejections = 0;
251: PetscBool stageok, accept = PETSC_TRUE;
252: PetscReal next_time_step = ts->time_step;
254: PetscFunctionBegin;
255: PetscCall(PetscCitationsRegister(citation, &cited));
257: if (!ts->steprollback && !ts->steprestart) {
258: bdf->k = PetscMin(bdf->k + 1, bdf->order);
259: PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol));
260: }
262: bdf->status = TS_STEP_INCOMPLETE;
263: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
264: if (ts->steprestart) {
265: PetscCall(TSBDF_Restart(ts, &stageok));
266: if (!stageok) goto reject_step;
267: }
269: bdf->time[0] = ts->ptime + ts->time_step;
270: PetscCall(TSBDF_Extrapolate(ts, bdf->k - (accept ? 0 : 1), bdf->time[0], bdf->work[0]));
271: PetscCall(TSPreStage(ts, bdf->time[0]));
272: PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
273: PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
274: PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], &stageok));
275: if (!stageok) goto reject_step;
277: bdf->status = TS_STEP_PENDING;
278: PetscCall(TSAdaptCandidatesClear(ts->adapt));
279: PetscCall(TSAdaptCandidateAdd(ts->adapt, BDF_SchemeName[bdf->k], bdf->k, 1, 1.0, 1.0, PETSC_TRUE));
280: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
281: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
282: if (!accept) {
283: ts->time_step = next_time_step;
284: goto reject_step;
285: }
287: PetscCall(VecCopy(bdf->work[0], ts->vec_sol));
288: ts->ptime += ts->time_step;
289: ts->time_step = next_time_step;
290: break;
292: reject_step:
293: ts->reject++;
294: accept = PETSC_FALSE;
295: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
296: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
297: ts->reason = TS_DIVERGED_STEP_REJECTED;
298: }
299: }
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode TSInterpolate_BDF(TS ts, PetscReal t, Vec X)
304: {
305: TS_BDF *bdf = (TS_BDF *)ts->data;
307: PetscFunctionBegin;
308: PetscCall(TSBDF_Interpolate(ts, bdf->k, t, X));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
313: {
314: TS_BDF *bdf = (TS_BDF *)ts->data;
315: PetscInt k = bdf->k;
316: PetscReal wltea, wlter;
317: Vec X = bdf->work[0], Y = bdf->vec_lte;
319: PetscFunctionBegin;
320: k = PetscMin(k, bdf->n - 1);
321: if (k == 0) {
322: *wlte = -1;
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
325: PetscCall(TSBDF_VecLTE(ts, k, Y));
326: PetscCall(VecAXPY(Y, 1, X));
327: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
328: if (order) *order = k + 1;
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode TSResizeRegister_BDF(TS ts, PetscBool reg)
333: {
334: TS_BDF *bdf = (TS_BDF *)ts->data;
335: const char *names[] = {"", "ts:bdf:1", "ts:bdf:2", "ts:bdf:3", "ts:bdf:4", "ts:bdf:5", "ts:bdf:6", "ts:bdf:7"};
336: PetscInt i, maxn = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
338: PetscFunctionBegin;
339: PetscAssert(maxn == 8, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "names need to be redefined");
340: if (reg) {
341: for (i = 1; i < PetscMin(bdf->n + 1, maxn); i++) { PetscCall(TSResizeRegisterVec(ts, names[i], bdf->work[i])); }
342: } else {
343: for (i = 1; i < maxn; i++) {
344: PetscCall(TSResizeRetrieveVec(ts, names[i], &bdf->work[i]));
345: if (!bdf->work[i]) break;
346: PetscCall(PetscObjectReference((PetscObject)bdf->work[i]));
347: if (bdf->transientvar) {
348: PetscCall(VecDuplicate(bdf->work[i], &bdf->tvwork[i]));
349: PetscCall(TSComputeTransientVariable(ts, bdf->work[i], bdf->tvwork[i]));
350: }
351: }
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes, Vec X, Vec F, TS ts)
357: {
358: TS_BDF *bdf = (TS_BDF *)ts->data;
359: DM dm, dmsave = ts->dm;
360: PetscReal t = bdf->time[0];
361: PetscReal shift = bdf->shift;
362: Vec V, V0;
364: PetscFunctionBegin;
365: PetscCall(SNESGetDM(snes, &dm));
366: PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
367: if (bdf->transientvar) { /* shift*C(X) + V0 */
368: PetscCall(TSComputeTransientVariable(ts, X, V));
369: PetscCall(VecAYPX(V, shift, V0));
370: } else { /* shift*X + V0 */
371: PetscCall(VecWAXPY(V, shift, X, V0));
372: }
374: /* F = Function(t,X,V) */
375: ts->dm = dm;
376: PetscCall(TSComputeIFunction(ts, t, X, V, F, PETSC_FALSE));
377: ts->dm = dmsave;
379: PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes, Vec X, Mat J, Mat P, TS ts)
384: {
385: TS_BDF *bdf = (TS_BDF *)ts->data;
386: DM dm, dmsave = ts->dm;
387: PetscReal t = bdf->time[0];
388: PetscReal shift = bdf->shift;
389: Vec V, V0;
391: PetscFunctionBegin;
392: PetscCall(SNESGetDM(snes, &dm));
393: PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
395: /* J,P = Jacobian(t,X,V) */
396: ts->dm = dm;
397: PetscCall(TSComputeIJacobian(ts, t, X, V, shift, J, P, PETSC_FALSE));
398: ts->dm = dmsave;
400: PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode TSReset_BDF(TS ts)
405: {
406: TS_BDF *bdf = (TS_BDF *)ts->data;
407: size_t i, n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
409: PetscFunctionBegin;
410: for (i = 0; i < n; i++) {
411: PetscCall(VecDestroy(&bdf->work[i]));
412: PetscCall(VecDestroy(&bdf->tvwork[i]));
413: }
414: PetscCall(VecDestroy(&bdf->vec_dot));
415: PetscCall(VecDestroy(&bdf->vec_wrk));
416: PetscCall(VecDestroy(&bdf->vec_lte));
417: if (ts->dm) PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: static PetscErrorCode TSDestroy_BDF(TS ts)
422: {
423: PetscFunctionBegin;
424: PetscCall(TSReset_BDF(ts));
425: PetscCall(PetscFree(ts->data));
426: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", NULL));
427: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", NULL));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode TSSetUp_BDF(TS ts)
432: {
433: TS_BDF *bdf = (TS_BDF *)ts->data;
434: size_t n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
435: PetscReal low, high, two = 2;
436: PetscInt cnt = 0;
438: PetscFunctionBegin;
439: PetscCall(TSHasTransientVariable(ts, &bdf->transientvar));
440: for (size_t i = 0; i < n; i++) {
441: if (!bdf->work[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->work[i]));
442: else cnt++;
443: if (i && bdf->transientvar && !bdf->tvwork[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->tvwork[i]));
444: }
445: if (!cnt) bdf->k = bdf->n = 0;
446: PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_dot));
447: PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_wrk));
448: PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_lte));
449: PetscCall(TSGetDM(ts, &ts->dm));
450: PetscCall(DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
452: PetscCall(TSGetAdapt(ts, &ts->adapt));
453: PetscCall(TSAdaptCandidatesClear(ts->adapt));
454: PetscCall(TSAdaptGetClip(ts->adapt, &low, &high));
455: PetscCall(TSAdaptSetClip(ts->adapt, low, PetscMin(high, two)));
457: PetscCall(TSGetSNES(ts, &ts->snes));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode TSSetFromOptions_BDF(TS ts, PetscOptionItems PetscOptionsObject)
462: {
463: PetscFunctionBegin;
464: PetscOptionsHeadBegin(PetscOptionsObject, "BDF ODE solver options");
465: {
466: PetscBool flg;
467: PetscInt order;
468: PetscCall(TSBDFGetOrder(ts, &order));
469: PetscCall(PetscOptionsInt("-ts_bdf_order", "Order of the BDF method", "TSBDFSetOrder", order, &order, &flg));
470: if (flg) PetscCall(TSBDFSetOrder(ts, order));
471: }
472: PetscOptionsHeadEnd();
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: static PetscErrorCode TSView_BDF(TS ts, PetscViewer viewer)
477: {
478: TS_BDF *bdf = (TS_BDF *)ts->data;
479: PetscBool iascii;
481: PetscFunctionBegin;
482: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
483: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Order=%" PetscInt_FMT "\n", bdf->order));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: /* ------------------------------------------------------------ */
489: static PetscErrorCode TSBDFSetOrder_BDF(TS ts, PetscInt order)
490: {
491: TS_BDF *bdf = (TS_BDF *)ts->data;
493: PetscFunctionBegin;
494: if (order == bdf->order) PetscFunctionReturn(PETSC_SUCCESS);
495: PetscCheck(order >= 1 && order <= 6, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "BDF Order %" PetscInt_FMT " not implemented", order);
496: bdf->order = order;
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: static PetscErrorCode TSBDFGetOrder_BDF(TS ts, PetscInt *order)
501: {
502: TS_BDF *bdf = (TS_BDF *)ts->data;
504: PetscFunctionBegin;
505: *order = bdf->order;
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: /* ------------------------------------------------------------ */
511: /*MC
512: TSBDF - DAE solver using BDF methods
514: Level: beginner
516: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSType`
517: M*/
518: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
519: {
520: TS_BDF *bdf;
522: PetscFunctionBegin;
523: ts->ops->reset = TSReset_BDF;
524: ts->ops->destroy = TSDestroy_BDF;
525: ts->ops->view = TSView_BDF;
526: ts->ops->setup = TSSetUp_BDF;
527: ts->ops->setfromoptions = TSSetFromOptions_BDF;
528: ts->ops->step = TSStep_BDF;
529: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
530: ts->ops->interpolate = TSInterpolate_BDF;
531: ts->ops->resizeregister = TSResizeRegister_BDF;
532: ts->ops->snesfunction = SNESTSFormFunction_BDF;
533: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
534: ts->default_adapt_type = TSADAPTBASIC;
536: ts->usessnes = PETSC_TRUE;
538: PetscCall(PetscNew(&bdf));
539: ts->data = (void *)bdf;
541: bdf->status = TS_STEP_COMPLETE;
542: for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(bdf->work); i++) { bdf->work[i] = bdf->tvwork[i] = NULL; }
544: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", TSBDFSetOrder_BDF));
545: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", TSBDFGetOrder_BDF));
546: PetscCall(TSBDFSetOrder(ts, 2));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: /* ------------------------------------------------------------ */
552: /*@
553: TSBDFSetOrder - Set the order of the `TSBDF` method
555: Logically Collective
557: Input Parameters:
558: + ts - timestepping context
559: - order - order of the method
561: Options Database Key:
562: . -ts_bdf_order <order> - select the order
564: Level: intermediate
566: .seealso: `TSBDFGetOrder()`, `TS`, `TSBDF`
567: @*/
568: PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order)
569: {
570: PetscFunctionBegin;
573: PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: TSBDFGetOrder - Get the order of the `TSBDF` method
580: Not Collective
582: Input Parameter:
583: . ts - timestepping context
585: Output Parameter:
586: . order - order of the method
588: Level: intermediate
590: .seealso: `TSBDFSetOrder()`, `TS`, `TSBDF`
591: @*/
592: PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order)
593: {
594: PetscFunctionBegin;
596: PetscAssertPointer(order, 2);
597: PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }