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: `TS`, `TSCreate()`, `TSSetType()`
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 BDF method

506:   Logically Collective on TS

508:   Input Parameters:
509: +  ts - timestepping context
510: -  order - order of the method

512:   Options Database:
513: .  -ts_bdf_order <order> - select the order

515:   Level: intermediate

517: @*/
518: PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order)
519: {
522:   PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order));
523:   return 0;
524: }

526: /*@
527:   TSBDFGetOrder - Get the order of the BDF method

529:   Not Collective

531:   Input Parameter:
532: .  ts - timestepping context

534:   Output Parameter:
535: .  order - order of the method

537:   Level: intermediate

539: @*/
540: PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order)
541: {
544:   PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order));
545:   return 0;
546: }