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>
6: #include <petsc/private/snesimpl.h>
8: PetscBool DGCite = PETSC_FALSE;
9: const char DGCitation[] = "@article{Gonzalez1996,\n"
10: " title = {Time integration and discrete Hamiltonian systems},\n"
11: " author = {Oscar Gonzalez},\n"
12: " journal = {Journal of Nonlinear Science},\n"
13: " volume = {6},\n"
14: " pages = {449--467},\n"
15: " doi = {10.1007/978-1-4612-1246-1_10},\n"
16: " year = {1996}\n}\n";
18: const char *DGTypes[] = {"gonzalez", "average", "none", "TSDGType", "DG_", NULL};
20: typedef struct {
21: PetscReal stage_time;
22: Vec X0, X, Xdot;
23: void *funcCtx;
24: TSDGType discgrad; /* Type of electrostatic model */
25: PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *);
26: PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *);
27: PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *);
28: PetscErrorCode (*IGfunc)(TS, PetscReal, Vec, Vec, Vec, void *);
29: PetscErrorCode (*IGjac)(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
30: } TS_DiscGrad;
32: /*@
33: TSDiscGradGetX0AndXdot - Gets the last solution and current time derivative held inside the `TS`
35: Collective
37: Input Parameters:
38: + ts - The `TS`
39: - dm - The `DM`, or `NULL` to use the embedded `DM`
41: Output Parameters:
42: + X0 - The solution from the last timestep
43: - Xdot - The current solution time derivative
45: Level: advanced
47: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradRestoreX0AndXdot()`
48: @*/
49: PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
50: {
51: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
53: PetscFunctionBegin;
54: if (X0) {
55: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0));
56: else *X0 = ts->vec_sol;
57: }
58: if (Xdot) {
59: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot));
60: else *Xdot = dg->Xdot;
61: }
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*@
66: TSDiscGradRestoreX0AndXdot - Restores the last solution and current time derivative held inside the `TS`
68: Collective
70: Input Parameters:
71: + ts - The `TS`
72: . dm - The `DM`, or `NULL` to use the embedded `DM`
73: . X0 - The solution from the last timestep
74: - Xdot - The current solution time derivative
76: Level: advanced
78: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradGetX0AndXdot()`
79: @*/
80: PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
81: {
82: PetscFunctionBegin;
83: if (X0) {
84: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0));
85: }
86: if (Xdot) {
87: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot));
88: }
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode DestroyInnerTS_Private(PetscCtxRt ptr)
93: {
94: PetscFunctionBegin;
95: PetscCall(TSDestroy((TS *)ptr));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, PetscCtx ctx)
100: {
101: DMSNES dmsnesFine, dmsnesCoarse;
103: PetscFunctionBegin;
104: PetscCall(DMGetDMSNES(fine, &dmsnesFine));
105: PetscCall(DMGetDMSNES(coarse, &dmsnesCoarse));
106: if (dmsnesCoarse->ops->computefunction == SNESTSFormFunction) {
107: // The context for this callback is a TS, which we need to recreate for the coarse problem
108: TS tsFine, tsCoarse;
109: MPI_Comm comm;
111: PetscCall(PetscContainerGetPointer(dmsnesFine->functionctxcontainer, &tsFine));
112: PetscCall(TSCreate(PetscObjectComm((PetscObject)tsFine), &tsCoarse));
113: PetscCall(TSSetType(tsCoarse, TSDISCGRAD));
114: PetscCall(TSSetDM(tsCoarse, coarse));
115: PetscCall(TSSetFromOptions(tsCoarse));
116: PetscCall(TSSetUp(tsCoarse));
117: {
118: TS_DiscGrad *dgFine = (TS_DiscGrad *)tsFine->data;
119: TS_DiscGrad *dgCoarse = (TS_DiscGrad *)tsCoarse->data;
121: tsCoarse->steps = tsFine->steps;
122: tsCoarse->ptime = tsFine->ptime;
123: tsCoarse->time_step = tsFine->time_step;
124: tsCoarse->time_step0 = tsFine->time_step0;
125: tsCoarse->ptime_prev = tsFine->ptime_prev;
126: tsCoarse->ptime_prev_rollback = tsFine->ptime_prev_rollback;
127: tsCoarse->solvetime = tsFine->solvetime;
129: dgCoarse->stage_time = dgFine->stage_time;
130: dgCoarse->funcCtx = dgFine->funcCtx;
131: dgCoarse->discgrad = dgFine->discgrad;
132: dgCoarse->Sfunc = dgFine->Sfunc;
133: dgCoarse->Gfunc = dgFine->Gfunc;
134: dgCoarse->Ffunc = dgFine->Ffunc;
135: dgCoarse->IGfunc = dgFine->IGfunc;
136: dgCoarse->IGjac = dgFine->IGjac;
137: }
138: PetscCall(PetscObjectGetComm((PetscObject)dmsnesCoarse->functionctxcontainer, &comm));
139: // Check this destruction for memory problems
140: PetscCall(PetscContainerDestroy(&dmsnesCoarse->functionctxcontainer));
141: PetscCall(PetscContainerCreate(comm, &dmsnesCoarse->functionctxcontainer));
142: PetscCall(PetscContainerSetPointer(dmsnesCoarse->functionctxcontainer, tsCoarse));
143: PetscCall(PetscContainerSetCtxDestroy(dmsnesCoarse->functionctxcontainer, DestroyInnerTS_Private));
144: }
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
149: {
150: TS ts = (TS)ctx;
151: Vec X0, Xdot, X0_c, Xdot_c;
153: PetscFunctionBegin;
154: PetscCall(TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot));
155: PetscCall(TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c));
156: PetscCall(MatRestrict(restrct, X0, X0_c));
157: PetscCall(MatRestrict(restrct, Xdot, Xdot_c));
158: PetscCall(VecPointwiseMult(X0_c, rscale, X0_c));
159: PetscCall(VecPointwiseMult(Xdot_c, rscale, Xdot_c));
160: PetscCall(TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot));
161: PetscCall(TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, PetscCtx ctx)
166: {
167: PetscFunctionBegin;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode DMSubDomainRestrictHook_TSDiscGrad(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
172: {
173: TS ts = (TS)ctx;
174: Vec X0, Xdot, X0_sub, Xdot_sub;
176: PetscFunctionBegin;
177: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot));
178: PetscCall(TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub));
180: PetscCall(VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD));
181: PetscCall(VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD));
183: PetscCall(VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD));
184: PetscCall(VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD));
186: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot));
187: PetscCall(TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode TSSetUp_DiscGrad(TS ts)
192: {
193: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
194: DM dm;
196: PetscFunctionBegin;
197: if (!dg->X) PetscCall(VecDuplicate(ts->vec_sol, &dg->X));
198: if (!dg->X0) PetscCall(VecDuplicate(ts->vec_sol, &dg->X0));
199: if (!dg->Xdot) PetscCall(VecDuplicate(ts->vec_sol, &dg->Xdot));
201: PetscCall(TSGetDM(ts, &dm));
202: PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts));
203: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode TSSetFromOptions_DiscGrad(TS ts, PetscOptionItems PetscOptionsObject)
208: {
209: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
211: PetscFunctionBegin;
212: PetscOptionsHeadBegin(PetscOptionsObject, "Discrete Gradients ODE solver options");
213: {
214: PetscCall(PetscOptionsEnum("-ts_discgrad_type", "Type of discrete gradient solver", "TSDiscGradSetDGType", DGTypes, (PetscEnum)dg->discgrad, (PetscEnum *)&dg->discgrad, NULL));
215: }
216: PetscOptionsHeadEnd();
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: static PetscErrorCode TSView_DiscGrad(TS ts, PetscViewer viewer)
221: {
222: PetscBool isascii;
224: PetscFunctionBegin;
225: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
226: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Discrete Gradients\n"));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode TSDiscGradGetType_DiscGrad(TS ts, TSDGType *dgtype)
231: {
232: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
234: PetscFunctionBegin;
235: *dgtype = dg->discgrad;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: static PetscErrorCode TSDiscGradSetType_DiscGrad(TS ts, TSDGType dgtype)
240: {
241: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
243: PetscFunctionBegin;
244: dg->discgrad = dgtype;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode TSReset_DiscGrad(TS ts)
249: {
250: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
252: PetscFunctionBegin;
253: PetscCall(VecDestroy(&dg->X));
254: PetscCall(VecDestroy(&dg->X0));
255: PetscCall(VecDestroy(&dg->Xdot));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: static PetscErrorCode TSDestroy_DiscGrad(TS ts)
260: {
261: DM dm;
263: PetscFunctionBegin;
264: PetscCall(TSReset_DiscGrad(ts));
265: PetscCall(TSGetDM(ts, &dm));
266: if (dm) {
267: PetscCall(DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts));
268: PetscCall(DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts));
269: }
270: PetscCall(PetscFree(ts->data));
271: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", NULL));
272: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", NULL));
273: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetType_C", NULL));
274: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetType_C", NULL));
275: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetImplicitFormulation_C", NULL));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X)
280: {
281: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
282: PetscReal dt = t - ts->ptime;
284: PetscFunctionBegin;
285: PetscCall(VecCopy(ts->vec_sol, dg->X));
286: PetscCall(VecWAXPY(X, dt, dg->Xdot, dg->X));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x)
291: {
292: SNES snes;
293: PetscInt nits, lits;
295: PetscFunctionBegin;
296: PetscCall(TSGetSNES(ts, &snes));
297: PetscCall(SNESSolve(snes, b, x));
298: PetscCall(SNESGetIterationNumber(snes, &nits));
299: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
300: ts->snes_its += nits;
301: ts->ksp_its += lits;
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode TSStep_DiscGrad(TS ts)
306: {
307: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
308: TSAdapt adapt;
309: TSStepStatus status = TS_STEP_INCOMPLETE;
310: PetscInt rejections = 0;
311: PetscBool stageok, accept = PETSC_TRUE;
312: PetscReal next_time_step = ts->time_step;
314: PetscFunctionBegin;
315: PetscCall(TSGetAdapt(ts, &adapt));
316: if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, dg->X0));
318: while (!ts->reason && status != TS_STEP_COMPLETE) {
319: PetscReal shift = 1 / (0.5 * ts->time_step);
321: dg->stage_time = ts->ptime + 0.5 * ts->time_step;
323: PetscCall(VecCopy(dg->X0, dg->X));
324: PetscCall(TSPreStage(ts, dg->stage_time));
325: PetscCall(TSDiscGrad_SNESSolve(ts, NULL, dg->X));
326: PetscCall(TSPostStage(ts, dg->stage_time, 0, &dg->X));
327: PetscCall(TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok));
328: if (!stageok) goto reject_step;
330: status = TS_STEP_PENDING;
331: PetscCall(VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X));
332: PetscCall(VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot));
333: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
334: status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
335: if (!accept) {
336: PetscCall(VecCopy(dg->X0, ts->vec_sol));
337: ts->time_step = next_time_step;
338: goto reject_step;
339: }
340: ts->ptime += ts->time_step;
341: ts->time_step = next_time_step;
342: break;
344: reject_step:
345: ts->reject++;
346: accept = PETSC_FALSE;
347: if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) {
348: ts->reason = TS_DIVERGED_STEP_REJECTED;
349: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
350: }
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y)
356: {
357: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
359: PetscFunctionBegin;
360: if (ns) *ns = 1;
361: if (Y) *Y = &dg->X;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*
366: This defines the nonlinear equation that is to be solved with SNES
367: G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0
368: */
370: /* x = (x+x')/2 */
371: /* NEED TO CALCULATE x_{n+1} from x and x_{n}*/
372: static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts)
373: {
374: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
375: PetscReal norm, shift = 1 / (0.5 * ts->time_step);
376: PetscInt n, dim;
377: Vec X0, Xdot, Xp, Xdiff;
378: Mat S;
379: PetscInt *S_prealloc_arr;
380: PetscReal Snorm;
381: PetscScalar F = 0, F0 = 0, Gp;
382: Vec G, SgF;
383: DM dm, dmsave;
384: MPI_Comm comm;
386: PetscFunctionBegin;
387: PetscCall(SNESGetDM(snes, &dm));
388: PetscCall(DMGetDimension(dm, &dim));
389: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
391: PetscCall(VecDuplicate(y, &Xp));
392: PetscCall(VecDuplicate(y, &Xdiff));
393: PetscCall(VecDuplicate(y, &SgF));
394: PetscCall(VecDuplicate(y, &G));
396: PetscCall(PetscObjectSetName((PetscObject)x, "x"));
397: PetscCall(VecViewFromOptions(x, NULL, "-x_view"));
399: PetscCall(VecGetLocalSize(y, &n));
400: PetscCall(MatCreate(comm, &S));
401: PetscCall(MatSetSizes(S, n, n, PETSC_DECIDE, PETSC_DECIDE));
402: PetscCall(MatSetFromOptions(S));
403: PetscCall(PetscMalloc1(n, &S_prealloc_arr));
404: for (PetscInt i = 0; i < n; ++i) S_prealloc_arr[i] = 2;
405: PetscCall(MatXAIJSetPreallocation(S, 1, S_prealloc_arr, NULL, NULL, NULL));
406: PetscCall(MatSetUp(S));
407: PetscCheck(dg->Sfunc, comm, PETSC_ERR_ARG_WRONGSTATE, "Sfunc is not set, call TSDiscGradSetFormulation()");
408: PetscCall((*dg->Sfunc)(ts, dg->stage_time, x, S, dg->funcCtx));
409: PetscCall(PetscFree(S_prealloc_arr));
410: PetscCall(PetscObjectSetName((PetscObject)S, "S"));
411: PetscCall(MatViewFromOptions(S, NULL, "-S_view"));
412: PetscCall(MatNorm(S, NORM_FROBENIUS, &Snorm));
413: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot));
414: PetscCall(VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x)); /* Xdot = shift (x - X0) */
416: PetscCall(VecAXPBYPCZ(Xp, -1, 2, 0, X0, x)); /* Xp = 2*x - X0 + (0)*Xp */
417: PetscCall(VecAXPBYPCZ(Xdiff, -1, 1, 0, X0, Xp)); /* Xdiff = xp - X0 + (0)*Xdiff */
419: PetscCall(PetscObjectSetName((PetscObject)X0, "X0"));
420: PetscCall(PetscObjectSetName((PetscObject)Xp, "Xp"));
421: PetscCall(VecViewFromOptions(X0, NULL, "-X0_view"));
422: PetscCall(VecViewFromOptions(Xp, NULL, "-Xp_view"));
424: if (Snorm == 0.) {
425: PetscCall(VecZeroEntries(G));
426: PetscCall(VecViewFromOptions(x, NULL, "-y_view"));
427: PetscCall(VecViewFromOptions(Xdot, NULL, "-Xdot_view"));
428: PetscCheck(dg->IGfunc, comm, PETSC_ERR_ARG_WRONGSTATE, "IGfunc is not set, call TSDiscGradSetImplicitFormulation()");
429: PetscCall((*dg->IGfunc)(ts, dg->stage_time, Xp, Xdot, y, dg->funcCtx));
430: goto end;
431: }
432: if (dg->discgrad == TS_DG_AVERAGE) {
433: /* Average Value DG:
434: \overline{\nabla} F (x_{n+1},x_{n}) = \int_0^1 \nabla F ((1-\xi)*x_{n+1} + \xi*x_{n}) d \xi */
435: PetscQuadrature quad;
436: PetscInt Nq;
437: const PetscReal *wq, *xq;
438: Vec Xquad, den;
440: PetscCheck(dg->Gfunc, comm, PETSC_ERR_ARG_WRONGSTATE, "Gfunc is not set, call TSDiscGradSetFormulation()");
441: PetscCall(PetscObjectSetName((PetscObject)G, "G"));
442: PetscCall(VecDuplicate(G, &Xquad));
443: PetscCall(VecDuplicate(G, &den));
444: PetscCall(VecZeroEntries(G));
446: /* \overline{\nabla} F = \nabla F ((1-\xi) x_{n} + \xi x_{n+1})*/
447: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 2, 0.0, 1.0, &quad));
448: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq));
449: for (PetscInt q = 0; q < Nq; ++q) {
450: PetscReal xi = xq[q], xim1 = 1 - xq[q];
451: PetscCall(VecZeroEntries(Xquad));
452: PetscCall(VecAXPBYPCZ(Xquad, xi, xim1, 1.0, X0, Xp));
453: PetscCall((*dg->Gfunc)(ts, dg->stage_time, Xquad, den, dg->funcCtx));
454: PetscCall(VecAXPY(G, wq[q], den));
455: PetscCall(PetscObjectSetName((PetscObject)den, "den"));
456: PetscCall(VecViewFromOptions(den, NULL, "-den_view"));
457: }
458: PetscCall(VecDestroy(&Xquad));
459: PetscCall(VecDestroy(&den));
460: PetscCall(PetscQuadratureDestroy(&quad));
461: } else if (dg->discgrad == TS_DG_GONZALEZ) {
462: PetscCheck(dg->Ffunc, comm, PETSC_ERR_ARG_WRONGSTATE, "Ffunc is not set, call TSDiscGradSetFormulation()");
463: PetscCheck(dg->Gfunc, comm, PETSC_ERR_ARG_WRONGSTATE, "Gfunc is not set, call TSDiscGradSetFormulation()");
464: PetscCall((*dg->Ffunc)(ts, dg->stage_time, Xp, &F, dg->funcCtx));
465: PetscCall((*dg->Ffunc)(ts, dg->stage_time, X0, &F0, dg->funcCtx));
466: PetscCall((*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx));
468: /* Adding Extra Gonzalez Term */
469: PetscCall(VecDot(Xdiff, G, &Gp));
470: PetscCall(VecNorm(Xdiff, NORM_2, &norm));
471: if (norm < PETSC_SQRT_MACHINE_EPSILON) {
472: Gp = 0;
473: } else {
474: /* Gp = (1/|xn+1 - xn|^2) * (F(xn+1) - F(xn) - Gp) */
475: Gp = (F - F0 - Gp) / PetscSqr(norm);
476: }
477: PetscCall(VecAXPY(G, Gp, Xdiff));
478: } else if (dg->discgrad == TS_DG_NONE) {
479: PetscCheck(dg->Gfunc, comm, PETSC_ERR_ARG_WRONGSTATE, "Gfunc is not set, call TSDiscGradSetFormulation()");
480: PetscCall((*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx));
481: } else {
482: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DG type not supported.");
483: }
484: PetscCall(MatMult(S, G, SgF)); /* Xdot = S*gradF */
486: PetscCall(PetscObjectSetName((PetscObject)G, "G"));
487: PetscCall(VecViewFromOptions(G, NULL, "-G_view"));
488: PetscCall(PetscObjectSetName((PetscObject)SgF, "SgF"));
489: PetscCall(VecViewFromOptions(SgF, NULL, "-SgF_view"));
490: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
491: dmsave = ts->dm;
492: ts->dm = dm;
493: PetscCall(VecAXPBYPCZ(y, 1, -1, 0, Xdot, SgF));
495: ts->dm = dmsave;
496: end:
497: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot));
499: PetscCall(VecDestroy(&Xp));
500: PetscCall(VecDestroy(&Xdiff));
501: PetscCall(VecDestroy(&SgF));
502: PetscCall(VecDestroy(&G));
503: PetscCall(MatDestroy(&S));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
508: {
509: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
510: PetscReal shift = 1 / (0.5 * ts->time_step);
511: Vec Xdot;
512: DM dm, dmsave;
514: PetscFunctionBegin;
515: PetscCall(SNESGetDM(snes, &dm));
516: /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
517: PetscCall(TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot));
519: dmsave = ts->dm;
520: ts->dm = dm;
521: if (dg->IGjac) PetscCall(TSComputeIJacobian_Internal(ts, dg->IGjac, NULL, dg->funcCtx, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE));
522: else PetscCall(TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE));
523: ts->dm = dmsave;
524: PetscCall(TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot));
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: 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 *), PetscCtx ctx)
529: {
530: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
532: PetscFunctionBegin;
533: *Sfunc = dg->Sfunc;
534: *Ffunc = dg->Ffunc;
535: *Gfunc = dg->Gfunc;
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: 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 *), PetscCtx ctx)
540: {
541: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
543: PetscFunctionBegin;
544: dg->Sfunc = Sfunc;
545: dg->Ffunc = Ffunc;
546: dg->Gfunc = Gfunc;
547: dg->funcCtx = ctx;
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: static PetscErrorCode TSDiscGradSetImplicitFormulation_DiscGrad(TS ts, PetscErrorCode (*IGfunc)(TS ts, PetscReal time, Vec u, Vec u_t, Vec G, void *ctx), PetscErrorCode (*IGjac)(TS ts, PetscReal time, Vec u, Vec u_t, PetscReal shift, Mat J, Mat Jp, void *ctx))
552: {
553: TS_DiscGrad *dg = (TS_DiscGrad *)ts->data;
555: PetscFunctionBegin;
556: dg->IGfunc = IGfunc;
557: dg->IGjac = IGjac;
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /*MC
562: TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method
564: Level: intermediate
566: Notes:
567: This is the implicit midpoint rule, with an optional term that guarantees the discrete
568: gradient property. This timestepper applies to systems of the form $u_t = S(u) \nabla F(u)$
569: where $S(u)$ is a linear operator, and $F$ is a functional of $u$.
571: For Hamiltonian systems designed to conserve the first integral (energy),
572: but also has the property for some systems of monotonicity in a functional.
574: .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TS`, `TSType`, `TSDiscGradSetFormulation()`, `TSDiscGradGetFormulation()`,
575: `TSDiscGradSetType()`, `TSDiscGradGetType()`, `TSDGType`
576: M*/
577: PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
578: {
579: TS_DiscGrad *th;
581: PetscFunctionBegin;
582: PetscCall(PetscCitationsRegister(DGCitation, &DGCite));
583: ts->ops->reset = TSReset_DiscGrad;
584: ts->ops->destroy = TSDestroy_DiscGrad;
585: ts->ops->view = TSView_DiscGrad;
586: ts->ops->setfromoptions = TSSetFromOptions_DiscGrad;
587: ts->ops->setup = TSSetUp_DiscGrad;
588: ts->ops->step = TSStep_DiscGrad;
589: ts->ops->interpolate = TSInterpolate_DiscGrad;
590: ts->ops->getstages = TSGetStages_DiscGrad;
591: ts->ops->snesfunction = SNESTSFormFunction_DiscGrad;
592: ts->ops->snesjacobian = SNESTSFormJacobian_DiscGrad;
593: ts->default_adapt_type = TSADAPTNONE;
595: ts->usessnes = PETSC_TRUE;
597: PetscCall(PetscNew(&th));
598: ts->data = (void *)th;
600: th->discgrad = TS_DG_NONE;
602: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetFormulation_C", TSDiscGradGetFormulation_DiscGrad));
603: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetFormulation_C", TSDiscGradSetFormulation_DiscGrad));
604: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradGetType_C", TSDiscGradGetType_DiscGrad));
605: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetType_C", TSDiscGradSetType_DiscGrad));
606: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDiscGradSetImplicitFormulation_C", TSDiscGradSetImplicitFormulation_DiscGrad));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /*@C
611: TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the
612: formulation $u_t = S \nabla F$ for `TSDISCGRAD`
614: Not Collective
616: Input Parameter:
617: . ts - timestepping context
619: Output Parameters:
620: + Sfunc - constructor for the S matrix from the formulation
621: . Ffunc - functional F from the formulation
622: . Gfunc - constructor for the gradient of F from the formulation
623: - ctx - the user context
625: Calling sequence of `Sfunc`:
626: + ts - the integrator
627: . time - the current time
628: . u - the solution
629: . S - the S-matrix from the formulation
630: - ctx - the user context
632: Calling sequence of `Ffunc`:
633: + ts - the integrator
634: . time - the current time
635: . u - the solution
636: . F - the computed function from the formulation
637: - ctx - the user context
639: Calling sequence of `Gfunc`:
640: + ts - the integrator
641: . time - the current time
642: . u - the solution
643: . G - the gradient of the computed function from the formulation
644: - ctx - the user context
646: Level: intermediate
648: .seealso: [](ch_ts), `TS`, `TSDISCGRAD`, `TSDiscGradSetFormulation()`, `TSDiscGradSetImplicitFormulation()`
649: @*/
650: PetscErrorCode TSDiscGradGetFormulation(TS ts, PetscErrorCode (**Sfunc)(TS ts, PetscReal time, Vec u, Mat S, PetscCtx ctx), PetscErrorCode (**Ffunc)(TS ts, PetscReal time, Vec u, PetscScalar *F, PetscCtx ctx), PetscErrorCode (**Gfunc)(TS ts, PetscReal time, Vec u, Vec G, PetscCtx ctx), PetscCtx ctx)
651: {
652: PetscFunctionBegin;
654: PetscAssertPointer(Sfunc, 2);
655: PetscAssertPointer(Ffunc, 3);
656: PetscAssertPointer(Gfunc, 4);
657: 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));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*@C
662: TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the
663: formulation $u_t = S(u) \nabla F(u)$ for `TSDISCGRAD`
665: Not Collective
667: Input Parameters:
668: + ts - timestepping context
669: . Sfunc - constructor for the S matrix from the formulation
670: . Ffunc - functional F from the formulation
671: . Gfunc - constructor for the gradient of F from the formulation
672: - ctx - optional context for the functions
674: Calling sequence of `Sfunc`:
675: + ts - the integrator
676: . time - the current time
677: . u - the solution
678: . S - the S-matrix from the formulation
679: - ctx - the user context
681: Calling sequence of `Ffunc`:
682: + ts - the integrator
683: . time - the current time
684: . u - the solution
685: . F - the computed function from the formulation
686: - ctx - the user context
688: Calling sequence of `Gfunc`:
689: + ts - the integrator
690: . time - the current time
691: . u - the solution
692: . G - the gradient of the computed function from the formulation
693: - ctx - the user context
695: Level: intermediate
697: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradGetFormulation()`, `TSDiscGradSetImplicitFormulation()`
698: @*/
699: PetscErrorCode TSDiscGradSetFormulation(TS ts, PetscErrorCode (*Sfunc)(TS ts, PetscReal time, Vec u, Mat S, PetscCtx ctx), PetscErrorCode (*Ffunc)(TS ts, PetscReal time, Vec u, PetscScalar *F, PetscCtx ctx), PetscErrorCode (*Gfunc)(TS ts, PetscReal time, Vec u, Vec G, PetscCtx ctx), PetscCtx ctx)
700: {
701: PetscFunctionBegin;
706: 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));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*@
711: TSDiscGradGetType - Checks for which discrete gradient to use in formulation for `TSDISCGRAD`
713: Not Collective
715: Input Parameter:
716: . ts - timestepping context
718: Output Parameter:
719: . dgtype - Discrete gradient type <none, gonzalez, average>
721: Level: advanced
723: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradSetType()`
724: @*/
725: PetscErrorCode TSDiscGradGetType(TS ts, TSDGType *dgtype)
726: {
727: PetscFunctionBegin;
729: PetscAssertPointer(dgtype, 2);
730: PetscUseMethod(ts, "TSDiscGradGetType_C", (TS, TSDGType *), (ts, dgtype));
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /*@
735: TSDiscGradSetType - Sets discrete gradient formulation.
737: Not Collective
739: Input Parameters:
740: + ts - timestepping context
741: - dgtype - Discrete gradient type <none, gonzalez, average>
743: Options Database Key:
744: . -ts_discgrad_type type - flag to choose discrete gradient type
746: Level: intermediate
748: Notes:
749: Without `dgtype` or with type `none`, the discrete gradients timestepper is just implicit midpoint.
751: .seealso: [](ch_ts), `TSDISCGRAD`
752: @*/
753: PetscErrorCode TSDiscGradSetType(TS ts, TSDGType dgtype)
754: {
755: PetscFunctionBegin;
757: PetscTryMethod(ts, "TSDiscGradSetType_C", (TS, TSDGType), (ts, dgtype));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*@C
762: TSDiscGradSetImplicitFormulation - Set the construction method for IG and its Jacobian from the
763: formulation $IG(t) = u_t - S(u) \nabla F(u)$ for `TSDISCGRAD`
765: Not Collective
767: Input Parameters:
768: + ts - timestepping context
769: . IGfunc - implicit formulation
770: - IGjac - implicit Jacobian
772: Calling sequence of `IGfunc`:
773: + ts - the integrator
774: . time - the current time
775: . u - the solution
776: . u_t - the time derivative
777: . G - the LHS of the formulation
778: - ctx - the user context
780: Calling sequence of `IGjac`:
781: + ts - the integrator
782: . time - the current time
783: . u - the solution
784: . u_t - the time derivative
785: . shift - the multiplier for the time derivative part
786: . J - the Jacobian of the LHS
787: . Jp - the Jacobian preconditioner of the LHS
788: - ctx - the user context
790: Level: intermediate
792: Note:
793: This allows the DG formulation to be given in the PETSc style for fully implicit solvers.
795: .seealso: [](ch_ts), `TSDISCGRAD`, `TSDiscGradGetFormulation()`, `TSDiscGradSetFormulation()`
796: @*/
797: PetscErrorCode TSDiscGradSetImplicitFormulation(TS ts, PetscErrorCode (*IGfunc)(TS ts, PetscReal time, Vec u, Vec u_t, Vec G, void *ctx), PetscErrorCode (*IGjac)(TS ts, PetscReal time, Vec u, Vec u_t, PetscReal shift, Mat J, Mat Jp, void *ctx))
798: {
799: PetscFunctionBegin;
803: PetscTryMethod(ts, "TSDiscGradSetImplicitFormulation_C", (TS, PetscErrorCode (*)(TS, PetscReal, Vec, Vec, Vec, void *), PetscErrorCode (*)(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *)), (ts, IGfunc, IGjac));
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }