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: [](chapter_ts), `TSBASICSYMPLECTIC`
38: M*/
40: /*MC
41: TSBASICSYMPLECTICVELVERLET - second order Velocity Verlet method (leapfrog method with starting process and determing velocity and position at the same time)
43: Level: intermediate
45: .seealso: [](chapter_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: [](chapter_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticRegisterDestroy()`
56: @*/
57: PetscErrorCode TSBasicSymplecticRegisterAll(void)
58: {
59: if (TSBasicSymplecticRegisterAllCalled) return 0;
60: TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
61: {
62: PetscReal c[1] = {1.0}, d[1] = {1.0};
63: TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER, 1, 1, c, d);
64: }
65: {
66: PetscReal c[2] = {0, 1.0}, d[2] = {0.5, 0.5};
67: TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET, 2, 2, c, d);
68: }
69: {
70: 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};
71: TSBasicSymplecticRegister(TSBASICSYMPLECTIC3, 3, 3, c, d);
72: }
73: {
74: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
75: PetscReal c[4] = {1.0 / 2.0 / (2.0 - CUBE../../../../..OFTWO), (1.0 - CUBE../../../../..OFTWO) / 2.0 / (2.0 - CUBE../../../../..OFTWO), (1.0 - CUBE../../../../..OFTWO) / 2.0 / (2.0 - CUBE../../../../..OFTWO), 1.0 / 2.0 / (2.0 - CUBE../../../../..OFTWO)}, d[4] = {1.0 / (2.0 - CUBE../../../../..OFTWO), -CUBE../../../../..OFTWO / (2.0 - CUBE../../../../..OFTWO), 1.0 / (2.0 - CUBE../../../../..OFTWO), 0};
76: TSBasicSymplecticRegister(TSBASICSYMPLECTIC4, 4, 4, c, d);
77: }
78: return 0;
79: }
81: /*@C
82: TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by `TSBasicSymplecticRegister()`.
84: Not Collective
86: Level: advanced
88: .seealso: [](chapter_ts), `TSBasicSymplecticRegister()`, `TSBasicSymplecticRegisterAll()`, `TSBASICSYMPLECTIC`
89: @*/
90: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
91: {
92: BasicSymplecticSchemeLink link;
94: while ((link = BasicSymplecticSchemeList)) {
95: BasicSymplecticScheme scheme = &link->sch;
96: BasicSymplecticSchemeList = link->next;
97: PetscFree2(scheme->c, scheme->d);
98: PetscFree(scheme->name);
99: PetscFree(link);
100: }
101: TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
102: return 0;
103: }
105: /*@C
106: TSBasicSymplecticInitializePackage - This function initializes everything in the `TSBASICSYMPLECTIC` package. It is called
107: from `TSInitializePackage()`.
109: Level: developer
111: .seealso: [](chapter_ts), `PetscInitialize()`, `TSBASICSYMPLECTIC`
112: @*/
113: PetscErrorCode TSBasicSymplecticInitializePackage(void)
114: {
115: if (TSBasicSymplecticPackageInitialized) return 0;
116: TSBasicSymplecticPackageInitialized = PETSC_TRUE;
117: TSBasicSymplecticRegisterAll();
118: PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
119: return 0;
120: }
122: /*@C
123: TSBasicSymplecticFinalizePackage - This function destroys everything in the `TSBASICSYMPLECTIC` package. It is
124: called from `PetscFinalize()`.
126: Level: developer
128: .seealso: [](chapter_ts), `PetscFinalize()`, `TSBASICSYMPLECTIC`
129: @*/
130: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
131: {
132: TSBasicSymplecticPackageInitialized = PETSC_FALSE;
133: TSBasicSymplecticRegisterDestroy();
134: return 0;
135: }
137: /*@C
138: TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.
140: Not Collective, but the same schemes should be registered on all processes on which they will be used
142: Input Parameters:
143: + name - identifier for method
144: . order - approximation order of method
145: . s - number of stages, this is the dimension of the matrices below
146: . c - coefficients for updating generalized position (dimension s)
147: - d - coefficients for updating generalized momentum (dimension s)
149: Level: advanced
151: Notes:
152: Several symplectic methods are provided, this function is only needed to create new methods.
154: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`
155: @*/
156: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name, PetscInt order, PetscInt s, PetscReal c[], PetscReal d[])
157: {
158: BasicSymplecticSchemeLink link;
159: BasicSymplecticScheme scheme;
165: TSBasicSymplecticInitializePackage();
166: PetscNew(&link);
167: scheme = &link->sch;
168: PetscStrallocpy(name, &scheme->name);
169: scheme->order = order;
170: scheme->s = s;
171: PetscMalloc2(s, &scheme->c, s, &scheme->d);
172: PetscArraycpy(scheme->c, c, s);
173: PetscArraycpy(scheme->d, d, s);
174: link->next = BasicSymplecticSchemeList;
175: BasicSymplecticSchemeList = link;
176: return 0;
177: }
179: /*
180: The simplified form of the equations are:
182: $ p_{i+1} = p_i + c_i*g(q_i)*h
183: $ q_{i+1} = q_i + d_i*f(p_{i+1},t_{i+1})*h
185: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.
187: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:
189: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
190: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
191: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
192: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2
194: */
195: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
196: {
197: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
198: BasicSymplecticScheme scheme = bsymp->scheme;
199: Vec solution = ts->vec_sol, update = bsymp->update, q, p, q_update, p_update;
200: IS is_q = bsymp->is_q, is_p = bsymp->is_p;
201: TS subts_q = bsymp->subts_q, subts_p = bsymp->subts_p;
202: PetscBool stageok;
203: PetscReal next_time_step = ts->time_step;
204: PetscInt iter;
206: VecGetSubVector(solution, is_q, &q);
207: VecGetSubVector(solution, is_p, &p);
208: VecGetSubVector(update, is_q, &q_update);
209: VecGetSubVector(update, is_p, &p_update);
211: for (iter = 0; iter < scheme->s; iter++) {
212: TSPreStage(ts, ts->ptime);
213: /* update velocity p */
214: if (scheme->c[iter]) {
215: TSComputeRHSFunction(subts_p, ts->ptime, q, p_update);
216: VecAXPY(p, scheme->c[iter] * ts->time_step, p_update);
217: }
218: /* update position q */
219: if (scheme->d[iter]) {
220: TSComputeRHSFunction(subts_q, ts->ptime, p, q_update);
221: VecAXPY(q, scheme->d[iter] * ts->time_step, q_update);
222: ts->ptime = ts->ptime + scheme->d[iter] * ts->time_step;
223: }
224: TSPostStage(ts, ts->ptime, 0, &solution);
225: TSAdaptCheckStage(ts->adapt, ts, ts->ptime, solution, &stageok);
226: if (!stageok) {
227: ts->reason = TS_DIVERGED_STEP_REJECTED;
228: return 0;
229: }
230: TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok);
231: if (!stageok) {
232: ts->reason = TS_DIVERGED_STEP_REJECTED;
233: return 0;
234: }
235: }
237: ts->time_step = next_time_step;
238: VecRestoreSubVector(solution, is_q, &q);
239: VecRestoreSubVector(solution, is_p, &p);
240: VecRestoreSubVector(update, is_q, &q_update);
241: VecRestoreSubVector(update, is_p, &p_update);
242: return 0;
243: }
245: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, void *ctx)
246: {
247: return 0;
248: }
250: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
251: {
252: return 0;
253: }
255: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm, DM subdm, void *ctx)
256: {
257: return 0;
258: }
260: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
261: {
262: return 0;
263: }
265: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
266: {
267: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
268: DM dm;
270: TSRHSSplitGetIS(ts, "position", &bsymp->is_q);
271: TSRHSSplitGetIS(ts, "momentum", &bsymp->is_p);
273: TSRHSSplitGetSubTS(ts, "position", &bsymp->subts_q);
274: TSRHSSplitGetSubTS(ts, "momentum", &bsymp->subts_p);
277: VecDuplicate(ts->vec_sol, &bsymp->update);
279: TSGetAdapt(ts, &ts->adapt);
280: TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
281: TSGetDM(ts, &dm);
282: if (dm) {
283: DMCoarsenHookAdd(dm, DMCoarsenHook_BasicSymplectic, DMRestrictHook_BasicSymplectic, ts);
284: DMSubDomainHookAdd(dm, DMSubDomainHook_BasicSymplectic, DMSubDomainRestrictHook_BasicSymplectic, ts);
285: }
286: return 0;
287: }
289: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
290: {
291: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
293: VecDestroy(&bsymp->update);
294: return 0;
295: }
297: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
298: {
299: TSReset_BasicSymplectic(ts);
300: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", NULL);
301: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", NULL);
302: PetscFree(ts->data);
303: return 0;
304: }
306: static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems *PetscOptionsObject)
307: {
308: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
310: PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
311: {
312: BasicSymplecticSchemeLink link;
313: PetscInt count, choice;
314: PetscBool flg;
315: const char **namelist;
317: for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++)
318: ;
319: PetscMalloc1(count, (char ***)&namelist);
320: for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++) namelist[count] = link->sch.name;
321: PetscOptionsEList("-ts_basicsymplectic_type", "Family of basic symplectic integration method", "TSBasicSymplecticSetType", (const char *const *)namelist, count, bsymp->scheme->name, &choice, &flg);
322: if (flg) TSBasicSymplecticSetType(ts, namelist[choice]);
323: PetscFree(namelist);
324: }
325: PetscOptionsHeadEnd();
326: return 0;
327: }
329: static PetscErrorCode TSView_BasicSymplectic(TS ts, PetscViewer viewer)
330: {
331: return 0;
332: }
334: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts, PetscReal t, Vec X)
335: {
336: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
337: Vec update = bsymp->update;
338: PetscReal alpha = (ts->ptime - t) / ts->time_step;
340: VecWAXPY(X, -ts->time_step, update, ts->vec_sol);
341: VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol);
342: return 0;
343: }
345: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
346: {
347: *yr = 1.0 + xr;
348: *yi = xi;
349: return 0;
350: }
352: /*@C
353: TSBasicSymplecticSetType - Set the type of the basic symplectic method
355: Logically Collective
357: Input Parameters:
358: + ts - timestepping context
359: - bsymptype - type of the symplectic scheme
361: Options Database Key:
362: . -ts_basicsymplectic_type <scheme> - select the scheme
364: Level: intermediate
366: Note:
367: The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum". Each split is associated with an `IS` object and a sub-`TS`
368: that is intended to store the user-provided RHS function.
370: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
371: @*/
372: PetscErrorCode TSBasicSymplecticSetType(TS ts, TSBasicSymplecticType bsymptype)
373: {
375: PetscTryMethod(ts, "TSBasicSymplecticSetType_C", (TS, TSBasicSymplecticType), (ts, bsymptype));
376: return 0;
377: }
379: /*@C
380: TSBasicSymplecticGetType - Get the type of the basic symplectic method
382: Logically Collective
384: Input Parameters:
385: + ts - timestepping context
386: - bsymptype - type of the basic symplectic scheme
388: Level: intermediate
390: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
391: @*/
392: PetscErrorCode TSBasicSymplecticGetType(TS ts, TSBasicSymplecticType *bsymptype)
393: {
395: PetscUseMethod(ts, "TSBasicSymplecticGetType_C", (TS, TSBasicSymplecticType *), (ts, bsymptype));
396: return 0;
397: }
399: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts, TSBasicSymplecticType bsymptype)
400: {
401: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
402: BasicSymplecticSchemeLink link;
403: PetscBool match;
405: if (bsymp->scheme) {
406: PetscStrcmp(bsymp->scheme->name, bsymptype, &match);
407: if (match) return 0;
408: }
409: for (link = BasicSymplecticSchemeList; link; link = link->next) {
410: PetscStrcmp(link->sch.name, bsymptype, &match);
411: if (match) {
412: bsymp->scheme = &link->sch;
413: return 0;
414: }
415: }
416: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", bsymptype);
417: }
419: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts, TSBasicSymplecticType *bsymptype)
420: {
421: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
423: *bsymptype = bsymp->scheme->name;
424: return 0;
425: }
427: /*MC
428: TSBASICSYMPLECTIC - ODE solver using basic symplectic integration schemes
430: These methods are intended for separable Hamiltonian systems
431: .vb
432: qdot = dH(q,p,t)/dp
433: pdot = -dH(q,p,t)/dq
434: .ve
436: where the Hamiltonian can be split into the sum of kinetic energy and potential energy
437: .vb
438: H(q,p,t) = T(p,t) + V(q,t).
439: .ve
441: As a result, the system can be genearlly represented by
442: .vb
443: qdot = f(p,t) = dT(p,t)/dp
444: pdot = g(q,t) = -dV(q,t)/dq
445: .ve
447: and solved iteratively with
448: .vb
449: q_new = q_old + d_i*h*f(p_old,t_old)
450: t_new = t_old + d_i*h
451: p_new = p_old + c_i*h*g(p_new,t_new)
452: i=0,1,...,n.
453: .ve
455: The solution vector should contain both q and p, which correspond to (generalized) position and momentum respectively. Note that the momentum component
456: could simply be velocity in some representations. The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum".
457: Each split is associated with an `IS` object and a sub-`TS` that is intended to store the user-provided RHS function.
459: Level: beginner
461: Reference:
462: . * - wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)
464: .seealso: [](chapter_ts), `TSCreate()`, `TSSetType()`, `TSRHSSplitSetIS()`, `TSRHSSplitSetRHSFunction()`, `TSType`
465: M*/
466: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
467: {
468: TS_BasicSymplectic *bsymp;
470: TSBasicSymplecticInitializePackage();
471: PetscNew(&bsymp);
472: ts->data = (void *)bsymp;
474: ts->ops->setup = TSSetUp_BasicSymplectic;
475: ts->ops->step = TSStep_BasicSymplectic;
476: ts->ops->reset = TSReset_BasicSymplectic;
477: ts->ops->destroy = TSDestroy_BasicSymplectic;
478: ts->ops->setfromoptions = TSSetFromOptions_BasicSymplectic;
479: ts->ops->view = TSView_BasicSymplectic;
480: ts->ops->interpolate = TSInterpolate_BasicSymplectic;
481: ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;
483: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic);
484: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic);
486: TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault);
487: return 0;
488: }