1 /*
2 Code for timestepping with BDF methods
3 */
4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
5 #include <petscdm.h>
6
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";
14
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;
24
25 PetscBool transientvar;
26 PetscBool extrapolate;
27 PetscInt order;
28 TSStepStatus status;
29 } TS_BDF;
30
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 */
LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])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 }
41
LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])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 }
54
TSBDF_GetVecs(TS ts,DM dm,Vec * Xdot,Vec * Ydot)55 static PetscErrorCode TSBDF_GetVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
56 {
57 TS_BDF *bdf = (TS_BDF *)ts->data;
58
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 }
69
TSBDF_RestoreVecs(TS ts,DM dm,Vec * Xdot,Vec * Ydot)70 static PetscErrorCode TSBDF_RestoreVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot)
71 {
72 TS_BDF *bdf = (TS_BDF *)ts->data;
73
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 }
86
DMCoarsenHook_TSBDF(DM fine,DM coarse,PetscCtx ctx)87 static PetscErrorCode DMCoarsenHook_TSBDF(DM fine, DM coarse, PetscCtx ctx)
88 {
89 PetscFunctionBegin;
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92
DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)93 static PetscErrorCode DMRestrictHook_TSBDF(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
94 {
95 TS ts = (TS)ctx;
96 Vec Ydot, Ydot_c;
97 Vec Xdot, Xdot_c;
98
99 PetscFunctionBegin;
100 PetscCall(TSBDF_GetVecs(ts, fine, &Xdot, &Ydot));
101 PetscCall(TSBDF_GetVecs(ts, coarse, &Xdot_c, &Ydot_c));
102
103 PetscCall(MatRestrict(restrct, Ydot, Ydot_c));
104 PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c));
105
106 PetscCall(TSBDF_RestoreVecs(ts, fine, &Xdot, &Ydot));
107 PetscCall(TSBDF_RestoreVecs(ts, coarse, &Xdot_c, &Ydot_c));
108 PetscFunctionReturn(PETSC_SUCCESS);
109 }
110
TSBDF_Advance(TS ts,PetscReal t,Vec X)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];
116
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 }
131
TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)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];
139
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 }
149
TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)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];
157
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 }
165
TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)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];
173
174 PetscFunctionBegin;
175 LagrangeBasisVals(n, t, time, alpha);
176 PetscCall(VecZeroEntries(X));
177 PetscCall(VecMAXPY(X, n, alpha, vecs));
178 PetscFunctionReturn(PETSC_SUCCESS);
179 }
180
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 */
TSBDF_PreSolve(TS ts)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];
192
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 }
203
TSBDF_SNESSolve(TS ts,Vec b,Vec x)204 static PetscErrorCode TSBDF_SNESSolve(TS ts, Vec b, Vec x)
205 {
206 PetscInt nits, lits;
207
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 }
217
TSBDF_Restart(TS ts,PetscBool * accept)218 static PetscErrorCode TSBDF_Restart(TS ts, PetscBool *accept)
219 {
220 TS_BDF *bdf = (TS_BDF *)ts->data;
221
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);
237
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 }
245
246 static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
247
TSStep_BDF(TS ts)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;
254
255 PetscFunctionBegin;
256 PetscCall(PetscCitationsRegister(citation, &cited));
257
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 }
262
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 }
269
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 else if (!accept) PetscCall(VecCopy(ts->vec_sol, bdf->work[0]));
273 PetscCall(TSPreStage(ts, bdf->time[0]));
274 PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
275 PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
276 PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], &stageok));
277 if (!stageok) goto reject_step;
278
279 bdf->status = TS_STEP_PENDING;
280 PetscCall(TSAdaptCandidatesClear(ts->adapt));
281 PetscCall(TSAdaptCandidateAdd(ts->adapt, BDF_SchemeName[bdf->k], bdf->k, 1, 1.0, 1.0, PETSC_TRUE));
282 PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
283 bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
284 if (!accept) {
285 ts->time_step = next_time_step;
286 goto reject_step;
287 }
288
289 PetscCall(VecCopy(bdf->work[0], ts->vec_sol));
290 ts->ptime += ts->time_step;
291 ts->time_step = next_time_step;
292 break;
293
294 reject_step:
295 ts->reject++;
296 accept = PETSC_FALSE;
297 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
298 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
299 ts->reason = TS_DIVERGED_STEP_REJECTED;
300 }
301 }
302 PetscFunctionReturn(PETSC_SUCCESS);
303 }
304
TSInterpolate_BDF(TS ts,PetscReal t,Vec X)305 static PetscErrorCode TSInterpolate_BDF(TS ts, PetscReal t, Vec X)
306 {
307 TS_BDF *bdf = (TS_BDF *)ts->data;
308
309 PetscFunctionBegin;
310 PetscCall(TSBDF_Interpolate(ts, bdf->k, t, X));
311 PetscFunctionReturn(PETSC_SUCCESS);
312 }
313
TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt * order,PetscReal * wlte)314 static PetscErrorCode TSEvaluateWLTE_BDF(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
315 {
316 TS_BDF *bdf = (TS_BDF *)ts->data;
317 PetscInt k = bdf->k;
318 PetscReal wltea, wlter;
319 Vec X = bdf->work[0], Y = bdf->vec_lte;
320
321 PetscFunctionBegin;
322 k = PetscMin(k, bdf->n - 1);
323 if (k == 0) {
324 *wlte = -1;
325 PetscFunctionReturn(PETSC_SUCCESS);
326 }
327 PetscCall(TSBDF_VecLTE(ts, k, Y));
328 PetscCall(VecAXPY(Y, 1, X));
329 PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
330 if (order) *order = k + 1;
331 PetscFunctionReturn(PETSC_SUCCESS);
332 }
333
TSResizeRegister_BDF(TS ts,PetscBool reg)334 static PetscErrorCode TSResizeRegister_BDF(TS ts, PetscBool reg)
335 {
336 TS_BDF *bdf = (TS_BDF *)ts->data;
337 const char *names[] = {"", "ts:bdf:1", "ts:bdf:2", "ts:bdf:3", "ts:bdf:4", "ts:bdf:5", "ts:bdf:6", "ts:bdf:7"};
338 PetscInt i, maxn = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
339
340 PetscFunctionBegin;
341 PetscAssert(maxn == 8, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "names need to be redefined");
342 if (reg) {
343 for (i = 1; i < PetscMin(bdf->n + 1, maxn); i++) PetscCall(TSResizeRegisterVec(ts, names[i], bdf->work[i]));
344 } else {
345 for (i = 1; i < maxn; i++) {
346 PetscCall(TSResizeRetrieveVec(ts, names[i], &bdf->work[i]));
347 if (!bdf->work[i]) break;
348 PetscCall(PetscObjectReference((PetscObject)bdf->work[i]));
349 if (bdf->transientvar) {
350 PetscCall(VecDuplicate(bdf->work[i], &bdf->tvwork[i]));
351 PetscCall(TSComputeTransientVariable(ts, bdf->work[i], bdf->tvwork[i]));
352 }
353 }
354 }
355 PetscFunctionReturn(PETSC_SUCCESS);
356 }
357
SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts)358 static PetscErrorCode SNESTSFormFunction_BDF(SNES snes, Vec X, Vec F, TS ts)
359 {
360 TS_BDF *bdf = (TS_BDF *)ts->data;
361 DM dm, dmsave = ts->dm;
362 PetscReal t = bdf->time[0];
363 PetscReal shift = bdf->shift;
364 Vec V, V0;
365
366 PetscFunctionBegin;
367 PetscCall(SNESGetDM(snes, &dm));
368 PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
369 if (bdf->transientvar) { /* shift*C(X) + V0 */
370 PetscCall(TSComputeTransientVariable(ts, X, V));
371 PetscCall(VecAYPX(V, shift, V0));
372 } else { /* shift*X + V0 */
373 PetscCall(VecWAXPY(V, shift, X, V0));
374 }
375
376 /* F = Function(t,X,V) */
377 ts->dm = dm;
378 PetscCall(TSComputeIFunction(ts, t, X, V, F, PETSC_FALSE));
379 ts->dm = dmsave;
380
381 PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
382 PetscFunctionReturn(PETSC_SUCCESS);
383 }
384
SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts)385 static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes, Vec X, Mat J, Mat P, TS ts)
386 {
387 TS_BDF *bdf = (TS_BDF *)ts->data;
388 DM dm, dmsave = ts->dm;
389 PetscReal t = bdf->time[0];
390 PetscReal shift = bdf->shift;
391 Vec V, V0;
392
393 PetscFunctionBegin;
394 PetscCall(SNESGetDM(snes, &dm));
395 PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
396
397 /* J,P = Jacobian(t,X,V) */
398 ts->dm = dm;
399 PetscCall(TSComputeIJacobian(ts, t, X, V, shift, J, P, PETSC_FALSE));
400 ts->dm = dmsave;
401
402 PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
403 PetscFunctionReturn(PETSC_SUCCESS);
404 }
405
TSReset_BDF(TS ts)406 static PetscErrorCode TSReset_BDF(TS ts)
407 {
408 TS_BDF *bdf = (TS_BDF *)ts->data;
409 size_t i, n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
410
411 PetscFunctionBegin;
412 for (i = 0; i < n; i++) {
413 PetscCall(VecDestroy(&bdf->work[i]));
414 PetscCall(VecDestroy(&bdf->tvwork[i]));
415 }
416 PetscCall(VecDestroy(&bdf->vec_dot));
417 PetscCall(VecDestroy(&bdf->vec_wrk));
418 PetscCall(VecDestroy(&bdf->vec_lte));
419 if (ts->dm) PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
420 PetscFunctionReturn(PETSC_SUCCESS);
421 }
422
TSDestroy_BDF(TS ts)423 static PetscErrorCode TSDestroy_BDF(TS ts)
424 {
425 PetscFunctionBegin;
426 PetscCall(TSReset_BDF(ts));
427 PetscCall(PetscFree(ts->data));
428 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", NULL));
429 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", NULL));
430 PetscFunctionReturn(PETSC_SUCCESS);
431 }
432
TSSetUp_BDF(TS ts)433 static PetscErrorCode TSSetUp_BDF(TS ts)
434 {
435 TS_BDF *bdf = (TS_BDF *)ts->data;
436 size_t n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
437 PetscReal low, high, two = 2;
438 PetscInt cnt = 0;
439
440 PetscFunctionBegin;
441 PetscCall(TSHasTransientVariable(ts, &bdf->transientvar));
442 for (size_t i = 0; i < n; i++) {
443 if (!bdf->work[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->work[i]));
444 else cnt++;
445 if (i && bdf->transientvar && !bdf->tvwork[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->tvwork[i]));
446 }
447 if (!cnt) bdf->k = bdf->n = 0;
448 PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_dot));
449 PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_wrk));
450 PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_lte));
451 PetscCall(TSGetDM(ts, &ts->dm));
452 PetscCall(DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
453
454 PetscCall(TSGetAdapt(ts, &ts->adapt));
455 PetscCall(TSAdaptCandidatesClear(ts->adapt));
456 PetscCall(TSAdaptGetClip(ts->adapt, &low, &high));
457 PetscCall(TSAdaptSetClip(ts->adapt, low, PetscMin(high, two)));
458
459 PetscCall(TSGetSNES(ts, &ts->snes));
460 PetscFunctionReturn(PETSC_SUCCESS);
461 }
462
TSSetFromOptions_BDF(TS ts,PetscOptionItems PetscOptionsObject)463 static PetscErrorCode TSSetFromOptions_BDF(TS ts, PetscOptionItems PetscOptionsObject)
464 {
465 TS_BDF *bdf = (TS_BDF *)ts->data;
466
467 PetscFunctionBegin;
468 PetscOptionsHeadBegin(PetscOptionsObject, "BDF ODE solver options");
469 {
470 PetscBool flg;
471 PetscInt order;
472 PetscCall(TSBDFGetOrder(ts, &order));
473 PetscCall(PetscOptionsInt("-ts_bdf_order", "Order of the BDF method", "TSBDFSetOrder", order, &order, &flg));
474 if (flg) PetscCall(TSBDFSetOrder(ts, order));
475 PetscCall(PetscOptionsBool("-ts_bdf_initial_guess_extrapolate", "Extrapolate the initial guess of the nonlinear solve from previous time steps", "", bdf->extrapolate, &bdf->extrapolate, NULL));
476 }
477 PetscOptionsHeadEnd();
478 PetscFunctionReturn(PETSC_SUCCESS);
479 }
480
TSView_BDF(TS ts,PetscViewer viewer)481 static PetscErrorCode TSView_BDF(TS ts, PetscViewer viewer)
482 {
483 TS_BDF *bdf = (TS_BDF *)ts->data;
484 PetscBool isascii;
485
486 PetscFunctionBegin;
487 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
488 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Order=%" PetscInt_FMT "\n", bdf->order));
489 PetscFunctionReturn(PETSC_SUCCESS);
490 }
491
TSBDFSetOrder_BDF(TS ts,PetscInt order)492 static PetscErrorCode TSBDFSetOrder_BDF(TS ts, PetscInt order)
493 {
494 TS_BDF *bdf = (TS_BDF *)ts->data;
495
496 PetscFunctionBegin;
497 if (order == bdf->order) PetscFunctionReturn(PETSC_SUCCESS);
498 PetscCheck(order >= 1 && order <= 6, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "BDF Order %" PetscInt_FMT " not implemented", order);
499 bdf->order = order;
500 PetscFunctionReturn(PETSC_SUCCESS);
501 }
502
TSBDFGetOrder_BDF(TS ts,PetscInt * order)503 static PetscErrorCode TSBDFGetOrder_BDF(TS ts, PetscInt *order)
504 {
505 TS_BDF *bdf = (TS_BDF *)ts->data;
506
507 PetscFunctionBegin;
508 *order = bdf->order;
509 PetscFunctionReturn(PETSC_SUCCESS);
510 }
511
512 /*MC
513 TSBDF - DAE solver using implicit backward differentiation formula (BDF) methods suitable for stiff ODEs.
514
515 Level: beginner
516
517 Options Database Keys:
518 + -ts_bdf_order <n> - Order of the BDF method
519 - -ts_bdf_initial_guess_extrapolate <true, false> - Extrapolate the initial guess of the nonlinear solve from previous time steps, defaults to true
520
521 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSType`, `TSBDFSetOrder()`
522 M*/
TSCreate_BDF(TS ts)523 PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
524 {
525 TS_BDF *bdf;
526
527 PetscFunctionBegin;
528 ts->ops->reset = TSReset_BDF;
529 ts->ops->destroy = TSDestroy_BDF;
530 ts->ops->view = TSView_BDF;
531 ts->ops->setup = TSSetUp_BDF;
532 ts->ops->setfromoptions = TSSetFromOptions_BDF;
533 ts->ops->step = TSStep_BDF;
534 ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
535 ts->ops->interpolate = TSInterpolate_BDF;
536 ts->ops->resizeregister = TSResizeRegister_BDF;
537 ts->ops->snesfunction = SNESTSFormFunction_BDF;
538 ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
539 ts->default_adapt_type = TSADAPTBASIC;
540
541 ts->usessnes = PETSC_TRUE;
542
543 PetscCall(PetscNew(&bdf));
544 ts->data = (void *)bdf;
545
546 bdf->extrapolate = PETSC_TRUE;
547 bdf->status = TS_STEP_COMPLETE;
548 for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(bdf->work); i++) bdf->work[i] = bdf->tvwork[i] = NULL;
549
550 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", TSBDFSetOrder_BDF));
551 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", TSBDFGetOrder_BDF));
552 PetscCall(TSBDFSetOrder(ts, 2));
553 PetscFunctionReturn(PETSC_SUCCESS);
554 }
555
556 /*@
557 TSBDFSetOrder - Set the order of the `TSBDF` method
558
559 Logically Collective
560
561 Input Parameters:
562 + ts - timestepping context
563 - order - order of the method
564
565 Options Database Key:
566 . -ts_bdf_order <order> - select the order
567
568 Level: intermediate
569
570 .seealso: `TSBDFGetOrder()`, `TS`, `TSBDF`
571 @*/
TSBDFSetOrder(TS ts,PetscInt order)572 PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order)
573 {
574 PetscFunctionBegin;
575 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
576 PetscValidLogicalCollectiveInt(ts, order, 2);
577 PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order));
578 PetscFunctionReturn(PETSC_SUCCESS);
579 }
580
581 /*@
582 TSBDFGetOrder - Get the order of the `TSBDF` method
583
584 Not Collective
585
586 Input Parameter:
587 . ts - timestepping context
588
589 Output Parameter:
590 . order - order of the method
591
592 Level: intermediate
593
594 .seealso: `TSBDFSetOrder()`, `TS`, `TSBDF`
595 @*/
TSBDFGetOrder(TS ts,PetscInt * order)596 PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order)
597 {
598 PetscFunctionBegin;
599 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
600 PetscAssertPointer(order, 2);
601 PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order));
602 PetscFunctionReturn(PETSC_SUCCESS);
603 }
604