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: const char *DGTypes[] = {"gonzalez", "average", "none", "TSDGType", "DG_", NULL};
19: typedef struct {
20: PetscReal stage_time;
21: Vec X0, X, Xdot;
22: void *funcCtx;
23: TSDGType discgrad; /* Type of electrostatic model */
24: PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *);
25: PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *);
26: PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *);
27: } TS_DiscGrad;
29: static PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
30: {
31: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
33: PetscFunctionBegin;
34: if (X0) {
35: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0));
36: else *X0 = ts->vec_sol;
37: }
38: if (Xdot) {
39: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot));
40: else *Xdot = dg->Xdot;
41: }
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
46: {
47: PetscFunctionBegin;
48: if (X0) {
49: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0));
50: }
51: if (Xdot) {
52: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot));
53: }
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, void *ctx)
58: {
59: PetscFunctionBegin;
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
64: {
65: TS ts = (TS)ctx;
66: Vec X0, Xdot, X0_c, Xdot_c;
68: PetscFunctionBegin;
69: PetscCall(TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot));
70: PetscCall(TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c));
71: PetscCall(MatRestrict(restrct, X0, X0_c));
72: PetscCall(MatRestrict(restrct, Xdot, Xdot_c));
73: PetscCall(VecPointwiseMult(X0_c, rscale, X0_c));
74: PetscCall(VecPointwiseMult(Xdot_c, rscale, Xdot_c));
75: PetscCall(TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot));
76: PetscCall(TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, void *ctx)
81: {
82: PetscFunctionBegin;
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode DMSubDomainRestrictHook_TSDiscGrad(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
87: {
88: TS ts = (TS)ctx;
89: Vec X0, Xdot, X0_sub, Xdot_sub;
91: PetscFunctionBegin;
92: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot));
93: PetscCall(TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub));
95: PetscCall(VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD));
96: PetscCall(VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD));
98: PetscCall(VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD));
99: PetscCall(VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD));
101: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot));
102: PetscCall(TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode TSSetUp_DiscGrad(TS ts)
107: {
108: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
109: DM dm;
111: PetscFunctionBegin;
112: if (!dg->X) PetscCall(VecDuplicate(ts->vec_sol, &dg->X));
113: if (!dg->X0) PetscCall(VecDuplicate(ts->vec_sol, &dg->X0));
114: if (!dg->Xdot) PetscCall(VecDuplicate(ts->vec_sol, &dg->Xdot));
116: PetscCall(TSGetDM(ts, &dm));
117: PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts));
118: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode TSSetFromOptions_DiscGrad(TS ts, PetscOptionItems PetscOptionsObject)
123: {
124: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
126: PetscFunctionBegin;
127: PetscOptionsHeadBegin(PetscOptionsObject, "Discrete Gradients ODE solver options");
128: {
129: PetscCall(PetscOptionsEnum("-ts_discgrad_type", "Type of discrete gradient solver", "TSDiscGradSetDGType", DGTypes, (PetscEnum)dg->discgrad, (PetscEnum *)&dg->discgrad, NULL));
130: }
131: PetscOptionsHeadEnd();
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: static PetscErrorCode TSView_DiscGrad(TS ts, PetscViewer viewer)
136: {
137: PetscBool isascii;
139: PetscFunctionBegin;
140: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
141: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Discrete Gradients\n"));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode TSDiscGradGetType_DiscGrad(TS ts, TSDGType *dgtype)
146: {
147: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
149: PetscFunctionBegin;
150: *dgtype = dg->discgrad;
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode TSDiscGradSetType_DiscGrad(TS ts, TSDGType dgtype)
155: {
156: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
158: PetscFunctionBegin;
159: dg->discgrad = dgtype;
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode TSReset_DiscGrad(TS ts)
164: {
165: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
167: PetscFunctionBegin;
168: PetscCall(VecDestroy(&dg->X));
169: PetscCall(VecDestroy(&dg->X0));
170: PetscCall(VecDestroy(&dg->Xdot));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode TSDestroy_DiscGrad(TS ts)
175: {
176: DM dm;
178: PetscFunctionBegin;
179: PetscCall(TSReset_DiscGrad(ts));
180: PetscCall(TSGetDM(ts, &dm));
181: if (dm) {
182: PetscCall(DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts));
183: PetscCall(DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts));
184: }
185: PetscCall(PetscFree(ts->data));
186: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", NULL));
187: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", NULL));
188: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetType_C", NULL));
189: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetType_C", NULL));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X)
194: {
195: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
196: PetscReal dt = t - ts->ptime;
198: PetscFunctionBegin;
199: PetscCall(VecCopy(ts->vec_sol, dg->X));
200: PetscCall(VecWAXPY(X, dt, dg->Xdot, dg->X));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x)
205: {
206: SNES snes;
207: PetscInt nits, lits;
209: PetscFunctionBegin;
210: PetscCall(TSGetSNES(ts, &snes));
211: PetscCall(SNESSolve(snes, b, x));
212: PetscCall(SNESGetIterationNumber(snes, &nits));
213: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
214: ts->snes_its += nits;
215: ts->ksp_its += lits;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: static PetscErrorCode TSStep_DiscGrad(TS ts)
220: {
221: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
222: TSAdapt adapt;
223: TSStepStatus status = TS_STEP_INCOMPLETE;
224: PetscInt rejections = 0;
225: PetscBool stageok, accept = PETSC_TRUE;
226: PetscReal next_time_step = ts->time_step;
228: PetscFunctionBegin;
229: PetscCall(TSGetAdapt(ts, &adapt));
230: if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, dg->X0));
232: while (!ts->reason && status != TS_STEP_COMPLETE) {
233: PetscReal shift = 1 / (0.5 * ts->time_step);
235: dg->stage_time = ts->ptime + 0.5 * ts->time_step;
237: PetscCall(VecCopy(dg->X0, dg->X));
238: PetscCall(TSPreStage(ts, dg->stage_time));
239: PetscCall(TSDiscGrad_SNESSolve(ts, NULL, dg->X));
240: PetscCall(TSPostStage(ts, dg->stage_time, 0, &dg->X));
241: PetscCall(TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok));
242: if (!stageok) goto reject_step;
244: status = TS_STEP_PENDING;
245: PetscCall(VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X));
246: PetscCall(VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot));
247: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
248: status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
249: if (!accept) {
250: PetscCall(VecCopy(dg->X0, ts->vec_sol));
251: ts->time_step = next_time_step;
252: goto reject_step;
253: }
254: ts->ptime += ts->time_step;
255: ts->time_step = next_time_step;
256: break;
258: reject_step:
259: ts->reject++;
260: accept = PETSC_FALSE;
261: if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) {
262: ts->reason = TS_DIVERGED_STEP_REJECTED;
263: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
264: }
265: }
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y)
270: {
271: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
273: PetscFunctionBegin;
274: if (ns) *ns = 1;
275: if (Y) *Y = &dg->X;
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*
280: This defines the nonlinear equation that is to be solved with SNES
281: G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0
282: */
284: /* x = (x+x')/2 */
285: /* NEED TO CALCULATE x_{n+1} from x and x_{n}*/
286: static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts)
287: {
288: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
289: PetscReal norm, shift = 1 / (0.5 * ts->time_step);
290: PetscInt n, dim;
291: Vec X0, Xdot, Xp, Xdiff;
292: Mat S;
293: PetscScalar F = 0, F0 = 0, Gp;
294: Vec G, SgF;
295: DM dm, dmsave;
297: PetscFunctionBegin;
298: PetscCall(SNESGetDM(snes, &dm));
299: PetscCall(DMGetDimension(dm, &dim));
301: PetscCall(VecDuplicate(y, &Xp));
302: PetscCall(VecDuplicate(y, &Xdiff));
303: PetscCall(VecDuplicate(y, &SgF));
304: PetscCall(VecDuplicate(y, &G));
306: PetscCall(PetscObjectSetName((PetscObject)x, "x"));
307: PetscCall(VecViewFromOptions(x, NULL, "-x_view"));
309: PetscCall(VecGetLocalSize(y, &n));
310: PetscCall(MatCreate(PETSC_COMM_WORLD, &S));
311: PetscCall(MatSetSizes(S, n, n, PETSC_DECIDE, PETSC_DECIDE));
312: PetscCall(MatSetFromOptions(S));
313: PetscInt *S_prealloc_arr;
314: PetscCall(PetscMalloc1(n, &S_prealloc_arr));
315: for (PetscInt i = 0; i < n; ++i) S_prealloc_arr[i] = 2;
316: PetscCall(MatXAIJSetPreallocation(S, 1, S_prealloc_arr, NULL, NULL, NULL));
317: PetscCall(MatSetUp(S));
318: PetscCall((*dg->Sfunc)(ts, dg->stage_time, x, S, dg->funcCtx));
319: PetscCall(PetscFree(S_prealloc_arr));
320: PetscCall(PetscObjectSetName((PetscObject)S, "S"));
321: PetscCall(MatViewFromOptions(S, NULL, "-S_view"));
322: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot));
323: PetscCall(VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x)); /* Xdot = shift (x - X0) */
325: PetscCall(VecAXPBYPCZ(Xp, -1, 2, 0, X0, x)); /* Xp = 2*x - X0 + (0)*Xp */
326: PetscCall(VecAXPBYPCZ(Xdiff, -1, 1, 0, X0, Xp)); /* Xdiff = xp - X0 + (0)*Xdiff */
328: PetscCall(PetscObjectSetName((PetscObject)X0, "X0"));
329: PetscCall(PetscObjectSetName((PetscObject)Xp, "Xp"));
330: PetscCall(VecViewFromOptions(X0, NULL, "-X0_view"));
331: PetscCall(VecViewFromOptions(Xp, NULL, "-Xp_view"));
333: if (dg->discgrad == TS_DG_AVERAGE) {
334: /* Average Value DG:
335: \overline{\nabla} F (x_{n+1},x_{n}) = \int_0^1 \nabla F ((1-\xi)*x_{n+1} + \xi*x_{n}) d \xi */
336: PetscQuadrature quad;
337: PetscInt Nq;
338: const PetscReal *wq, *xq;
339: Vec Xquad, den;
341: PetscCall(PetscObjectSetName((PetscObject)G, "G"));
342: PetscCall(VecDuplicate(G, &Xquad));
343: PetscCall(VecDuplicate(G, &den));
344: PetscCall(VecZeroEntries(G));
346: /* \overline{\nabla} F = \nabla F ((1-\xi) x_{n} + \xi x_{n+1})*/
347: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 2, 0.0, 1.0, &quad));
348: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq));
349: for (PetscInt q = 0; q < Nq; ++q) {
350: PetscReal xi = xq[q], xim1 = 1 - xq[q];
351: PetscCall(VecZeroEntries(Xquad));
352: PetscCall(VecAXPBYPCZ(Xquad, xi, xim1, 1.0, X0, Xp));
353: PetscCall((*dg->Gfunc)(ts, dg->stage_time, Xquad, den, dg->funcCtx));
354: PetscCall(VecAXPY(G, wq[q], den));
355: PetscCall(PetscObjectSetName((PetscObject)den, "den"));
356: PetscCall(VecViewFromOptions(den, NULL, "-den_view"));
357: }
358: PetscCall(VecDestroy(&Xquad));
359: PetscCall(VecDestroy(&den));
360: PetscCall(PetscQuadratureDestroy(&quad));
361: } else if (dg->discgrad == TS_DG_GONZALEZ) {
362: PetscCall((*dg->Ffunc)(ts, dg->stage_time, Xp, &F, dg->funcCtx));
363: PetscCall((*dg->Ffunc)(ts, dg->stage_time, X0, &F0, dg->funcCtx));
364: PetscCall((*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx));
366: /* Adding Extra Gonzalez Term */
367: PetscCall(VecDot(Xdiff, G, &Gp));
368: PetscCall(VecNorm(Xdiff, NORM_2, &norm));
369: if (norm < PETSC_SQRT_MACHINE_EPSILON) {
370: Gp = 0;
371: } else {
372: /* Gp = (1/|xn+1 - xn|^2) * (F(xn+1) - F(xn) - Gp) */
373: Gp = (F - F0 - Gp) / PetscSqr(norm);
374: }
375: PetscCall(VecAXPY(G, Gp, Xdiff));
376: } else if (dg->discgrad == TS_DG_NONE) {
377: PetscCall((*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx));
378: } else {
379: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DG type not supported.");
380: }
381: PetscCall(MatMult(S, G, SgF)); /* Xdot = S*gradF */
383: PetscCall(PetscObjectSetName((PetscObject)G, "G"));
384: PetscCall(VecViewFromOptions(G, NULL, "-G_view"));
385: PetscCall(PetscObjectSetName((PetscObject)SgF, "SgF"));
386: PetscCall(VecViewFromOptions(SgF, NULL, "-SgF_view"));
387: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
388: dmsave = ts->dm;
389: ts->dm = dm;
390: PetscCall(VecAXPBYPCZ(y, 1, -1, 0, Xdot, SgF));
392: ts->dm = dmsave;
393: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot));
395: PetscCall(VecDestroy(&Xp));
396: PetscCall(VecDestroy(&Xdiff));
397: PetscCall(VecDestroy(&SgF));
398: PetscCall(VecDestroy(&G));
399: PetscCall(MatDestroy(&S));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
404: {
405: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
406: PetscReal shift = 1 / (0.5 * ts->time_step);
407: Vec Xdot;
408: DM dm, dmsave;
410: PetscFunctionBegin;
411: PetscCall(SNESGetDM(snes, &dm));
412: /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
413: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot));
415: dmsave = ts->dm;
416: ts->dm = dm;
417: PetscCall(TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE));
418: ts->dm = dmsave;
419: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: 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)
424: {
425: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
427: PetscFunctionBegin;
428: *Sfunc = dg->Sfunc;
429: *Ffunc = dg->Ffunc;
430: *Gfunc = dg->Gfunc;
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: 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)
435: {
436: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
438: PetscFunctionBegin;
439: dg->Sfunc = Sfunc;
440: dg->Ffunc = Ffunc;
441: dg->Gfunc = Gfunc;
442: dg->funcCtx = ctx;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*MC
447: TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method
449: Level: intermediate
451: Notes:
452: This is the implicit midpoint rule, with an optional term that guarantees the discrete
453: gradient property. This timestepper applies to systems of the form $u_t = S(u) \nabla F(u)$
454: where $S(u)$ is a linear operator, and $F$ is a functional of $u$.
456: .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`
457: M*/
458: PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
459: {
460: TS_DiscGrad *th;
462: PetscFunctionBegin;
463: PetscCall(PetscCitationsRegister(DGCitation, &DGCite));
464: ts->ops->reset = TSReset_DiscGrad;
465: ts->ops->destroy = TSDestroy_DiscGrad;
466: ts->ops->view = TSView_DiscGrad;
467: ts->ops->setfromoptions = TSSetFromOptions_DiscGrad;
468: ts->ops->setup = TSSetUp_DiscGrad;
469: ts->ops->step = TSStep_DiscGrad;
470: ts->ops->interpolate = TSInterpolate_DiscGrad;
471: ts->ops->getstages = TSGetStages_DiscGrad;
472: ts->ops->snesfunction = SNESTSFormFunction_DiscGrad;
473: ts->ops->snesjacobian = SNESTSFormJacobian_DiscGrad;
474: ts->default_adapt_type = TSADAPTNONE;
476: ts->usessnes = PETSC_TRUE;
478: PetscCall(PetscNew(&th));
479: ts->data = (void *)th;
481: th->discgrad = TS_DG_NONE;
483: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", TSDiscGradGetFormulation_DiscGrad));
484: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", TSDiscGradSetFormulation_DiscGrad));
485: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetType_C", TSDiscGradGetType_DiscGrad));
486: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetType_C", TSDiscGradSetType_DiscGrad));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@C
491: TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the
492: formulation $u_t = S \nabla F$ for `TSDISCGRAD`
494: Not Collective
496: Input Parameter:
497: . ts - timestepping context
499: Output Parameters:
500: + Sfunc - constructor for the S matrix from the formulation
501: . Ffunc - functional F from the formulation
502: . Gfunc - constructor for the gradient of F from the formulation
503: - ctx - the user context
505: Calling sequence of `Sfunc`:
506: + ts - the integrator
507: . time - the current time
508: . u - the solution
509: . S - the S-matrix from the formulation
510: - ctx - the user context
512: Calling sequence of `Ffunc`:
513: + ts - the integrator
514: . time - the current time
515: . u - the solution
516: . F - the computed function from the formulation
517: - ctx - the user context
519: Calling sequence of `Gfunc`:
520: + ts - the integrator
521: . time - the current time
522: . u - the solution
523: . G - the gradient of the computed function from the formulation
524: - ctx - the user context
526: Level: intermediate
528: .seealso: [](ch_ts), `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`
529: @*/
530: 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)
531: {
532: PetscFunctionBegin;
534: PetscAssertPointer(Sfunc, 2);
535: PetscAssertPointer(Ffunc, 3);
536: PetscAssertPointer(Gfunc, 4);
537: 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));
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@C
542: TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the
543: formulation $u_t = S(u) \nabla F(u)$ for `TSDISCGRAD`
545: Not Collective
547: Input Parameters:
548: + ts - timestepping context
549: . Sfunc - constructor for the S matrix from the formulation
550: . Ffunc - functional F from the formulation
551: . Gfunc - constructor for the gradient of F from the formulation
552: - ctx - optional context for the functions
554: Calling sequence of `Sfunc`:
555: + ts - the integrator
556: . time - the current time
557: . u - the solution
558: . S - the S-matrix from the formulation
559: - ctx - the user context
561: Calling sequence of `Ffunc`:
562: + ts - the integrator
563: . time - the current time
564: . u - the solution
565: . F - the computed function from the formulation
566: - ctx - the user context
568: Calling sequence of `Gfunc`:
569: + ts - the integrator
570: . time - the current time
571: . u - the solution
572: . G - the gradient of the computed function from the formulation
573: - ctx - the user context
575: Level: intermediate
577: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradGetFormulation()`
578: @*/
579: 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)
580: {
581: PetscFunctionBegin;
586: 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));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /*@
591: TSDiscGradGetType - Checks for which discrete gradient to use in formulation for `TSDISCGRAD`
593: Not Collective
595: Input Parameter:
596: . ts - timestepping context
598: Output Parameter:
599: . dgtype - Discrete gradient type <none, gonzalez, average>
601: Level: advanced
603: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradSetType()`
604: @*/
605: PetscErrorCode TSDiscGradGetType(TS ts, TSDGType *dgtype)
606: {
607: PetscFunctionBegin;
609: PetscAssertPointer(dgtype, 2);
610: PetscUseMethod(ts, "TSDiscGradGetType_C", (TS, TSDGType *), (ts, dgtype));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /*@
615: TSDiscGradSetType - Sets discrete gradient formulation.
617: Not Collective
619: Input Parameters:
620: + ts - timestepping context
621: - dgtype - Discrete gradient type <none, gonzalez, average>
623: Options Database Key:
624: . -ts_discgrad_type <type> - flag to choose discrete gradient type
626: Level: intermediate
628: Notes:
629: Without `dgtype` or with type `none`, the discrete gradients timestepper is just implicit midpoint.
631: .seealso: [](ch_ts), `TSDISCGRAD`
632: @*/
633: PetscErrorCode TSDiscGradSetType(TS ts, TSDGType dgtype)
634: {
635: PetscFunctionBegin;
637: PetscTryMethod(ts, "TSDiscGradSetType_C", (TS, TSDGType), (ts, dgtype));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }