Actual source code: alpha1.c

      1: /*
      2:   Code for timestepping with implicit generalized-\alpha method
      3:   for first order systems.
      4: */
      5: #include <petsc/private/tsimpl.h>

      7: static PetscBool  cited      = PETSC_FALSE;
      8: static const char citation[] = "@article{Jansen2000,\n"
      9:                                "  title   = {A generalized-$\\alpha$ method for integrating the filtered {N}avier--{S}tokes equations with a stabilized finite element method},\n"
     10:                                "  author  = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert},\n"
     11:                                "  journal = {Computer Methods in Applied Mechanics and Engineering},\n"
     12:                                "  volume  = {190},\n"
     13:                                "  number  = {3--4},\n"
     14:                                "  pages   = {305--319},\n"
     15:                                "  year    = {2000},\n"
     16:                                "  issn    = {0045-7825},\n"
     17:                                "  doi     = {http://dx.doi.org/10.1016/S0045-7825(00)00203-6}\n}\n";

     19: typedef struct {
     20:   PetscReal stage_time;
     21:   PetscReal shift_V;
     22:   PetscReal scale_F;
     23:   Vec       X0, Xa, X1;
     24:   Vec       V0, Va, V1;

     26:   PetscReal Alpha_m;
     27:   PetscReal Alpha_f;
     28:   PetscReal Gamma;
     29:   PetscInt  order;

     31:   Vec vec_sol_prev;
     32:   Vec vec_lte_work;

     34:   TSStepStatus status;
     35: } TS_Alpha;

     37: static PetscErrorCode TSResizeRegister_Alpha(TS ts, PetscBool reg)
     38: {
     39:   TS_Alpha *th = (TS_Alpha *)ts->data;

     41:   PetscFunctionBegin;
     42:   if (reg) {
     43:     PetscCall(TSResizeRegisterVec(ts, "ts:theta:sol_prev", th->vec_sol_prev));
     44:     PetscCall(TSResizeRegisterVec(ts, "ts:theta:X0", th->X0));
     45:   } else {
     46:     PetscCall(TSResizeRetrieveVec(ts, "ts:theta:sol_prev", &th->vec_sol_prev));
     47:     PetscCall(PetscObjectReference((PetscObject)th->vec_sol_prev));
     48:     PetscCall(TSResizeRetrieveVec(ts, "ts:theta:X0", &th->X0));
     49:     PetscCall(PetscObjectReference((PetscObject)th->X0));
     50:   }
     51:   PetscFunctionReturn(PETSC_SUCCESS);
     52: }

     54: static PetscErrorCode TSAlpha_StageTime(TS ts)
     55: {
     56:   TS_Alpha *th      = (TS_Alpha *)ts->data;
     57:   PetscReal t       = ts->ptime;
     58:   PetscReal dt      = ts->time_step;
     59:   PetscReal Alpha_m = th->Alpha_m;
     60:   PetscReal Alpha_f = th->Alpha_f;
     61:   PetscReal Gamma   = th->Gamma;

     63:   PetscFunctionBegin;
     64:   th->stage_time = t + Alpha_f * dt;
     65:   th->shift_V    = Alpha_m / (Alpha_f * Gamma * dt);
     66:   th->scale_F    = 1 / Alpha_f;
     67:   PetscFunctionReturn(PETSC_SUCCESS);
     68: }

     70: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
     71: {
     72:   TS_Alpha *th = (TS_Alpha *)ts->data;
     73:   Vec       X1 = X, V1 = th->V1;
     74:   Vec       Xa = th->Xa, Va = th->Va;
     75:   Vec       X0 = th->X0, V0 = th->V0;
     76:   PetscReal dt      = ts->time_step;
     77:   PetscReal Alpha_m = th->Alpha_m;
     78:   PetscReal Alpha_f = th->Alpha_f;
     79:   PetscReal Gamma   = th->Gamma;

     81:   PetscFunctionBegin;
     82:   /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
     83:   PetscCall(VecWAXPY(V1, -1.0, X0, X1));
     84:   PetscCall(VecAXPBY(V1, 1 - 1 / Gamma, 1 / (Gamma * dt), V0));
     85:   /* Xa = X0 + Alpha_f*(X1-X0) */
     86:   PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
     87:   PetscCall(VecAYPX(Xa, Alpha_f, X0));
     88:   /* Va = V0 + Alpha_m*(V1-V0) */
     89:   PetscCall(VecWAXPY(Va, -1.0, V0, V1));
     90:   PetscCall(VecAYPX(Va, Alpha_m, V0));
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
     95: {
     96:   PetscInt nits, lits;

     98:   PetscFunctionBegin;
     99:   PetscCall(SNESSolve(ts->snes, b, x));
    100:   PetscCall(SNESGetIterationNumber(ts->snes, &nits));
    101:   PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
    102:   ts->snes_its += nits;
    103:   ts->ksp_its += lits;
    104:   PetscFunctionReturn(PETSC_SUCCESS);
    105: }

    107: /*
    108:   Compute a consistent initial state for the generalized-alpha method.
    109:   - Solve two successive backward Euler steps with halved time step.
    110:   - Compute the initial time derivative using backward differences.
    111:   - If using adaptivity, estimate the LTE of the initial step.
    112: */
    113: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
    114: {
    115:   TS_Alpha *th = (TS_Alpha *)ts->data;
    116:   PetscReal time_step;
    117:   PetscReal alpha_m, alpha_f, gamma;
    118:   Vec       X0 = ts->vec_sol, X1, X2 = th->X1;
    119:   PetscBool stageok;

    121:   PetscFunctionBegin;
    122:   PetscCall(VecDuplicate(X0, &X1));

    124:   /* Setup backward Euler with halved time step */
    125:   PetscCall(TSAlphaGetParams(ts, &alpha_m, &alpha_f, &gamma));
    126:   PetscCall(TSAlphaSetParams(ts, 1, 1, 1));
    127:   PetscCall(TSGetTimeStep(ts, &time_step));
    128:   ts->time_step = time_step / 2;
    129:   PetscCall(TSAlpha_StageTime(ts));
    130:   th->stage_time = ts->ptime;
    131:   PetscCall(VecZeroEntries(th->V0));

    133:   /* First BE step, (t0,X0) -> (t1,X1) */
    134:   th->stage_time += ts->time_step;
    135:   PetscCall(VecCopy(X0, th->X0));
    136:   PetscCall(TSPreStage(ts, th->stage_time));
    137:   PetscCall(VecCopy(th->X0, X1));
    138:   PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
    139:   PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
    140:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
    141:   if (!stageok) goto finally;

    143:   /* Second BE step, (t1,X1) -> (t2,X2) */
    144:   th->stage_time += ts->time_step;
    145:   PetscCall(VecCopy(X1, th->X0));
    146:   PetscCall(TSPreStage(ts, th->stage_time));
    147:   PetscCall(VecCopy(th->X0, X2));
    148:   PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
    149:   PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
    150:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X2, &stageok));
    151:   if (!stageok) goto finally;

    153:   /* Compute V0 ~ dX/dt at t0 with backward differences */
    154:   PetscCall(VecZeroEntries(th->V0));
    155:   PetscCall(VecAXPY(th->V0, -3 / time_step, X0));
    156:   PetscCall(VecAXPY(th->V0, +4 / time_step, X1));
    157:   PetscCall(VecAXPY(th->V0, -1 / time_step, X2));

    159:   /* Rough, lower-order estimate LTE of the initial step */
    160:   if (th->vec_lte_work) {
    161:     PetscCall(VecZeroEntries(th->vec_lte_work));
    162:     PetscCall(VecAXPY(th->vec_lte_work, +2, X2));
    163:     PetscCall(VecAXPY(th->vec_lte_work, -4, X1));
    164:     PetscCall(VecAXPY(th->vec_lte_work, +2, X0));
    165:   }

    167: finally:
    168:   /* Revert TSAlpha to the initial state (t0,X0), but retain
    169:      potential time step reduction by factor after failed solve. */
    170:   if (initok) *initok = stageok;
    171:   PetscCall(TSSetTimeStep(ts, 2 * ts->time_step));
    172:   PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
    173:   PetscCall(VecCopy(ts->vec_sol, th->X0));

    175:   PetscCall(VecDestroy(&X1));
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: static PetscErrorCode TSStep_Alpha(TS ts)
    180: {
    181:   TS_Alpha *th         = (TS_Alpha *)ts->data;
    182:   PetscInt  rejections = 0;
    183:   PetscBool stageok, accept = PETSC_TRUE;
    184:   PetscReal next_time_step = ts->time_step;

    186:   PetscFunctionBegin;
    187:   PetscCall(PetscCitationsRegister(citation, &cited));

    189:   if (!ts->steprollback) {
    190:     if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
    191:     PetscCall(VecCopy(ts->vec_sol, th->X0));
    192:     PetscCall(VecCopy(th->V1, th->V0));
    193:   }

    195:   th->status = TS_STEP_INCOMPLETE;
    196:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
    197:     if (ts->steprestart) {
    198:       PetscCall(TSAlpha_Restart(ts, &stageok));
    199:       if (!stageok) goto reject_step;
    200:     }

    202:     PetscCall(TSAlpha_StageTime(ts));
    203:     PetscCall(VecCopy(th->X0, th->X1));
    204:     PetscCall(TSPreStage(ts, th->stage_time));
    205:     PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
    206:     PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
    207:     PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
    208:     if (!stageok) goto reject_step;

    210:     th->status = TS_STEP_PENDING;
    211:     PetscCall(VecCopy(th->X1, ts->vec_sol));
    212:     PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
    213:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
    214:     if (!accept) {
    215:       PetscCall(VecCopy(th->X0, ts->vec_sol));
    216:       ts->time_step = next_time_step;
    217:       goto reject_step;
    218:     }

    220:     ts->ptime += ts->time_step;
    221:     ts->time_step = next_time_step;
    222:     break;

    224:   reject_step:
    225:     ts->reject++;
    226:     accept = PETSC_FALSE;
    227:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
    228:       ts->reason = TS_DIVERGED_STEP_REJECTED;
    229:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
    230:     }
    231:   }
    232:   PetscFunctionReturn(PETSC_SUCCESS);
    233: }

    235: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
    236: {
    237:   TS_Alpha *th = (TS_Alpha *)ts->data;
    238:   Vec       X  = th->X1;           /* X = solution */
    239:   Vec       Y  = th->vec_lte_work; /* Y = X + LTE  */
    240:   PetscReal wltea, wlter;

    242:   PetscFunctionBegin;
    243:   if (!th->vec_sol_prev) {
    244:     *wlte = -1;
    245:     PetscFunctionReturn(PETSC_SUCCESS);
    246:   }
    247:   if (!th->vec_lte_work) {
    248:     *wlte = -1;
    249:     PetscFunctionReturn(PETSC_SUCCESS);
    250:   }
    251:   if (ts->steprestart) {
    252:     /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
    253:     PetscCall(VecAXPY(Y, 1, X));
    254:   } else {
    255:     /* Compute LTE using backward differences with non-constant time step */
    256:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
    257:     PetscReal   a = 1 + h_prev / h;
    258:     PetscScalar scal[3];
    259:     Vec         vecs[3];
    260:     scal[0] = +1 / a;
    261:     scal[1] = -1 / (a - 1);
    262:     scal[2] = +1 / (a * (a - 1));
    263:     vecs[0] = th->X1;
    264:     vecs[1] = th->X0;
    265:     vecs[2] = th->vec_sol_prev;
    266:     PetscCall(VecCopy(X, Y));
    267:     PetscCall(VecMAXPY(Y, 3, scal, vecs));
    268:   }
    269:   PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
    270:   if (order) *order = 2;
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
    275: {
    276:   TS_Alpha *th    = (TS_Alpha *)ts->data;
    277:   PetscReal dt    = t - ts->ptime;
    278:   PetscReal Gamma = th->Gamma;

    280:   PetscFunctionBegin;
    281:   PetscCall(VecWAXPY(th->V1, -1.0, th->X0, ts->vec_sol));
    282:   PetscCall(VecAXPBY(th->V1, 1 - 1 / Gamma, 1 / (Gamma * ts->time_step), th->V0));
    283:   PetscCall(VecCopy(ts->vec_sol, X));
    284:   /* X = X + Gamma*dT*V1 */
    285:   PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
    286:   /* X = X + (1-Gamma)*dT*V0 */
    287:   PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    289: }

    291: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
    292: {
    293:   TS_Alpha *th = (TS_Alpha *)ts->data;
    294:   PetscReal ta = th->stage_time;
    295:   Vec       Xa = th->Xa, Va = th->Va;

    297:   PetscFunctionBegin;
    298:   PetscCall(TSAlpha_StageVecs(ts, X));
    299:   /* F = Function(ta,Xa,Va) */
    300:   PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
    301:   PetscCall(VecScale(F, th->scale_F));
    302:   PetscFunctionReturn(PETSC_SUCCESS);
    303: }

    305: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
    306: {
    307:   TS_Alpha *th = (TS_Alpha *)ts->data;
    308:   PetscReal ta = th->stage_time;
    309:   Vec       Xa = th->Xa, Va = th->Va;
    310:   PetscReal dVdX = th->shift_V;

    312:   PetscFunctionBegin;
    313:   /* J,P = Jacobian(ta,Xa,Va) */
    314:   PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
    315:   PetscFunctionReturn(PETSC_SUCCESS);
    316: }

    318: static PetscErrorCode TSReset_Alpha(TS ts)
    319: {
    320:   TS_Alpha *th = (TS_Alpha *)ts->data;

    322:   PetscFunctionBegin;
    323:   PetscCall(VecDestroy(&th->X0));
    324:   PetscCall(VecDestroy(&th->Xa));
    325:   PetscCall(VecDestroy(&th->X1));
    326:   PetscCall(VecDestroy(&th->V0));
    327:   PetscCall(VecDestroy(&th->Va));
    328:   PetscCall(VecDestroy(&th->V1));
    329:   PetscCall(VecDestroy(&th->vec_sol_prev));
    330:   PetscCall(VecDestroy(&th->vec_lte_work));
    331:   PetscFunctionReturn(PETSC_SUCCESS);
    332: }

    334: static PetscErrorCode TSDestroy_Alpha(TS ts)
    335: {
    336:   PetscFunctionBegin;
    337:   PetscCall(TSReset_Alpha(ts));
    338:   PetscCall(PetscFree(ts->data));

    340:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
    341:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
    342:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
    343:   PetscFunctionReturn(PETSC_SUCCESS);
    344: }

    346: static PetscErrorCode TSSetUp_Alpha(TS ts)
    347: {
    348:   TS_Alpha *th = (TS_Alpha *)ts->data;
    349:   PetscBool match;

    351:   PetscFunctionBegin;
    352:   if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
    353:   PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
    354:   PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
    355:   PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
    356:   PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
    357:   PetscCall(VecDuplicate(ts->vec_sol, &th->V1));

    359:   PetscCall(TSGetAdapt(ts, &ts->adapt));
    360:   PetscCall(TSAdaptCandidatesClear(ts->adapt));
    361:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
    362:   if (!match) {
    363:     if (!th->vec_sol_prev) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
    364:     if (!th->vec_lte_work) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
    365:   }

    367:   PetscCall(TSGetSNES(ts, &ts->snes));
    368:   PetscFunctionReturn(PETSC_SUCCESS);
    369: }

    371: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems PetscOptionsObject)
    372: {
    373:   TS_Alpha *th = (TS_Alpha *)ts->data;

    375:   PetscFunctionBegin;
    376:   PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
    377:   {
    378:     PetscBool flg;
    379:     PetscReal radius = 1;
    380:     PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
    381:     if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
    382:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
    383:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
    384:     PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
    385:     PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
    386:   }
    387:   PetscOptionsHeadEnd();
    388:   PetscFunctionReturn(PETSC_SUCCESS);
    389: }

    391: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
    392: {
    393:   TS_Alpha *th = (TS_Alpha *)ts->data;
    394:   PetscBool isascii;

    396:   PetscFunctionBegin;
    397:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    398:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  Alpha_m=%g, Alpha_f=%g, Gamma=%g\n", (double)th->Alpha_m, (double)th->Alpha_f, (double)th->Gamma));
    399:   PetscFunctionReturn(PETSC_SUCCESS);
    400: }

    402: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
    403: {
    404:   PetscReal alpha_m, alpha_f, gamma;

    406:   PetscFunctionBegin;
    407:   PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
    408:   alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
    409:   alpha_f = 1 / (1 + radius);
    410:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
    411:   PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
    412:   PetscFunctionReturn(PETSC_SUCCESS);
    413: }

    415: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
    416: {
    417:   TS_Alpha *th  = (TS_Alpha *)ts->data;
    418:   PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
    419:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

    421:   PetscFunctionBegin;
    422:   th->Alpha_m = alpha_m;
    423:   th->Alpha_f = alpha_f;
    424:   th->Gamma   = gamma;
    425:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }

    429: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
    430: {
    431:   TS_Alpha *th = (TS_Alpha *)ts->data;

    433:   PetscFunctionBegin;
    434:   if (alpha_m) *alpha_m = th->Alpha_m;
    435:   if (alpha_f) *alpha_f = th->Alpha_f;
    436:   if (gamma) *gamma = th->Gamma;
    437:   PetscFunctionReturn(PETSC_SUCCESS);
    438: }

    440: /*MC
    441:   TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method {cite}`jansen_2000` {cite}`chung1993` for first-order systems

    443:   Level: beginner

    445: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
    446: M*/
    447: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
    448: {
    449:   TS_Alpha *th;

    451:   PetscFunctionBegin;
    452:   ts->ops->reset          = TSReset_Alpha;
    453:   ts->ops->destroy        = TSDestroy_Alpha;
    454:   ts->ops->view           = TSView_Alpha;
    455:   ts->ops->setup          = TSSetUp_Alpha;
    456:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
    457:   ts->ops->step           = TSStep_Alpha;
    458:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
    459:   ts->ops->interpolate    = TSInterpolate_Alpha;
    460:   ts->ops->resizeregister = TSResizeRegister_Alpha;
    461:   ts->ops->snesfunction   = SNESTSFormFunction_Alpha;
    462:   ts->ops->snesjacobian   = SNESTSFormJacobian_Alpha;
    463:   ts->default_adapt_type  = TSADAPTNONE;

    465:   ts->usessnes = PETSC_TRUE;

    467:   PetscCall(PetscNew(&th));
    468:   ts->data = (void *)th;

    470:   th->Alpha_m = 0.5;
    471:   th->Alpha_f = 0.5;
    472:   th->Gamma   = 0.5;
    473:   th->order   = 2;

    475:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
    476:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
    477:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
    478:   PetscFunctionReturn(PETSC_SUCCESS);
    479: }

    481: /*@
    482:   TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
    483:   (i.e. high-frequency numerical damping)

    485:   Logically Collective

    487:   Input Parameters:
    488: + ts     - timestepping context
    489: - radius - the desired spectral radius

    491:   Options Database Key:
    492: . -ts_alpha_radius <radius> - set alpha radius

    494:   Level: intermediate

    496:   Notes:
    497:   The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
    498:   be computed in terms of a specified spectral radius $\rho$ in [0, 1] for infinite time step
    499:   in order to control high-frequency numerical damping\:

    501:   $$
    502:   \begin{align*}
    503:   \alpha_m = 0.5*(3-\rho)/(1+\rho) \\
    504:   \alpha_f = 1/(1+\rho)
    505:   \end{align*}
    506:   $$

    508: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
    509: @*/
    510: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
    511: {
    512:   PetscFunctionBegin;
    515:   PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
    516:   PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
    517:   PetscFunctionReturn(PETSC_SUCCESS);
    518: }

    520: /*@
    521:   TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`

    523:   Logically Collective

    525:   Input Parameters:
    526: + ts      - timestepping context
    527: . alpha_m - algorithmic parameter
    528: . alpha_f - algorithmic parameter
    529: - gamma   - algorithmic parameter

    531:   Options Database Keys:
    532: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
    533: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
    534: - -ts_alpha_gamma   <gamma>   - set gamma

    536:   Level: advanced

    538:   Note:
    539:   Second-order accuracy can be obtained so long as\:  $\gamma = 0.5 + \alpha_m - \alpha_f$

    541:   Unconditional stability requires\: $\alpha_m >= \alpha_f >= 0.5$

    543:   Backward Euler method is recovered with\: $\alpha_m = \alpha_f = \gamma = 1$

    545:   Use of this function is normally only required to hack `TSALPHA` to use a modified
    546:   integration scheme. Users should call `TSAlphaSetRadius()` to set the desired spectral radius
    547:   of the methods (i.e. high-frequency damping) in order so select optimal values for these
    548:   parameters.

    550: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
    551: @*/
    552: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
    553: {
    554:   PetscFunctionBegin;
    559:   PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
    560:   PetscFunctionReturn(PETSC_SUCCESS);
    561: }

    563: /*@
    564:   TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`

    566:   Not Collective

    568:   Input Parameter:
    569: . ts - timestepping context

    571:   Output Parameters:
    572: + alpha_m - algorithmic parameter
    573: . alpha_f - algorithmic parameter
    574: - gamma   - algorithmic parameter

    576:   Level: advanced

    578:   Note:
    579:   Use of this function is normally only required to hack `TSALPHA` to use a modified
    580:   integration scheme. Users should call `TSAlphaSetRadius()` to set the high-frequency damping
    581:   (i.e. spectral radius of the method) in order so select optimal values for these parameters.

    583: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
    584: @*/
    585: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
    586: {
    587:   PetscFunctionBegin;
    589:   if (alpha_m) PetscAssertPointer(alpha_m, 2);
    590:   if (alpha_f) PetscAssertPointer(alpha_f, 3);
    591:   if (gamma) PetscAssertPointer(gamma, 4);
    592:   PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
    593:   PetscFunctionReturn(PETSC_SUCCESS);
    594: }