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: PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *);
21: PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *);
22: PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *);
23: } TS_DiscGrad;
25: static PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
26: {
27: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
31: if (X0) {
32: if (dm && dm != ts->dm) {DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0);}
33: else {*X0 = ts->vec_sol;}
34: }
35: if (Xdot) {
36: if (dm && dm != ts->dm) {DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);}
37: else {*Xdot = dg->Xdot;}
38: }
39: return(0);
40: }
42: static PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
43: {
47: if (X0) {
48: if (dm && dm != ts->dm) {DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0);}
49: }
50: if (Xdot) {
51: if (dm && dm != ts->dm) {DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);}
52: }
53: return(0);
54: }
56: static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, void *ctx)
57: {
59: return(0);
60: }
62: static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
63: {
64: TS ts = (TS) ctx;
65: Vec X0, Xdot, X0_c, Xdot_c;
69: TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot);
70: TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c);
71: MatRestrict(restrct, X0, X0_c);
72: MatRestrict(restrct, Xdot, Xdot_c);
73: VecPointwiseMult(X0_c, rscale, X0_c);
74: VecPointwiseMult(Xdot_c, rscale, Xdot_c);
75: TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot);
76: TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c);
77: return(0);
78: }
80: static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, void *ctx)
81: {
83: return(0);
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;
93: TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);
94: TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);
96: VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);
97: VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);
99: VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);
100: VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);
102: TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);
103: TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);
104: return(0);
105: }
107: static PetscErrorCode TSSetUp_DiscGrad(TS ts)
108: {
109: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
110: DM dm;
114: if (!dg->X) {VecDuplicate(ts->vec_sol, &dg->X);}
115: if (!dg->X0) {VecDuplicate(ts->vec_sol, &dg->X0);}
116: if (!dg->Xdot) {VecDuplicate(ts->vec_sol, &dg->Xdot);}
118: TSGetDM(ts, &dm);
119: DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);
120: DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);
121: return(0);
122: }
124: static PetscErrorCode TSSetFromOptions_DiscGrad(PetscOptionItems *PetscOptionsObject, TS ts)
125: {
129: PetscOptionsHead(PetscOptionsObject, "Discrete Gradients ODE solver options");
130: {
131: //PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSDiscGradSetTheta",th->Theta,&th->Theta,NULL);
132: }
133: PetscOptionsTail();
134: return(0);
135: }
137: static PetscErrorCode TSView_DiscGrad(TS ts,PetscViewer viewer)
138: {
139: PetscBool iascii;
143: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
144: if (iascii) {
145: PetscViewerASCIIPrintf(viewer," Discrete Gradients\n");
146: }
147: return(0);
148: }
150: static PetscErrorCode TSReset_DiscGrad(TS ts)
151: {
152: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
156: VecDestroy(&dg->X);
157: VecDestroy(&dg->X0);
158: VecDestroy(&dg->Xdot);
159: return(0);
160: }
162: static PetscErrorCode TSDestroy_DiscGrad(TS ts)
163: {
164: DM dm;
168: TSReset_DiscGrad(ts);
169: TSGetDM(ts, &dm);
170: if (dm) {
171: DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);
172: DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);
173: }
174: PetscFree(ts->data);
175: PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",NULL);
176: PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",NULL);
177: return(0);
178: }
180: static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X)
181: {
182: TS_DiscGrad *dg = (TS_DiscGrad*)ts->data;
183: PetscReal dt = t - ts->ptime;
187: VecCopy(ts->vec_sol, dg->X);
188: VecWAXPY(X, dt, dg->Xdot, dg->X);
189: return(0);
190: }
192: static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x)
193: {
194: SNES snes;
195: PetscInt nits, lits;
199: TSGetSNES(ts, &snes);
200: SNESSolve(snes, b, x);
201: SNESGetIterationNumber(snes, &nits);
202: SNESGetLinearSolveIterations(snes, &lits);
203: ts->snes_its += nits;
204: ts->ksp_its += lits;
205: return(0);
206: }
208: static PetscErrorCode TSStep_DiscGrad(TS ts)
209: {
210: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
211: TSAdapt adapt;
212: TSStepStatus status = TS_STEP_INCOMPLETE;
213: PetscInt rejections = 0;
214: PetscBool stageok, accept = PETSC_TRUE;
215: PetscReal next_time_step = ts->time_step;
219: TSGetAdapt(ts, &adapt);
220: if (!ts->steprollback) {VecCopy(ts->vec_sol, dg->X0);}
222: while (!ts->reason && status != TS_STEP_COMPLETE) {
223: PetscReal shift = 1/(0.5*ts->time_step);
225: dg->stage_time = ts->ptime + 0.5*ts->time_step;
227: VecCopy(dg->X0, dg->X);
228: TSPreStage(ts, dg->stage_time);
229: TSDiscGrad_SNESSolve(ts, NULL, dg->X);
230: TSPostStage(ts, dg->stage_time, 0, &dg->X);
231: TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok);
232: if (!stageok) goto reject_step;
234: status = TS_STEP_PENDING;
235: VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X);
236: VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot);
237: TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
238: status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
239: if (!accept) {
240: VecCopy(dg->X0, ts->vec_sol);
241: ts->time_step = next_time_step;
242: goto reject_step;
243: }
244: ts->ptime += ts->time_step;
245: ts->time_step = next_time_step;
246: break;
248: reject_step:
249: ts->reject++; accept = PETSC_FALSE;
250: if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) {
251: ts->reason = TS_DIVERGED_STEP_REJECTED;
252: PetscInfo2(ts, "Step=%D, step rejections %D greater than current TS allowed, stopping solve\n", ts->steps, rejections);
253: }
254: }
255: return(0);
256: }
258: static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y)
259: {
260: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
263: if (ns) *ns = 1;
264: if (Y) *Y = &(dg->X);
265: return(0);
266: }
268: /*
269: This defines the nonlinear equation that is to be solved with SNES
270: G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0
271: */
272: static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts)
273: {
274: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
275: PetscReal shift = 1/(0.5*ts->time_step);
276: Vec X0, Xdot;
277: DM dm, dmsave;
281: SNESGetDM(snes, &dm);
282: TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);
283: VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x);
285: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
286: dmsave = ts->dm;
287: ts->dm = dm;
288: TSComputeIFunction(ts, dg->stage_time, x, Xdot, y, PETSC_FALSE);
289: ts->dm = dmsave;
290: TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);
291: return(0);
292: }
294: static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
295: {
296: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
297: PetscReal shift = 1/(0.5*ts->time_step);
298: Vec Xdot;
299: DM dm,dmsave;
303: SNESGetDM(snes, &dm);
304: /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
305: TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot);
307: dmsave = ts->dm;
308: ts->dm = dm;
309: TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE);
310: ts->dm = dmsave;
311: TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot);
312: return(0);
313: }
315: 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 *))
316: {
317: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
320: *Sfunc = dg->Sfunc;
321: *Ffunc = dg->Ffunc;
322: *Gfunc = dg->Gfunc;
323: return(0);
324: }
326: 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 *))
327: {
328: TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
331: dg->Sfunc = Sfunc;
332: dg->Ffunc = Ffunc;
333: dg->Gfunc = Gfunc;
334: return(0);
335: }
337: /*MC
338: TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method
340: Level: beginner
342: Notes: This is the implicit midpoint rule, with an optional term that guarantees the discrete gradient property. This
343: timestepper applies to systems of the form
344: $ u_t = S(u) grad F(u)
345: where S(u) is a linear operator, and F is a functional of u.
347: .seealso: TSCreate(), TSSetType(), TS, TSTHETA, TSDiscGradSetFormulation()
348: M*/
349: PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
350: {
351: TS_DiscGrad *th;
355: PetscCitationsRegister(DGCitation, &DGCite);
356: ts->ops->reset = TSReset_DiscGrad;
357: ts->ops->destroy = TSDestroy_DiscGrad;
358: ts->ops->view = TSView_DiscGrad;
359: ts->ops->setfromoptions = TSSetFromOptions_DiscGrad;
360: ts->ops->setup = TSSetUp_DiscGrad;
361: ts->ops->step = TSStep_DiscGrad;
362: ts->ops->interpolate = TSInterpolate_DiscGrad;
363: ts->ops->getstages = TSGetStages_DiscGrad;
364: ts->ops->snesfunction = SNESTSFormFunction_DiscGrad;
365: ts->ops->snesjacobian = SNESTSFormJacobian_DiscGrad;
366: ts->default_adapt_type = TSADAPTNONE;
368: ts->usessnes = PETSC_TRUE;
370: PetscNewLog(ts,&th);
371: ts->data = (void*)th;
372: PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",TSDiscGradGetFormulation_DiscGrad);
373: PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",TSDiscGradSetFormulation_DiscGrad);
374: return(0);
375: }
377: /*@C
378: TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the formulation u_t = S grad F
380: Not Collective
382: Input Parameter:
383: . ts - timestepping context
385: Output Parameters:
386: + Sfunc - constructor for the S matrix from the formulation
387: . Ffunc - functional F from the formulation
388: - Gfunc - constructor for the gradient of F from the formulation
391: Calling sequence of Sfunc:
392: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *)
394: Calling sequence of Ffunc:
395: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *)
397: Calling sequence of Gfunc:
398: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *)
400: Level: intermediate
402: .seealso: TSDiscGradSetFormulation()
403: @*/
404: PetscErrorCode TSDiscGradGetFormulation(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *))
405: {
413: 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*)),(ts,Sfunc,Ffunc,Gfunc));
414: return(0);
415: }
417: /*@C
418: TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the formulation u_t = S(u) grad F(u)
420: Not Collective
422: Input Parameters:
423: + ts - timestepping context
424: . Sfunc - constructor for the S matrix from the formulation
425: . Ffunc - functional F from the formulation
426: - Gfunc - constructor for the gradient of F from the formulation
428: Calling sequence of Sfunc:
429: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *)
431: Calling sequence of Ffunc:
432: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *)
434: Calling sequence of Gfunc:
435: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *)
437: Level: Intermediate
439: .seealso: TSDiscGradGetFormulation()
440: @*/
441: PetscErrorCode TSDiscGradSetFormulation(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)
442: {
450: 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*)),(ts,Sfunc,Ffunc,Gfunc));
451: return(0);
452: }