Actual source code: basicsymplectic.c
1: /*
2: Code for Timestepping with basic symplectic integrators for separable Hamiltonian systems
3: */
4: #include <petsc/private/tsimpl.h>
5: #include <petscdm.h>
7: static TSBasicSymplecticType TSBasicSymplecticDefault = TSBASICSYMPLECTICSIEULER;
8: static PetscBool TSBasicSymplecticRegisterAllCalled;
9: static PetscBool TSBasicSymplecticPackageInitialized;
11: typedef struct _BasicSymplecticScheme *BasicSymplecticScheme;
12: typedef struct _BasicSymplecticSchemeLink *BasicSymplecticSchemeLink;
14: struct _BasicSymplecticScheme {
15: char *name;
16: PetscInt order;
17: PetscInt s; /* number of stages */
18: PetscReal *c, *d;
19: };
20: struct _BasicSymplecticSchemeLink {
21: struct _BasicSymplecticScheme sch;
22: BasicSymplecticSchemeLink next;
23: };
24: static BasicSymplecticSchemeLink BasicSymplecticSchemeList;
25: typedef struct {
26: TS subts_p, subts_q; /* sub TS contexts that holds the RHSFunction pointers */
27: IS is_p, is_q; /* IS sets for position and momentum respectively */
28: Vec update; /* a nest work vector for generalized coordinates */
29: BasicSymplecticScheme scheme;
30: } TS_BasicSymplectic;
32: /*MC
33: TSBASICSYMPLECTICSIEULER - first order semi-implicit Euler method
35: Level: intermediate
37: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
38: M*/
40: /*MC
41: TSBASICSYMPLECTICVELVERLET - second order Velocity Verlet method (leapfrog method with starting process and determining velocity and position at the same time)
43: Level: intermediate
45: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
46: M*/
48: /*@C
49: TSBasicSymplecticRegisterAll - Registers all of the basic symplectic integration methods in `TSBASICSYMPLECTIC`
51: Not Collective, but should be called by all processes which will need the schemes to be registered
53: Level: advanced
55: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticRegisterDestroy()`
56: @*/
57: PetscErrorCode TSBasicSymplecticRegisterAll(void)
58: {
59: PetscFunctionBegin;
60: if (TSBasicSymplecticRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
61: TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
62: {
63: PetscReal c[1] = {1.0}, d[1] = {1.0};
64: PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER, 1, 1, c, d));
65: }
66: {
67: PetscReal c[2] = {0, 1.0}, d[2] = {0.5, 0.5};
68: PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET, 2, 2, c, d));
69: }
70: {
71: PetscReal c[3] = {1, -2.0 / 3.0, 2.0 / 3.0}, d[3] = {-1.0 / 24.0, 3.0 / 4.0, 7.0 / 24.0};
72: PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTIC3, 3, 3, c, d));
73: }
74: {
75: #define CUBEROOTOFTWO 1.2599210498948731647672106
76: PetscReal c[4] = {1.0 / 2.0 / (2.0 - CUBEROOTOFTWO), (1.0 - CUBEROOTOFTWO) / 2.0 / (2.0 - CUBEROOTOFTWO), (1.0 - CUBEROOTOFTWO) / 2.0 / (2.0 - CUBEROOTOFTWO), 1.0 / 2.0 / (2.0 - CUBEROOTOFTWO)}, d[4] = {1.0 / (2.0 - CUBEROOTOFTWO), -CUBEROOTOFTWO / (2.0 - CUBEROOTOFTWO), 1.0 / (2.0 - CUBEROOTOFTWO), 0};
77: PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTIC4, 4, 4, c, d));
78: }
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*@C
83: TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by `TSBasicSymplecticRegister()`.
85: Not Collective
87: Level: advanced
89: .seealso: [](ch_ts), `TSBasicSymplecticRegister()`, `TSBasicSymplecticRegisterAll()`, `TSBASICSYMPLECTIC`
90: @*/
91: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
92: {
93: BasicSymplecticSchemeLink link;
95: PetscFunctionBegin;
96: while ((link = BasicSymplecticSchemeList)) {
97: BasicSymplecticScheme scheme = &link->sch;
98: BasicSymplecticSchemeList = link->next;
99: PetscCall(PetscFree2(scheme->c, scheme->d));
100: PetscCall(PetscFree(scheme->name));
101: PetscCall(PetscFree(link));
102: }
103: TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /*@C
108: TSBasicSymplecticInitializePackage - This function initializes everything in the `TSBASICSYMPLECTIC` package. It is called
109: from `TSInitializePackage()`.
111: Level: developer
113: .seealso: [](ch_ts), `PetscInitialize()`, `TSBASICSYMPLECTIC`
114: @*/
115: PetscErrorCode TSBasicSymplecticInitializePackage(void)
116: {
117: PetscFunctionBegin;
118: if (TSBasicSymplecticPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
119: TSBasicSymplecticPackageInitialized = PETSC_TRUE;
120: PetscCall(TSBasicSymplecticRegisterAll());
121: PetscCall(PetscRegisterFinalize(TSBasicSymplecticFinalizePackage));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: /*@C
126: TSBasicSymplecticFinalizePackage - This function destroys everything in the `TSBASICSYMPLECTIC` package. It is
127: called from `PetscFinalize()`.
129: Level: developer
131: .seealso: [](ch_ts), `PetscFinalize()`, `TSBASICSYMPLECTIC`
132: @*/
133: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
134: {
135: PetscFunctionBegin;
136: TSBasicSymplecticPackageInitialized = PETSC_FALSE;
137: PetscCall(TSBasicSymplecticRegisterDestroy());
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@C
142: TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.
144: Not Collective, but the same schemes should be registered on all processes on which they will be used
146: Input Parameters:
147: + name - identifier for method
148: . order - approximation order of method
149: . s - number of stages, this is the dimension of the matrices below
150: . c - coefficients for updating generalized position (dimension s)
151: - d - coefficients for updating generalized momentum (dimension s)
153: Level: advanced
155: Note:
156: Several symplectic methods are provided, this function is only needed to create new methods.
158: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
159: @*/
160: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name, PetscInt order, PetscInt s, PetscReal c[], PetscReal d[])
161: {
162: BasicSymplecticSchemeLink link;
163: BasicSymplecticScheme scheme;
165: PetscFunctionBegin;
166: PetscAssertPointer(name, 1);
167: PetscAssertPointer(c, 4);
168: PetscAssertPointer(d, 5);
170: PetscCall(TSBasicSymplecticInitializePackage());
171: PetscCall(PetscNew(&link));
172: scheme = &link->sch;
173: PetscCall(PetscStrallocpy(name, &scheme->name));
174: scheme->order = order;
175: scheme->s = s;
176: PetscCall(PetscMalloc2(s, &scheme->c, s, &scheme->d));
177: PetscCall(PetscArraycpy(scheme->c, c, s));
178: PetscCall(PetscArraycpy(scheme->d, d, s));
179: link->next = BasicSymplecticSchemeList;
180: BasicSymplecticSchemeList = link;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*
185: The simplified form of the equations are:
187: .vb
188: q_{i+1} = q_i + c_i*g(p_i)*h
189: p_{i+1} = p_i + d_i*f(q_{i+1})*h
190: .ve
192: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.
194: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:
195: .vb
196: - Update the position of the particle by adding to it its velocity multiplied by c_1
197: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by d_1
198: - Update the position of the particle by adding to it its (updated) velocity multiplied by c_2
199: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by d_2
200: .ve
202: */
203: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
204: {
205: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
206: BasicSymplecticScheme scheme = bsymp->scheme;
207: Vec solution = ts->vec_sol, update = bsymp->update, q, p, q_update, p_update;
208: IS is_q = bsymp->is_q, is_p = bsymp->is_p;
209: TS subts_q = bsymp->subts_q, subts_p = bsymp->subts_p;
210: PetscBool stageok = PETSC_TRUE;
211: PetscReal ptime, next_time_step = ts->time_step;
212: PetscInt n;
214: PetscFunctionBegin;
215: PetscCall(TSGetStepNumber(ts, &n));
216: PetscCall(TSSetStepNumber(subts_p, n));
217: PetscCall(TSSetStepNumber(subts_q, n));
218: PetscCall(TSGetTime(ts, &ptime));
219: PetscCall(TSSetTime(subts_p, ptime));
220: PetscCall(TSSetTime(subts_q, ptime));
221: PetscCall(VecGetSubVector(update, is_q, &q_update));
222: PetscCall(VecGetSubVector(update, is_p, &p_update));
223: for (PetscInt iter = 0; iter < scheme->s; iter++) {
224: PetscCall(TSPreStage(ts, ptime));
225: PetscCall(VecGetSubVector(solution, is_q, &q));
226: PetscCall(VecGetSubVector(solution, is_p, &p));
227: /* update position q */
228: if (scheme->c[iter]) {
229: PetscCall(TSComputeRHSFunction(subts_q, ptime, p, q_update));
230: PetscCall(VecAXPY(q, scheme->c[iter] * ts->time_step, q_update));
231: }
232: /* update velocity p */
233: if (scheme->d[iter]) {
234: ptime = ptime + scheme->d[iter] * ts->time_step;
235: PetscCall(TSComputeRHSFunction(subts_p, ptime, q, p_update));
236: PetscCall(VecAXPY(p, scheme->d[iter] * ts->time_step, p_update));
237: }
238: PetscCall(VecRestoreSubVector(solution, is_q, &q));
239: PetscCall(VecRestoreSubVector(solution, is_p, &p));
240: PetscCall(TSPostStage(ts, ptime, 0, &solution));
241: PetscCall(TSAdaptCheckStage(ts->adapt, ts, ptime, solution, &stageok));
242: if (!stageok) goto finally;
243: PetscCall(TSFunctionDomainError(ts, ptime, solution, &stageok));
244: if (!stageok) goto finally;
245: }
247: finally:
248: if (!stageok) ts->reason = TS_DIVERGED_STEP_REJECTED;
249: else ts->ptime += next_time_step;
250: PetscCall(VecRestoreSubVector(update, is_q, &q_update));
251: PetscCall(VecRestoreSubVector(update, is_p, &p_update));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, void *ctx)
256: {
257: PetscFunctionBegin;
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
262: {
263: PetscFunctionBegin;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm, DM subdm, void *ctx)
268: {
269: PetscFunctionBegin;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
274: {
275: PetscFunctionBegin;
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
280: {
281: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
282: DM dm;
284: PetscFunctionBegin;
285: PetscCall(TSRHSSplitGetIS(ts, "position", &bsymp->is_q));
286: PetscCall(TSRHSSplitGetIS(ts, "momentum", &bsymp->is_p));
287: PetscCheck(bsymp->is_q && bsymp->is_p, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up RHSSplits with TSRHSSplitSetIS() using split names position and momentum respectively in order to use -ts_type basicsymplectic");
288: PetscCall(TSRHSSplitGetSubTS(ts, "position", &bsymp->subts_q));
289: PetscCall(TSRHSSplitGetSubTS(ts, "momentum", &bsymp->subts_p));
290: PetscCheck(bsymp->subts_q && bsymp->subts_p, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up the RHSFunctions for position and momentum using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS");
292: PetscCall(VecDuplicate(ts->vec_sol, &bsymp->update));
294: PetscCall(TSGetAdapt(ts, &ts->adapt));
295: PetscCall(TSAdaptCandidatesClear(ts->adapt)); /* make sure to use fixed time stepping */
296: PetscCall(TSGetDM(ts, &dm));
297: if (dm) {
298: PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_BasicSymplectic, DMRestrictHook_BasicSymplectic, ts));
299: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_BasicSymplectic, DMSubDomainRestrictHook_BasicSymplectic, ts));
300: }
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
305: {
306: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
308: PetscFunctionBegin;
309: PetscCall(VecDestroy(&bsymp->update));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
314: {
315: PetscFunctionBegin;
316: PetscCall(TSReset_BasicSymplectic(ts));
317: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", NULL));
318: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", NULL));
319: PetscCall(PetscFree(ts->data));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems *PetscOptionsObject)
324: {
325: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
327: PetscFunctionBegin;
328: PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
329: {
330: BasicSymplecticSchemeLink link;
331: PetscInt count, choice;
332: PetscBool flg;
333: const char **namelist;
335: for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++);
336: PetscCall(PetscMalloc1(count, (char ***)&namelist));
337: for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++) namelist[count] = link->sch.name;
338: PetscCall(PetscOptionsEList("-ts_basicsymplectic_type", "Family of basic symplectic integration method", "TSBasicSymplecticSetType", (const char *const *)namelist, count, bsymp->scheme->name, &choice, &flg));
339: if (flg) PetscCall(TSBasicSymplecticSetType(ts, namelist[choice]));
340: PetscCall(PetscFree(namelist));
341: }
342: PetscOptionsHeadEnd();
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode TSView_BasicSymplectic(TS ts, PetscViewer viewer)
347: {
348: PetscFunctionBegin;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts, PetscReal t, Vec X)
353: {
354: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
355: Vec update = bsymp->update;
356: PetscReal alpha = (ts->ptime - t) / ts->time_step;
358: PetscFunctionBegin;
359: PetscCall(VecWAXPY(X, -ts->time_step, update, ts->vec_sol));
360: PetscCall(VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
365: {
366: PetscFunctionBegin;
367: *yr = 1.0 + xr;
368: *yi = xi;
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: /*@
373: TSBasicSymplecticSetType - Set the type of the basic symplectic method
375: Logically Collective
377: Input Parameters:
378: + ts - timestepping context
379: - bsymptype - type of the symplectic scheme
381: Options Database Key:
382: . -ts_basicsymplectic_type <scheme> - select the scheme
384: Level: intermediate
386: Note:
387: The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum".
388: Each split is associated with an `IS` object and a sub-`TS`
389: that is intended to store the user-provided RHS function.
391: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`
392: @*/
393: PetscErrorCode TSBasicSymplecticSetType(TS ts, TSBasicSymplecticType bsymptype)
394: {
395: PetscFunctionBegin;
397: PetscTryMethod(ts, "TSBasicSymplecticSetType_C", (TS, TSBasicSymplecticType), (ts, bsymptype));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: /*@
402: TSBasicSymplecticGetType - Get the type of the basic symplectic method
404: Logically Collective
406: Input Parameters:
407: + ts - timestepping context
408: - bsymptype - type of the basic symplectic scheme
410: Level: intermediate
412: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
413: @*/
414: PetscErrorCode TSBasicSymplecticGetType(TS ts, TSBasicSymplecticType *bsymptype)
415: {
416: PetscFunctionBegin;
418: PetscUseMethod(ts, "TSBasicSymplecticGetType_C", (TS, TSBasicSymplecticType *), (ts, bsymptype));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts, TSBasicSymplecticType bsymptype)
423: {
424: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
425: BasicSymplecticSchemeLink link;
426: PetscBool match;
428: PetscFunctionBegin;
429: if (bsymp->scheme) {
430: PetscCall(PetscStrcmp(bsymp->scheme->name, bsymptype, &match));
431: if (match) PetscFunctionReturn(PETSC_SUCCESS);
432: }
433: for (link = BasicSymplecticSchemeList; link; link = link->next) {
434: PetscCall(PetscStrcmp(link->sch.name, bsymptype, &match));
435: if (match) {
436: bsymp->scheme = &link->sch;
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
439: }
440: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", bsymptype);
441: }
443: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts, TSBasicSymplecticType *bsymptype)
444: {
445: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
447: PetscFunctionBegin;
448: *bsymptype = bsymp->scheme->name;
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*MC
453: TSBASICSYMPLECTIC - ODE solver using basic symplectic integration schemes <https://en.wikipedia.org/wiki/Symplectic_integrator>
455: These methods are intended for separable Hamiltonian systems
457: $$
458: \begin{align*}
459: qdot = dH(q,p,t)/dp \\
460: pdot = -dH(q,p,t)/dq
461: \end{align*}
462: $$
464: where the Hamiltonian can be split into the sum of kinetic energy and potential energy
466: $$
467: H(q,p,t) = T(p,t) + V(q,t).
468: $$
470: As a result, the system can be generally represented by
472: $$
473: \begin{align*}
474: qdot = f(p,t) = dT(p,t)/dp \\
475: pdot = g(q,t) = -dV(q,t)/dq
476: \end{align*}
477: $$
479: and solved iteratively with
481: $$
482: \begin{align*}
483: q_new = q_old + d_i*h*f(p_old,t_old) \\
484: t_new = t_old + d_i*h \\
485: p_new = p_old + c_i*h*g(p_new,t_new) \\
486: i = 0,1,...,n.
487: \end{align*}
488: $$
490: The solution vector should contain both q and p, which correspond to (generalized) position and momentum respectively. Note that the momentum component
491: could simply be velocity in some representations. The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum".
492: Each split is associated with an `IS` object and a sub-`TS` that is intended to store the user-provided RHS function.
494: Level: beginner
496: .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TSRHSSplitSetIS()`, `TSRHSSplitSetRHSFunction()`, `TSType`
497: M*/
498: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
499: {
500: TS_BasicSymplectic *bsymp;
502: PetscFunctionBegin;
503: PetscCall(TSBasicSymplecticInitializePackage());
504: PetscCall(PetscNew(&bsymp));
505: ts->data = (void *)bsymp;
507: ts->ops->setup = TSSetUp_BasicSymplectic;
508: ts->ops->step = TSStep_BasicSymplectic;
509: ts->ops->reset = TSReset_BasicSymplectic;
510: ts->ops->destroy = TSDestroy_BasicSymplectic;
511: ts->ops->setfromoptions = TSSetFromOptions_BasicSymplectic;
512: ts->ops->view = TSView_BasicSymplectic;
513: ts->ops->interpolate = TSInterpolate_BasicSymplectic;
514: ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;
516: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic));
517: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic));
519: PetscCall(TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }