Actual source code: tsdiscgrad.c

      1: /*
      2:   Code for timestepping with discrete gradient integrators
      3: */
      4: #include <petsc/private/tsimpl.h>
      5: #include <petscdm.h>

      7: PetscBool  DGCite       = PETSC_FALSE;
      8: const char DGCitation[] = "@article{Gonzalez1996,\n"
      9:                           "  title   = {Time integration and discrete Hamiltonian systems},\n"
     10:                           "  author  = {Oscar Gonzalez},\n"
     11:                           "  journal = {Journal of Nonlinear Science},\n"
     12:                           "  volume  = {6},\n"
     13:                           "  pages   = {449--467},\n"
     14:                           "  doi     = {10.1007/978-1-4612-1246-1_10},\n"
     15:                           "  year    = {1996}\n}\n";

     17: typedef struct {
     18:   PetscReal stage_time;
     19:   Vec       X0, X, Xdot;
     20:   void     *funcCtx;
     21:   PetscBool gonzalez;
     22:   PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *);
     23:   PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *);
     24:   PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *);
     25: } TS_DiscGrad;

     27: static PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
     28: {
     29:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

     31:   PetscFunctionBegin;
     32:   if (X0) {
     33:     if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0));
     34:     else *X0 = ts->vec_sol;
     35:   }
     36:   if (Xdot) {
     37:     if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot));
     38:     else *Xdot = dg->Xdot;
     39:   }
     40:   PetscFunctionReturn(PETSC_SUCCESS);
     41: }

     43: static PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
     44: {
     45:   PetscFunctionBegin;
     46:   if (X0) {
     47:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0));
     48:   }
     49:   if (Xdot) {
     50:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot));
     51:   }
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, void *ctx)
     56: {
     57:   PetscFunctionBegin;
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

     61: static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
     62: {
     63:   TS  ts = (TS)ctx;
     64:   Vec X0, Xdot, X0_c, Xdot_c;

     66:   PetscFunctionBegin;
     67:   PetscCall(TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot));
     68:   PetscCall(TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c));
     69:   PetscCall(MatRestrict(restrct, X0, X0_c));
     70:   PetscCall(MatRestrict(restrct, Xdot, Xdot_c));
     71:   PetscCall(VecPointwiseMult(X0_c, rscale, X0_c));
     72:   PetscCall(VecPointwiseMult(Xdot_c, rscale, Xdot_c));
     73:   PetscCall(TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot));
     74:   PetscCall(TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c));
     75:   PetscFunctionReturn(PETSC_SUCCESS);
     76: }

     78: static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, void *ctx)
     79: {
     80:   PetscFunctionBegin;
     81:   PetscFunctionReturn(PETSC_SUCCESS);
     82: }

     84: static PetscErrorCode DMSubDomainRestrictHook_TSDiscGrad(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
     85: {
     86:   TS  ts = (TS)ctx;
     87:   Vec X0, Xdot, X0_sub, Xdot_sub;

     89:   PetscFunctionBegin;
     90:   PetscCall(TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot));
     91:   PetscCall(TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub));

     93:   PetscCall(VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD));
     94:   PetscCall(VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD));

     96:   PetscCall(VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD));
     97:   PetscCall(VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD));

     99:   PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot));
    100:   PetscCall(TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub));
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: static PetscErrorCode TSSetUp_DiscGrad(TS ts)
    105: {
    106:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
    107:   DM           dm;

    109:   PetscFunctionBegin;
    110:   if (!dg->X) PetscCall(VecDuplicate(ts->vec_sol, &dg->X));
    111:   if (!dg->X0) PetscCall(VecDuplicate(ts->vec_sol, &dg->X0));
    112:   if (!dg->Xdot) PetscCall(VecDuplicate(ts->vec_sol, &dg->Xdot));

    114:   PetscCall(TSGetDM(ts, &dm));
    115:   PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts));
    116:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts));
    117:   PetscFunctionReturn(PETSC_SUCCESS);
    118: }

    120: static PetscErrorCode TSSetFromOptions_DiscGrad(TS ts, PetscOptionItems *PetscOptionsObject)
    121: {
    122:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    124:   PetscFunctionBegin;
    125:   PetscOptionsHeadBegin(PetscOptionsObject, "Discrete Gradients ODE solver options");
    126:   {
    127:     PetscCall(PetscOptionsBool("-ts_discgrad_gonzalez", "Use Gonzalez term in discrete gradients formulation", "TSDiscGradUseGonzalez", dg->gonzalez, &dg->gonzalez, NULL));
    128:   }
    129:   PetscOptionsHeadEnd();
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: static PetscErrorCode TSView_DiscGrad(TS ts, PetscViewer viewer)
    134: {
    135:   PetscBool iascii;

    137:   PetscFunctionBegin;
    138:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    139:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  Discrete Gradients\n"));
    140:   PetscFunctionReturn(PETSC_SUCCESS);
    141: }

    143: static PetscErrorCode TSDiscGradIsGonzalez_DiscGrad(TS ts, PetscBool *gonzalez)
    144: {
    145:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    147:   PetscFunctionBegin;
    148:   *gonzalez = dg->gonzalez;
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: static PetscErrorCode TSDiscGradUseGonzalez_DiscGrad(TS ts, PetscBool flg)
    153: {
    154:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    156:   PetscFunctionBegin;
    157:   dg->gonzalez = flg;
    158:   PetscFunctionReturn(PETSC_SUCCESS);
    159: }

    161: static PetscErrorCode TSReset_DiscGrad(TS ts)
    162: {
    163:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    165:   PetscFunctionBegin;
    166:   PetscCall(VecDestroy(&dg->X));
    167:   PetscCall(VecDestroy(&dg->X0));
    168:   PetscCall(VecDestroy(&dg->Xdot));
    169:   PetscFunctionReturn(PETSC_SUCCESS);
    170: }

    172: static PetscErrorCode TSDestroy_DiscGrad(TS ts)
    173: {
    174:   DM dm;

    176:   PetscFunctionBegin;
    177:   PetscCall(TSReset_DiscGrad(ts));
    178:   PetscCall(TSGetDM(ts, &dm));
    179:   if (dm) {
    180:     PetscCall(DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts));
    181:     PetscCall(DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts));
    182:   }
    183:   PetscCall(PetscFree(ts->data));
    184:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", NULL));
    185:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", NULL));
    186:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradIsGonzalez_C", NULL));
    187:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradUseGonzalez_C", NULL));
    188:   PetscFunctionReturn(PETSC_SUCCESS);
    189: }

    191: static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X)
    192: {
    193:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
    194:   PetscReal    dt = t - ts->ptime;

    196:   PetscFunctionBegin;
    197:   PetscCall(VecCopy(ts->vec_sol, dg->X));
    198:   PetscCall(VecWAXPY(X, dt, dg->Xdot, dg->X));
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x)
    203: {
    204:   SNES     snes;
    205:   PetscInt nits, lits;

    207:   PetscFunctionBegin;
    208:   PetscCall(TSGetSNES(ts, &snes));
    209:   PetscCall(SNESSolve(snes, b, x));
    210:   PetscCall(SNESGetIterationNumber(snes, &nits));
    211:   PetscCall(SNESGetLinearSolveIterations(snes, &lits));
    212:   ts->snes_its += nits;
    213:   ts->ksp_its += lits;
    214:   PetscFunctionReturn(PETSC_SUCCESS);
    215: }

    217: static PetscErrorCode TSStep_DiscGrad(TS ts)
    218: {
    219:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
    220:   TSAdapt      adapt;
    221:   TSStepStatus status     = TS_STEP_INCOMPLETE;
    222:   PetscInt     rejections = 0;
    223:   PetscBool    stageok, accept = PETSC_TRUE;
    224:   PetscReal    next_time_step = ts->time_step;

    226:   PetscFunctionBegin;
    227:   PetscCall(TSGetAdapt(ts, &adapt));
    228:   if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, dg->X0));

    230:   while (!ts->reason && status != TS_STEP_COMPLETE) {
    231:     PetscReal shift = 1 / (0.5 * ts->time_step);

    233:     dg->stage_time = ts->ptime + 0.5 * ts->time_step;

    235:     PetscCall(VecCopy(dg->X0, dg->X));
    236:     PetscCall(TSPreStage(ts, dg->stage_time));
    237:     PetscCall(TSDiscGrad_SNESSolve(ts, NULL, dg->X));
    238:     PetscCall(TSPostStage(ts, dg->stage_time, 0, &dg->X));
    239:     PetscCall(TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok));
    240:     if (!stageok) goto reject_step;

    242:     status = TS_STEP_PENDING;
    243:     PetscCall(VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X));
    244:     PetscCall(VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot));
    245:     PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
    246:     status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
    247:     if (!accept) {
    248:       PetscCall(VecCopy(dg->X0, ts->vec_sol));
    249:       ts->time_step = next_time_step;
    250:       goto reject_step;
    251:     }
    252:     ts->ptime += ts->time_step;
    253:     ts->time_step = next_time_step;
    254:     break;

    256:   reject_step:
    257:     ts->reject++;
    258:     accept = PETSC_FALSE;
    259:     if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) {
    260:       ts->reason = TS_DIVERGED_STEP_REJECTED;
    261:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
    262:     }
    263:   }
    264:   PetscFunctionReturn(PETSC_SUCCESS);
    265: }

    267: static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y)
    268: {
    269:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    271:   PetscFunctionBegin;
    272:   if (ns) *ns = 1;
    273:   if (Y) *Y = &dg->X;
    274:   PetscFunctionReturn(PETSC_SUCCESS);
    275: }

    277: /*
    278:   This defines the nonlinear equation that is to be solved with SNES
    279:     G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0
    280: */

    282: /* x = (x+x')/2 */
    283: /* NEED TO CALCULATE x_{n+1} from x and x_{n}*/
    284: static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts)
    285: {
    286:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
    287:   PetscReal    norm, shift = 1 / (0.5 * ts->time_step);
    288:   PetscInt     n;
    289:   Vec          X0, Xdot, Xp, Xdiff;
    290:   Mat          S;
    291:   PetscScalar  F = 0, F0 = 0, Gp;
    292:   Vec          G, SgF;
    293:   DM           dm, dmsave;

    295:   PetscFunctionBegin;
    296:   PetscCall(SNESGetDM(snes, &dm));

    298:   PetscCall(VecDuplicate(y, &Xp));
    299:   PetscCall(VecDuplicate(y, &Xdiff));
    300:   PetscCall(VecDuplicate(y, &SgF));
    301:   PetscCall(VecDuplicate(y, &G));

    303:   PetscCall(VecGetLocalSize(y, &n));
    304:   PetscCall(MatCreate(PETSC_COMM_WORLD, &S));
    305:   PetscCall(MatSetSizes(S, PETSC_DECIDE, PETSC_DECIDE, n, n));
    306:   PetscCall(MatSetFromOptions(S));
    307:   PetscCall(MatSetUp(S));
    308:   PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
    309:   PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));

    311:   PetscCall(TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot));
    312:   PetscCall(VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x)); /* Xdot = shift (x - X0) */

    314:   PetscCall(VecAXPBYPCZ(Xp, -1, 2, 0, X0, x));     /* Xp = 2*x - X0 + (0)*Xmid */
    315:   PetscCall(VecAXPBYPCZ(Xdiff, -1, 1, 0, X0, Xp)); /* Xdiff = xp - X0 + (0)*Xdiff */

    317:   if (dg->gonzalez) {
    318:     PetscCall((*dg->Sfunc)(ts, dg->stage_time, x, S, dg->funcCtx));
    319:     PetscCall((*dg->Ffunc)(ts, dg->stage_time, Xp, &F, dg->funcCtx));
    320:     PetscCall((*dg->Ffunc)(ts, dg->stage_time, X0, &F0, dg->funcCtx));
    321:     PetscCall((*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx));

    323:     /* Adding Extra Gonzalez Term */
    324:     PetscCall(VecDot(Xdiff, G, &Gp));
    325:     PetscCall(VecNorm(Xdiff, NORM_2, &norm));
    326:     if (norm < PETSC_SQRT_MACHINE_EPSILON) {
    327:       Gp = 0;
    328:     } else {
    329:       /* Gp = (1/|xn+1 - xn|^2) * (F(xn+1) - F(xn) - Gp) */
    330:       Gp = (F - F0 - Gp) / PetscSqr(norm);
    331:     }
    332:     PetscCall(VecAXPY(G, Gp, Xdiff));
    333:     PetscCall(MatMult(S, G, SgF)); /* S*gradF */

    335:   } else {
    336:     PetscCall((*dg->Sfunc)(ts, dg->stage_time, x, S, dg->funcCtx));
    337:     PetscCall((*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx));

    339:     PetscCall(MatMult(S, G, SgF)); /* Xdot = S*gradF */
    340:   }
    341:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
    342:   dmsave = ts->dm;
    343:   ts->dm = dm;
    344:   PetscCall(VecAXPBYPCZ(y, 1, -1, 0, Xdot, SgF));
    345:   ts->dm = dmsave;
    346:   PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot));

    348:   PetscCall(VecDestroy(&Xp));
    349:   PetscCall(VecDestroy(&Xdiff));
    350:   PetscCall(VecDestroy(&SgF));
    351:   PetscCall(VecDestroy(&G));
    352:   PetscCall(MatDestroy(&S));
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

    356: static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
    357: {
    358:   TS_DiscGrad *dg    = (TS_DiscGrad *)ts->data;
    359:   PetscReal    shift = 1 / (0.5 * ts->time_step);
    360:   Vec          Xdot;
    361:   DM           dm, dmsave;

    363:   PetscFunctionBegin;
    364:   PetscCall(SNESGetDM(snes, &dm));
    365:   /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
    366:   PetscCall(TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot));

    368:   dmsave = ts->dm;
    369:   ts->dm = dm;
    370:   PetscCall(TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE));
    371:   ts->dm = dmsave;
    372:   PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot));
    373:   PetscFunctionReturn(PETSC_SUCCESS);
    374: }

    376: static PetscErrorCode TSDiscGradGetFormulation_DiscGrad(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx)
    377: {
    378:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    380:   PetscFunctionBegin;
    381:   *Sfunc = dg->Sfunc;
    382:   *Ffunc = dg->Ffunc;
    383:   *Gfunc = dg->Gfunc;
    384:   PetscFunctionReturn(PETSC_SUCCESS);
    385: }

    387: static PetscErrorCode TSDiscGradSetFormulation_DiscGrad(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx)
    388: {
    389:   TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;

    391:   PetscFunctionBegin;
    392:   dg->Sfunc   = Sfunc;
    393:   dg->Ffunc   = Ffunc;
    394:   dg->Gfunc   = Gfunc;
    395:   dg->funcCtx = ctx;
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: /*MC
    400:   TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method

    402:   Level: intermediate

    404:   Notes:
    405:   This is the implicit midpoint rule, with an optional term that guarantees the discrete
    406:   gradient property. This timestepper applies to systems of the form $u_t = S(u) \nabla F(u)$
    407:   where $S(u)$ is a linear operator, and $F$ is a functional of $u$.

    409: .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`
    410: M*/
    411: PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
    412: {
    413:   TS_DiscGrad *th;

    415:   PetscFunctionBegin;
    416:   PetscCall(PetscCitationsRegister(DGCitation, &DGCite));
    417:   ts->ops->reset          = TSReset_DiscGrad;
    418:   ts->ops->destroy        = TSDestroy_DiscGrad;
    419:   ts->ops->view           = TSView_DiscGrad;
    420:   ts->ops->setfromoptions = TSSetFromOptions_DiscGrad;
    421:   ts->ops->setup          = TSSetUp_DiscGrad;
    422:   ts->ops->step           = TSStep_DiscGrad;
    423:   ts->ops->interpolate    = TSInterpolate_DiscGrad;
    424:   ts->ops->getstages      = TSGetStages_DiscGrad;
    425:   ts->ops->snesfunction   = SNESTSFormFunction_DiscGrad;
    426:   ts->ops->snesjacobian   = SNESTSFormJacobian_DiscGrad;
    427:   ts->default_adapt_type  = TSADAPTNONE;

    429:   ts->usessnes = PETSC_TRUE;

    431:   PetscCall(PetscNew(&th));
    432:   ts->data = (void *)th;

    434:   th->gonzalez = PETSC_FALSE;

    436:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", TSDiscGradGetFormulation_DiscGrad));
    437:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", TSDiscGradSetFormulation_DiscGrad));
    438:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradIsGonzalez_C", TSDiscGradIsGonzalez_DiscGrad));
    439:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradUseGonzalez_C", TSDiscGradUseGonzalez_DiscGrad));
    440:   PetscFunctionReturn(PETSC_SUCCESS);
    441: }

    443: /*@C
    444:   TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the
    445:   formulation $u_t = S \nabla F$ for `TSDISCGRAD`

    447:   Not Collective

    449:   Input Parameter:
    450: . ts - timestepping context

    452:   Output Parameters:
    453: + Sfunc - constructor for the S matrix from the formulation
    454: . Ffunc - functional F from the formulation
    455: . Gfunc - constructor for the gradient of F from the formulation
    456: - ctx   - the user context

    458:   Calling sequence of `Sfunc`:
    459: + ts   - the integrator
    460: . time - the current time
    461: . u    - the solution
    462: . S    - the S-matrix from the formulation
    463: - ctx  - the user context

    465:   Calling sequence of `Ffunc`:
    466: + ts   - the integrator
    467: . time - the current time
    468: . u    - the solution
    469: . F    - the computed function from the formulation
    470: - ctx  - the user context

    472:   Calling sequence of `Gfunc`:
    473: + ts   - the integrator
    474: . time - the current time
    475: . u    - the solution
    476: . G    - the gradient of the computed function from the formulation
    477: - ctx  - the user context

    479:   Level: intermediate

    481: .seealso: [](ch_ts), `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`
    482: @*/
    483: PetscErrorCode TSDiscGradGetFormulation(TS ts, PetscErrorCode (**Sfunc)(TS ts, PetscReal time, Vec u, Mat S, void *ctx), PetscErrorCode (**Ffunc)(TS ts, PetscReal time, Vec u, PetscScalar *F, void *ctx), PetscErrorCode (**Gfunc)(TS ts, PetscReal time, Vec u, Vec G, void *ctx), void *ctx)
    484: {
    485:   PetscFunctionBegin;
    487:   PetscAssertPointer(Sfunc, 2);
    488:   PetscAssertPointer(Ffunc, 3);
    489:   PetscAssertPointer(Gfunc, 4);
    490:   PetscUseMethod(ts, "TSDiscGradGetFormulation_C", (TS, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *), void *), (ts, Sfunc, Ffunc, Gfunc, ctx));
    491:   PetscFunctionReturn(PETSC_SUCCESS);
    492: }

    494: /*@C
    495:   TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the
    496:   formulation $u_t = S(u) \nabla F(u)$ for `TSDISCGRAD`

    498:   Not Collective

    500:   Input Parameters:
    501: + ts    - timestepping context
    502: . Sfunc - constructor for the S matrix from the formulation
    503: . Ffunc - functional F from the formulation
    504: . Gfunc - constructor for the gradient of F from the formulation
    505: - ctx   - optional context for the functions

    507:   Calling sequence of `Sfunc`:
    508: + ts   - the integrator
    509: . time - the current time
    510: . u    - the solution
    511: . S    - the S-matrix from the formulation
    512: - ctx  - the user context

    514:   Calling sequence of `Ffunc`:
    515: + ts   - the integrator
    516: . time - the current time
    517: . u    - the solution
    518: . F    - the computed function from the formulation
    519: - ctx  - the user context

    521:   Calling sequence of `Gfunc`:
    522: + ts   - the integrator
    523: . time - the current time
    524: . u    - the solution
    525: . G    - the gradient of the computed function from the formulation
    526: - ctx  - the user context

    528:   Level: intermediate

    530: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradGetFormulation()`
    531: @*/
    532: PetscErrorCode TSDiscGradSetFormulation(TS ts, PetscErrorCode (*Sfunc)(TS ts, PetscReal time, Vec u, Mat S, void *ctx), PetscErrorCode (*Ffunc)(TS ts, PetscReal time, Vec u, PetscScalar *F, void *ctx), PetscErrorCode (*Gfunc)(TS ts, PetscReal time, Vec u, Vec G, void *ctx), void *ctx)
    533: {
    534:   PetscFunctionBegin;
    539:   PetscTryMethod(ts, "TSDiscGradSetFormulation_C", (TS, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *), void *), (ts, Sfunc, Ffunc, Gfunc, ctx));
    540:   PetscFunctionReturn(PETSC_SUCCESS);
    541: }

    543: /*@
    544:   TSDiscGradIsGonzalez - Checks flag for whether to use additional conservative terms in
    545:   discrete gradient formulation for `TSDISCGRAD`

    547:   Not Collective

    549:   Input Parameter:
    550: . ts - timestepping context

    552:   Output Parameter:
    553: . gonzalez - `PETSC_TRUE` when using the Gonzalez term

    555:   Level: advanced

    557: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradUseGonzalez()`
    558: @*/
    559: PetscErrorCode TSDiscGradIsGonzalez(TS ts, PetscBool *gonzalez)
    560: {
    561:   PetscFunctionBegin;
    563:   PetscAssertPointer(gonzalez, 2);
    564:   PetscUseMethod(ts, "TSDiscGradIsGonzalez_C", (TS, PetscBool *), (ts, gonzalez));
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566: }

    568: /*@
    569:   TSDiscGradUseGonzalez - Sets discrete gradient formulation with or without additional
    570:   conservative terms.

    572:   Not Collective

    574:   Input Parameters:
    575: + ts  - timestepping context
    576: - flg - `PETSC_TRUE` to use the Gonzalez term

    578:   Options Database Key:
    579: . -ts_discgrad_gonzalez <flg> - use the Gonzalez term for the discrete gradient formulation

    581:   Level: intermediate

    583:   Notes:
    584:   Without `flg`, the discrete gradients timestepper is just backwards Euler.

    586: .seealso: [](ch_ts), `TSDISCGRAD`
    587: @*/
    588: PetscErrorCode TSDiscGradUseGonzalez(TS ts, PetscBool flg)
    589: {
    590:   PetscFunctionBegin;
    592:   PetscTryMethod(ts, "TSDiscGradUseGonzalez_C", (TS, PetscBool), (ts, flg));
    593:   PetscFunctionReturn(PETSC_SUCCESS);
    594: }