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));

    226:   bdf->time[0] = ts->ptime + ts->time_step / 2;
    227:   PetscCall(VecCopy(bdf->work[1], bdf->work[0]));
    228:   PetscCall(TSPreStage(ts, bdf->time[0]));
    229:   PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
    230:   PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
    231:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], accept));
    232:   if (!*accept) PetscFunctionReturn(PETSC_SUCCESS);

    234:   bdf->k = PetscMin(2, bdf->order);
    235:   bdf->n++;
    236:   PetscCall(VecCopy(bdf->work[0], bdf->work[2]));
    237:   bdf->time[2] = bdf->time[0];
    238:   PetscCall(TSComputeTransientVariable(ts, bdf->work[2], bdf->tvwork[2]));
    239:   PetscFunctionReturn(PETSC_SUCCESS);
    240: }

    242: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};

    244: static PetscErrorCode TSStep_BDF(TS ts)
    245: {
    246:   TS_BDF   *bdf        = (TS_BDF *)ts->data;
    247:   PetscInt  rejections = 0;
    248:   PetscBool stageok, accept = PETSC_TRUE;
    249:   PetscReal next_time_step = ts->time_step;

    251:   PetscFunctionBegin;
    252:   PetscCall(PetscCitationsRegister(citation, &cited));

    254:   if (!ts->steprollback && !ts->steprestart) {
    255:     bdf->k = PetscMin(bdf->k + 1, bdf->order);
    256:     PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol));
    257:   }

    259:   bdf->status = TS_STEP_INCOMPLETE;
    260:   while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
    261:     if (ts->steprestart) {
    262:       PetscCall(TSBDF_Restart(ts, &stageok));
    263:       if (!stageok) goto reject_step;
    264:     }

    266:     bdf->time[0] = ts->ptime + ts->time_step;
    267:     PetscCall(TSBDF_Extrapolate(ts, bdf->k - (accept ? 0 : 1), bdf->time[0], bdf->work[0]));
    268:     PetscCall(TSPreStage(ts, bdf->time[0]));
    269:     PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0]));
    270:     PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0]));
    271:     PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], &stageok));
    272:     if (!stageok) goto reject_step;

    274:     bdf->status = TS_STEP_PENDING;
    275:     PetscCall(TSAdaptCandidatesClear(ts->adapt));
    276:     PetscCall(TSAdaptCandidateAdd(ts->adapt, BDF_SchemeName[bdf->k], bdf->k, 1, 1.0, 1.0, PETSC_TRUE));
    277:     PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
    278:     bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
    279:     if (!accept) {
    280:       ts->time_step = next_time_step;
    281:       goto reject_step;
    282:     }

    284:     PetscCall(VecCopy(bdf->work[0], ts->vec_sol));
    285:     ts->ptime += ts->time_step;
    286:     ts->time_step = next_time_step;
    287:     break;

    289:   reject_step:
    290:     ts->reject++;
    291:     accept = PETSC_FALSE;
    292:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
    293:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
    294:       ts->reason = TS_DIVERGED_STEP_REJECTED;
    295:     }
    296:   }
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: static PetscErrorCode TSInterpolate_BDF(TS ts, PetscReal t, Vec X)
    301: {
    302:   TS_BDF *bdf = (TS_BDF *)ts->data;

    304:   PetscFunctionBegin;
    305:   PetscCall(TSBDF_Interpolate(ts, bdf->k, t, X));
    306:   PetscFunctionReturn(PETSC_SUCCESS);
    307: }

    309: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
    310: {
    311:   TS_BDF   *bdf = (TS_BDF *)ts->data;
    312:   PetscInt  k   = bdf->k;
    313:   PetscReal wltea, wlter;
    314:   Vec       X = bdf->work[0], Y = bdf->vec_lte;

    316:   PetscFunctionBegin;
    317:   k = PetscMin(k, bdf->n - 1);
    318:   PetscCall(TSBDF_VecLTE(ts, k, Y));
    319:   PetscCall(VecAXPY(Y, 1, X));
    320:   PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
    321:   if (order) *order = k + 1;
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }

    325: static PetscErrorCode TSResizeRegister_BDF(TS ts, PetscBool reg)
    326: {
    327:   TS_BDF     *bdf     = (TS_BDF *)ts->data;
    328:   const char *names[] = {"", "ts:bdf:1", "ts:bdf:2", "ts:bdf:3", "ts:bdf:4", "ts:bdf:5", "ts:bdf:6", "ts:bdf:7"};
    329:   PetscInt    i, maxn = PETSC_STATIC_ARRAY_LENGTH(bdf->work);

    331:   PetscFunctionBegin;
    332:   PetscAssert(maxn == 8, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "names need to be redefined");
    333:   if (reg) {
    334:     for (i = 1; i < PetscMin(bdf->n + 1, maxn); i++) { PetscCall(TSResizeRegisterVec(ts, names[i], bdf->work[i])); }
    335:   } else {
    336:     for (i = 1; i < maxn; i++) {
    337:       PetscCall(TSResizeRetrieveVec(ts, names[i], &bdf->work[i]));
    338:       if (!bdf->work[i]) break;
    339:       PetscCall(PetscObjectReference((PetscObject)bdf->work[i]));
    340:       if (bdf->transientvar) {
    341:         PetscCall(VecDuplicate(bdf->work[i], &bdf->tvwork[i]));
    342:         PetscCall(TSComputeTransientVariable(ts, bdf->work[i], bdf->tvwork[i]));
    343:       }
    344:     }
    345:   }
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes, Vec X, Vec F, TS ts)
    350: {
    351:   TS_BDF   *bdf = (TS_BDF *)ts->data;
    352:   DM        dm, dmsave = ts->dm;
    353:   PetscReal t     = bdf->time[0];
    354:   PetscReal shift = bdf->shift;
    355:   Vec       V, V0;

    357:   PetscFunctionBegin;
    358:   PetscCall(SNESGetDM(snes, &dm));
    359:   PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));
    360:   if (bdf->transientvar) { /* shift*C(X) + V0 */
    361:     PetscCall(TSComputeTransientVariable(ts, X, V));
    362:     PetscCall(VecAYPX(V, shift, V0));
    363:   } else { /* shift*X + V0 */
    364:     PetscCall(VecWAXPY(V, shift, X, V0));
    365:   }

    367:   /* F = Function(t,X,V) */
    368:   ts->dm = dm;
    369:   PetscCall(TSComputeIFunction(ts, t, X, V, F, PETSC_FALSE));
    370:   ts->dm = dmsave;

    372:   PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
    373:   PetscFunctionReturn(PETSC_SUCCESS);
    374: }

    376: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes, Vec X, Mat J, Mat P, TS ts)
    377: {
    378:   TS_BDF   *bdf = (TS_BDF *)ts->data;
    379:   DM        dm, dmsave = ts->dm;
    380:   PetscReal t     = bdf->time[0];
    381:   PetscReal shift = bdf->shift;
    382:   Vec       V, V0;

    384:   PetscFunctionBegin;
    385:   PetscCall(SNESGetDM(snes, &dm));
    386:   PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0));

    388:   /* J,P = Jacobian(t,X,V) */
    389:   ts->dm = dm;
    390:   PetscCall(TSComputeIJacobian(ts, t, X, V, shift, J, P, PETSC_FALSE));
    391:   ts->dm = dmsave;

    393:   PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0));
    394:   PetscFunctionReturn(PETSC_SUCCESS);
    395: }

    397: static PetscErrorCode TSReset_BDF(TS ts)
    398: {
    399:   TS_BDF *bdf = (TS_BDF *)ts->data;
    400:   size_t  i, n = PETSC_STATIC_ARRAY_LENGTH(bdf->work);

    402:   PetscFunctionBegin;
    403:   for (i = 0; i < n; i++) {
    404:     PetscCall(VecDestroy(&bdf->work[i]));
    405:     PetscCall(VecDestroy(&bdf->tvwork[i]));
    406:   }
    407:   PetscCall(VecDestroy(&bdf->vec_dot));
    408:   PetscCall(VecDestroy(&bdf->vec_wrk));
    409:   PetscCall(VecDestroy(&bdf->vec_lte));
    410:   if (ts->dm) PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));
    411:   PetscFunctionReturn(PETSC_SUCCESS);
    412: }

    414: static PetscErrorCode TSDestroy_BDF(TS ts)
    415: {
    416:   PetscFunctionBegin;
    417:   PetscCall(TSReset_BDF(ts));
    418:   PetscCall(PetscFree(ts->data));
    419:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", NULL));
    420:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", NULL));
    421:   PetscFunctionReturn(PETSC_SUCCESS);
    422: }

    424: static PetscErrorCode TSSetUp_BDF(TS ts)
    425: {
    426:   TS_BDF   *bdf = (TS_BDF *)ts->data;
    427:   size_t    n   = PETSC_STATIC_ARRAY_LENGTH(bdf->work);
    428:   PetscReal low, high, two = 2;
    429:   PetscInt  cnt = 0;

    431:   PetscFunctionBegin;
    432:   PetscCall(TSHasTransientVariable(ts, &bdf->transientvar));
    433:   for (size_t i = 0; i < n; i++) {
    434:     if (!bdf->work[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->work[i]));
    435:     else cnt++;
    436:     if (i && bdf->transientvar && !bdf->tvwork[i]) PetscCall(VecDuplicate(ts->vec_sol, &bdf->tvwork[i]));
    437:   }
    438:   if (!cnt) bdf->k = bdf->n = 0;
    439:   PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_dot));
    440:   PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_wrk));
    441:   PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_lte));
    442:   PetscCall(TSGetDM(ts, &ts->dm));
    443:   PetscCall(DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts));

    445:   PetscCall(TSGetAdapt(ts, &ts->adapt));
    446:   PetscCall(TSAdaptCandidatesClear(ts->adapt));
    447:   PetscCall(TSAdaptGetClip(ts->adapt, &low, &high));
    448:   PetscCall(TSAdaptSetClip(ts->adapt, low, PetscMin(high, two)));

    450:   PetscCall(TSGetSNES(ts, &ts->snes));
    451:   PetscFunctionReturn(PETSC_SUCCESS);
    452: }

    454: static PetscErrorCode TSSetFromOptions_BDF(TS ts, PetscOptionItems *PetscOptionsObject)
    455: {
    456:   PetscFunctionBegin;
    457:   PetscOptionsHeadBegin(PetscOptionsObject, "BDF ODE solver options");
    458:   {
    459:     PetscBool flg;
    460:     PetscInt  order;
    461:     PetscCall(TSBDFGetOrder(ts, &order));
    462:     PetscCall(PetscOptionsInt("-ts_bdf_order", "Order of the BDF method", "TSBDFSetOrder", order, &order, &flg));
    463:     if (flg) PetscCall(TSBDFSetOrder(ts, order));
    464:   }
    465:   PetscOptionsHeadEnd();
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: static PetscErrorCode TSView_BDF(TS ts, PetscViewer viewer)
    470: {
    471:   TS_BDF   *bdf = (TS_BDF *)ts->data;
    472:   PetscBool iascii;

    474:   PetscFunctionBegin;
    475:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    476:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  Order=%" PetscInt_FMT "\n", bdf->order));
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }

    480: /* ------------------------------------------------------------ */

    482: static PetscErrorCode TSBDFSetOrder_BDF(TS ts, PetscInt order)
    483: {
    484:   TS_BDF *bdf = (TS_BDF *)ts->data;

    486:   PetscFunctionBegin;
    487:   if (order == bdf->order) PetscFunctionReturn(PETSC_SUCCESS);
    488:   PetscCheck(order >= 1 && order <= 6, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "BDF Order %" PetscInt_FMT " not implemented", order);
    489:   bdf->order = order;
    490:   PetscFunctionReturn(PETSC_SUCCESS);
    491: }

    493: static PetscErrorCode TSBDFGetOrder_BDF(TS ts, PetscInt *order)
    494: {
    495:   TS_BDF *bdf = (TS_BDF *)ts->data;

    497:   PetscFunctionBegin;
    498:   *order = bdf->order;
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: /* ------------------------------------------------------------ */

    504: /*MC
    505:       TSBDF - DAE solver using BDF methods

    507:   Level: beginner

    509: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSType`
    510: M*/
    511: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
    512: {
    513:   TS_BDF *bdf;

    515:   PetscFunctionBegin;
    516:   ts->ops->reset          = TSReset_BDF;
    517:   ts->ops->destroy        = TSDestroy_BDF;
    518:   ts->ops->view           = TSView_BDF;
    519:   ts->ops->setup          = TSSetUp_BDF;
    520:   ts->ops->setfromoptions = TSSetFromOptions_BDF;
    521:   ts->ops->step           = TSStep_BDF;
    522:   ts->ops->evaluatewlte   = TSEvaluateWLTE_BDF;
    523:   ts->ops->interpolate    = TSInterpolate_BDF;
    524:   ts->ops->resizeregister = TSResizeRegister_BDF;
    525:   ts->ops->snesfunction   = SNESTSFormFunction_BDF;
    526:   ts->ops->snesjacobian   = SNESTSFormJacobian_BDF;
    527:   ts->default_adapt_type  = TSADAPTBASIC;

    529:   ts->usessnes = PETSC_TRUE;

    531:   PetscCall(PetscNew(&bdf));
    532:   ts->data = (void *)bdf;

    534:   bdf->status = TS_STEP_COMPLETE;
    535:   for (size_t i = 0; i < PETSC_STATIC_ARRAY_LENGTH(bdf->work); i++) { bdf->work[i] = bdf->tvwork[i] = NULL; }

    537:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", TSBDFSetOrder_BDF));
    538:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", TSBDFGetOrder_BDF));
    539:   PetscCall(TSBDFSetOrder(ts, 2));
    540:   PetscFunctionReturn(PETSC_SUCCESS);
    541: }

    543: /* ------------------------------------------------------------ */

    545: /*@
    546:   TSBDFSetOrder - Set the order of the `TSBDF` method

    548:   Logically Collective

    550:   Input Parameters:
    551: + ts    - timestepping context
    552: - order - order of the method

    554:   Options Database Key:
    555: . -ts_bdf_order <order> - select the order

    557:   Level: intermediate

    559: .seealso: `TSBDFGetOrder()`, `TS`, `TSBDF`
    560: @*/
    561: PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order)
    562: {
    563:   PetscFunctionBegin;
    566:   PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order));
    567:   PetscFunctionReturn(PETSC_SUCCESS);
    568: }

    570: /*@
    571:   TSBDFGetOrder - Get the order of the `TSBDF` method

    573:   Not Collective

    575:   Input Parameter:
    576: . ts - timestepping context

    578:   Output Parameter:
    579: . order - order of the method

    581:   Level: intermediate

    583: .seealso: `TSBDFSetOrder()`, `TS`, `TSBDF`
    584: @*/
    585: PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order)
    586: {
    587:   PetscFunctionBegin;
    589:   PetscAssertPointer(order, 2);
    590:   PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order));
    591:   PetscFunctionReturn(PETSC_SUCCESS);
    592: }