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));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
358: {
359: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
360: PetscReal shift = 1 / (0.5 * ts->time_step);
361: Vec Xdot;
362: DM dm, dmsave;
364: PetscFunctionBegin;
365: PetscCall(SNESGetDM(snes, &dm));
366: /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
367: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot));
369: dmsave = ts->dm;
370: ts->dm = dm;
371: PetscCall(TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE));
372: ts->dm = dmsave;
373: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot));
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: 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)
378: {
379: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
381: PetscFunctionBegin;
382: *Sfunc = dg->Sfunc;
383: *Ffunc = dg->Ffunc;
384: *Gfunc = dg->Gfunc;
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: 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)
389: {
390: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
392: PetscFunctionBegin;
393: dg->Sfunc = Sfunc;
394: dg->Ffunc = Ffunc;
395: dg->Gfunc = Gfunc;
396: dg->funcCtx = ctx;
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*MC
401: TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method
403: Level: intermediate
405: Notes:
406: This is the implicit midpoint rule, with an optional term that guarantees the discrete
407: gradient property. This timestepper applies to systems of the form $u_t = S(u) \nabla F(u)$
408: where $S(u)$ is a linear operator, and $F$ is a functional of $u$.
410: .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`
411: M*/
412: PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
413: {
414: TS_DiscGrad *th;
416: PetscFunctionBegin;
417: PetscCall(PetscCitationsRegister(DGCitation, &DGCite));
418: ts->ops->reset = TSReset_DiscGrad;
419: ts->ops->destroy = TSDestroy_DiscGrad;
420: ts->ops->view = TSView_DiscGrad;
421: ts->ops->setfromoptions = TSSetFromOptions_DiscGrad;
422: ts->ops->setup = TSSetUp_DiscGrad;
423: ts->ops->step = TSStep_DiscGrad;
424: ts->ops->interpolate = TSInterpolate_DiscGrad;
425: ts->ops->getstages = TSGetStages_DiscGrad;
426: ts->ops->snesfunction = SNESTSFormFunction_DiscGrad;
427: ts->ops->snesjacobian = SNESTSFormJacobian_DiscGrad;
428: ts->default_adapt_type = TSADAPTNONE;
430: ts->usessnes = PETSC_TRUE;
432: PetscCall(PetscNew(&th));
433: ts->data = (void *)th;
435: th->gonzalez = PETSC_FALSE;
437: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", TSDiscGradGetFormulation_DiscGrad));
438: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", TSDiscGradSetFormulation_DiscGrad));
439: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradIsGonzalez_C", TSDiscGradIsGonzalez_DiscGrad));
440: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradUseGonzalez_C", TSDiscGradUseGonzalez_DiscGrad));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@C
445: TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the
446: formulation $u_t = S \nabla F$ for `TSDISCGRAD`
448: Not Collective
450: Input Parameter:
451: . ts - timestepping context
453: Output Parameters:
454: + Sfunc - constructor for the S matrix from the formulation
455: . Ffunc - functional F from the formulation
456: . Gfunc - constructor for the gradient of F from the formulation
457: - ctx - the user context
459: Calling sequence of `Sfunc`:
460: + ts - the integrator
461: . time - the current time
462: . u - the solution
463: . S - the S-matrix from the formulation
464: - ctx - the user context
466: Calling sequence of `Ffunc`:
467: + ts - the integrator
468: . time - the current time
469: . u - the solution
470: . F - the computed function from the formulation
471: - ctx - the user context
473: Calling sequence of `Gfunc`:
474: + ts - the integrator
475: . time - the current time
476: . u - the solution
477: . G - the gradient of the computed function from the formulation
478: - ctx - the user context
480: Level: intermediate
482: .seealso: [](ch_ts), `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`
483: @*/
484: 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)
485: {
486: PetscFunctionBegin;
488: PetscAssertPointer(Sfunc, 2);
489: PetscAssertPointer(Ffunc, 3);
490: PetscAssertPointer(Gfunc, 4);
491: 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));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@C
496: TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the
497: formulation $u_t = S(u) \nabla F(u)$ for `TSDISCGRAD`
499: Not Collective
501: Input Parameters:
502: + ts - timestepping context
503: . Sfunc - constructor for the S matrix from the formulation
504: . Ffunc - functional F from the formulation
505: . Gfunc - constructor for the gradient of F from the formulation
506: - ctx - optional context for the functions
508: Calling sequence of `Sfunc`:
509: + ts - the integrator
510: . time - the current time
511: . u - the solution
512: . S - the S-matrix from the formulation
513: - ctx - the user context
515: Calling sequence of `Ffunc`:
516: + ts - the integrator
517: . time - the current time
518: . u - the solution
519: . F - the computed function from the formulation
520: - ctx - the user context
522: Calling sequence of `Gfunc`:
523: + ts - the integrator
524: . time - the current time
525: . u - the solution
526: . G - the gradient of the computed function from the formulation
527: - ctx - the user context
529: Level: intermediate
531: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradGetFormulation()`
532: @*/
533: 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)
534: {
535: PetscFunctionBegin;
540: 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));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: TSDiscGradIsGonzalez - Checks flag for whether to use additional conservative terms in
546: discrete gradient formulation for `TSDISCGRAD`
548: Not Collective
550: Input Parameter:
551: . ts - timestepping context
553: Output Parameter:
554: . gonzalez - `PETSC_TRUE` when using the Gonzalez term
556: Level: advanced
558: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradUseGonzalez()`
559: @*/
560: PetscErrorCode TSDiscGradIsGonzalez(TS ts, PetscBool *gonzalez)
561: {
562: PetscFunctionBegin;
564: PetscAssertPointer(gonzalez, 2);
565: PetscUseMethod(ts, "TSDiscGradIsGonzalez_C", (TS, PetscBool *), (ts, gonzalez));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: /*@
570: TSDiscGradUseGonzalez - Sets discrete gradient formulation with or without additional
571: conservative terms.
573: Not Collective
575: Input Parameters:
576: + ts - timestepping context
577: - flg - `PETSC_TRUE` to use the Gonzalez term
579: Options Database Key:
580: . -ts_discgrad_gonzalez <flg> - use the Gonzalez term for the discrete gradient formulation
582: Level: intermediate
584: Notes:
585: Without `flg`, the discrete gradients timestepper is just backwards Euler.
587: .seealso: [](ch_ts), `TSDISCGRAD`
588: @*/
589: PetscErrorCode TSDiscGradUseGonzalez(TS ts, PetscBool flg)
590: {
591: PetscFunctionBegin;
593: PetscTryMethod(ts, "TSDiscGradUseGonzalez_C", (TS, PetscBool), (ts, flg));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }