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
34: Level: intermediate
35: .seealso: `TSBASICSYMPLECTIC`
36: M*/
38: /*MC
39: TSBASICSYMPLECTICVELVERLET - second order Velocity Verlet method (leapfrog method with starting process and determing velocity and position at the same time)
40: Level: intermediate
41: .seealso: `TSBASICSYMPLECTIC`
42: M*/
44: /*@C
45: TSBasicSymplecticRegisterAll - Registers all of the basic symplectic integration methods in TSBasicSymplectic
47: Not Collective, but should be called by all processes which will need the schemes to be registered
49: Level: advanced
51: .seealso: `TSBasicSymplecticRegisterDestroy()`
52: @*/
53: PetscErrorCode TSBasicSymplecticRegisterAll(void)
54: {
55: if (TSBasicSymplecticRegisterAllCalled) return 0;
56: TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
57: {
58: PetscReal c[1] = {1.0}, d[1] = {1.0};
59: TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER, 1, 1, c, d);
60: }
61: {
62: PetscReal c[2] = {0, 1.0}, d[2] = {0.5, 0.5};
63: TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET, 2, 2, c, d);
64: }
65: {
66: 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};
67: TSBasicSymplecticRegister(TSBASICSYMPLECTIC3, 3, 3, c, d);
68: }
69: {
70: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
71: 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};
72: TSBasicSymplecticRegister(TSBASICSYMPLECTIC4, 4, 4, c, d);
73: }
74: return 0;
75: }
77: /*@C
78: TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by TSBasicSymplecticRegister().
80: Not Collective
82: Level: advanced
84: .seealso: `TSBasicSymplecticRegister()`, `TSBasicSymplecticRegisterAll()`
85: @*/
86: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
87: {
88: BasicSymplecticSchemeLink link;
90: while ((link = BasicSymplecticSchemeList)) {
91: BasicSymplecticScheme scheme = &link->sch;
92: BasicSymplecticSchemeList = link->next;
93: PetscFree2(scheme->c, scheme->d);
94: PetscFree(scheme->name);
95: PetscFree(link);
96: }
97: TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
98: return 0;
99: }
101: /*@C
102: TSBasicSymplecticInitializePackage - This function initializes everything in the TSBasicSymplectic package. It is called
103: from TSInitializePackage().
105: Level: developer
107: .seealso: `PetscInitialize()`
108: @*/
109: PetscErrorCode TSBasicSymplecticInitializePackage(void)
110: {
111: if (TSBasicSymplecticPackageInitialized) return 0;
112: TSBasicSymplecticPackageInitialized = PETSC_TRUE;
113: TSBasicSymplecticRegisterAll();
114: PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
115: return 0;
116: }
118: /*@C
119: TSBasicSymplecticFinalizePackage - This function destroys everything in the TSBasicSymplectic package. It is
120: called from PetscFinalize().
122: Level: developer
124: .seealso: `PetscFinalize()`
125: @*/
126: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
127: {
128: TSBasicSymplecticPackageInitialized = PETSC_FALSE;
129: TSBasicSymplecticRegisterDestroy();
130: return 0;
131: }
133: /*@C
134: TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.
136: Not Collective, but the same schemes should be registered on all processes on which they will be used
138: Input Parameters:
139: + name - identifier for method
140: . order - approximation order of method
141: . s - number of stages, this is the dimension of the matrices below
142: . c - coefficients for updating generalized position (dimension s)
143: - d - coefficients for updating generalized momentum (dimension s)
145: Notes:
146: Several symplectic methods are provided, this function is only needed to create new methods.
148: Level: advanced
150: .seealso: `TSBasicSymplectic`
151: @*/
152: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name, PetscInt order, PetscInt s, PetscReal c[], PetscReal d[])
153: {
154: BasicSymplecticSchemeLink link;
155: BasicSymplecticScheme scheme;
161: TSBasicSymplecticInitializePackage();
162: PetscNew(&link);
163: scheme = &link->sch;
164: PetscStrallocpy(name, &scheme->name);
165: scheme->order = order;
166: scheme->s = s;
167: PetscMalloc2(s, &scheme->c, s, &scheme->d);
168: PetscArraycpy(scheme->c, c, s);
169: PetscArraycpy(scheme->d, d, s);
170: link->next = BasicSymplecticSchemeList;
171: BasicSymplecticSchemeList = link;
172: return 0;
173: }
175: /*
176: The simplified form of the equations are:
178: $ p_{i+1} = p_i + c_i*g(q_i)*h
179: $ q_{i+1} = q_i + d_i*f(p_{i+1},t_{i+1})*h
181: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.
183: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:
185: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
186: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
187: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
188: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2
190: */
191: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
192: {
193: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
194: BasicSymplecticScheme scheme = bsymp->scheme;
195: Vec solution = ts->vec_sol, update = bsymp->update, q, p, q_update, p_update;
196: IS is_q = bsymp->is_q, is_p = bsymp->is_p;
197: TS subts_q = bsymp->subts_q, subts_p = bsymp->subts_p;
198: PetscBool stageok;
199: PetscReal next_time_step = ts->time_step;
200: PetscInt iter;
202: VecGetSubVector(solution, is_q, &q);
203: VecGetSubVector(solution, is_p, &p);
204: VecGetSubVector(update, is_q, &q_update);
205: VecGetSubVector(update, is_p, &p_update);
207: for (iter = 0; iter < scheme->s; iter++) {
208: TSPreStage(ts, ts->ptime);
209: /* update velocity p */
210: if (scheme->c[iter]) {
211: TSComputeRHSFunction(subts_p, ts->ptime, q, p_update);
212: VecAXPY(p, scheme->c[iter] * ts->time_step, p_update);
213: }
214: /* update position q */
215: if (scheme->d[iter]) {
216: TSComputeRHSFunction(subts_q, ts->ptime, p, q_update);
217: VecAXPY(q, scheme->d[iter] * ts->time_step, q_update);
218: ts->ptime = ts->ptime + scheme->d[iter] * ts->time_step;
219: }
220: TSPostStage(ts, ts->ptime, 0, &solution);
221: TSAdaptCheckStage(ts->adapt, ts, ts->ptime, solution, &stageok);
222: if (!stageok) {
223: ts->reason = TS_DIVERGED_STEP_REJECTED;
224: return 0;
225: }
226: TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok);
227: if (!stageok) {
228: ts->reason = TS_DIVERGED_STEP_REJECTED;
229: return 0;
230: }
231: }
233: ts->time_step = next_time_step;
234: VecRestoreSubVector(solution, is_q, &q);
235: VecRestoreSubVector(solution, is_p, &p);
236: VecRestoreSubVector(update, is_q, &q_update);
237: VecRestoreSubVector(update, is_p, &p_update);
238: return 0;
239: }
241: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, void *ctx)
242: {
243: return 0;
244: }
246: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
247: {
248: return 0;
249: }
251: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm, DM subdm, void *ctx)
252: {
253: return 0;
254: }
256: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
257: {
258: return 0;
259: }
261: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
262: {
263: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
264: DM dm;
266: TSRHSSplitGetIS(ts, "position", &bsymp->is_q);
267: TSRHSSplitGetIS(ts, "momentum", &bsymp->is_p);
269: TSRHSSplitGetSubTS(ts, "position", &bsymp->subts_q);
270: TSRHSSplitGetSubTS(ts, "momentum", &bsymp->subts_p);
273: VecDuplicate(ts->vec_sol, &bsymp->update);
275: TSGetAdapt(ts, &ts->adapt);
276: TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
277: TSGetDM(ts, &dm);
278: if (dm) {
279: DMCoarsenHookAdd(dm, DMCoarsenHook_BasicSymplectic, DMRestrictHook_BasicSymplectic, ts);
280: DMSubDomainHookAdd(dm, DMSubDomainHook_BasicSymplectic, DMSubDomainRestrictHook_BasicSymplectic, ts);
281: }
282: return 0;
283: }
285: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
286: {
287: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
289: VecDestroy(&bsymp->update);
290: return 0;
291: }
293: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
294: {
295: TSReset_BasicSymplectic(ts);
296: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", NULL);
297: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", NULL);
298: PetscFree(ts->data);
299: return 0;
300: }
302: static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems *PetscOptionsObject)
303: {
304: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
306: PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
307: {
308: BasicSymplecticSchemeLink link;
309: PetscInt count, choice;
310: PetscBool flg;
311: const char **namelist;
313: for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++)
314: ;
315: PetscMalloc1(count, (char ***)&namelist);
316: for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++) namelist[count] = link->sch.name;
317: PetscOptionsEList("-ts_basicsymplectic_type", "Family of basic symplectic integration method", "TSBasicSymplecticSetType", (const char *const *)namelist, count, bsymp->scheme->name, &choice, &flg);
318: if (flg) TSBasicSymplecticSetType(ts, namelist[choice]);
319: PetscFree(namelist);
320: }
321: PetscOptionsHeadEnd();
322: return 0;
323: }
325: static PetscErrorCode TSView_BasicSymplectic(TS ts, PetscViewer viewer)
326: {
327: return 0;
328: }
330: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts, PetscReal t, Vec X)
331: {
332: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
333: Vec update = bsymp->update;
334: PetscReal alpha = (ts->ptime - t) / ts->time_step;
336: VecWAXPY(X, -ts->time_step, update, ts->vec_sol);
337: VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol);
338: return 0;
339: }
341: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
342: {
343: *yr = 1.0 + xr;
344: *yi = xi;
345: return 0;
346: }
348: /*@C
349: TSBasicSymplecticSetType - Set the type of the basic symplectic method
351: Logically Collective on TS
353: Input Parameters:
354: + ts - timestepping context
355: - bsymptype - type of the symplectic scheme
357: Options Database:
358: . -ts_basicsymplectic_type <scheme> - select the scheme
360: Notes:
361: 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
362: that is intended to store the user-provided RHS function.
364: Level: intermediate
365: @*/
366: PetscErrorCode TSBasicSymplecticSetType(TS ts, TSBasicSymplecticType bsymptype)
367: {
369: PetscTryMethod(ts, "TSBasicSymplecticSetType_C", (TS, TSBasicSymplecticType), (ts, bsymptype));
370: return 0;
371: }
373: /*@C
374: TSBasicSymplecticGetType - Get the type of the basic symplectic method
376: Logically Collective on TS
378: Input Parameters:
379: + ts - timestepping context
380: - bsymptype - type of the basic symplectic scheme
382: Level: intermediate
383: @*/
384: PetscErrorCode TSBasicSymplecticGetType(TS ts, TSBasicSymplecticType *bsymptype)
385: {
387: PetscUseMethod(ts, "TSBasicSymplecticGetType_C", (TS, TSBasicSymplecticType *), (ts, bsymptype));
388: return 0;
389: }
391: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts, TSBasicSymplecticType bsymptype)
392: {
393: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
394: BasicSymplecticSchemeLink link;
395: PetscBool match;
397: if (bsymp->scheme) {
398: PetscStrcmp(bsymp->scheme->name, bsymptype, &match);
399: if (match) return 0;
400: }
401: for (link = BasicSymplecticSchemeList; link; link = link->next) {
402: PetscStrcmp(link->sch.name, bsymptype, &match);
403: if (match) {
404: bsymp->scheme = &link->sch;
405: return 0;
406: }
407: }
408: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", bsymptype);
409: }
411: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts, TSBasicSymplecticType *bsymptype)
412: {
413: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
415: *bsymptype = bsymp->scheme->name;
416: return 0;
417: }
419: /*MC
420: TSBasicSymplectic - ODE solver using basic symplectic integration schemes
422: These methods are intened for separable Hamiltonian systems
424: $ qdot = dH(q,p,t)/dp
425: $ pdot = -dH(q,p,t)/dq
427: where the Hamiltonian can be split into the sum of kinetic energy and potential energy
429: $ H(q,p,t) = T(p,t) + V(q,t).
431: As a result, the system can be genearlly represented by
433: $ qdot = f(p,t) = dT(p,t)/dp
434: $ pdot = g(q,t) = -dV(q,t)/dq
436: and solved iteratively with
438: $ q_new = q_old + d_i*h*f(p_old,t_old)
439: $ t_new = t_old + d_i*h
440: $ p_new = p_old + c_i*h*g(p_new,t_new)
441: $ i=0,1,...,n.
443: The solution vector should contain both q and p, which correspond to (generalized) position and momentum respectively. Note that the momentum component could simply be velocity in some representations.
444: 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 that is intended to store the user-provided RHS function.
446: Reference: wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)
448: Level: beginner
450: .seealso: `TSCreate()`, `TSSetType()`, `TSRHSSplitSetIS()`, `TSRHSSplitSetRHSFunction()`
452: M*/
453: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
454: {
455: TS_BasicSymplectic *bsymp;
457: TSBasicSymplecticInitializePackage();
458: PetscNew(&bsymp);
459: ts->data = (void *)bsymp;
461: ts->ops->setup = TSSetUp_BasicSymplectic;
462: ts->ops->step = TSStep_BasicSymplectic;
463: ts->ops->reset = TSReset_BasicSymplectic;
464: ts->ops->destroy = TSDestroy_BasicSymplectic;
465: ts->ops->setfromoptions = TSSetFromOptions_BasicSymplectic;
466: ts->ops->view = TSView_BasicSymplectic;
467: ts->ops->interpolate = TSInterpolate_BasicSymplectic;
468: ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;
470: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic);
471: PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic);
473: TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault);
474: return 0;
475: }