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: }