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) {ts->reason = TS_DIVERGED_STEP_REJECTED; return 0;}
223: TSFunctionDomainError(ts,ts->ptime+ts->time_step,update,&stageok);
224: if (!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return 0;}
225: }
227: ts->time_step = next_time_step;
228: VecRestoreSubVector(solution,is_q,&q);
229: VecRestoreSubVector(solution,is_p,&p);
230: VecRestoreSubVector(update,is_q,&q_update);
231: VecRestoreSubVector(update,is_p,&p_update);
232: return 0;
233: }
235: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine,DM coarse,void *ctx)
236: {
237: return 0;
238: }
240: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
241: {
242: return 0;
243: }
245: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm,DM subdm,void *ctx)
246: {
247: return 0;
248: }
250: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
251: {
252: return 0;
253: }
255: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
256: {
257: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
258: DM dm;
260: TSRHSSplitGetIS(ts,"position",&bsymp->is_q);
261: TSRHSSplitGetIS(ts,"momentum",&bsymp->is_p);
263: TSRHSSplitGetSubTS(ts,"position",&bsymp->subts_q);
264: TSRHSSplitGetSubTS(ts,"momentum",&bsymp->subts_p);
267: VecDuplicate(ts->vec_sol,&bsymp->update);
269: TSGetAdapt(ts,&ts->adapt);
270: TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
271: TSGetDM(ts,&dm);
272: if (dm) {
273: DMCoarsenHookAdd(dm,DMCoarsenHook_BasicSymplectic,DMRestrictHook_BasicSymplectic,ts);
274: DMSubDomainHookAdd(dm,DMSubDomainHook_BasicSymplectic,DMSubDomainRestrictHook_BasicSymplectic,ts);
275: }
276: return 0;
277: }
279: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
280: {
281: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
283: VecDestroy(&bsymp->update);
284: return 0;
285: }
287: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
288: {
289: TSReset_BasicSymplectic(ts);
290: PetscFree(ts->data);
291: return 0;
292: }
294: static PetscErrorCode TSSetFromOptions_BasicSymplectic(PetscOptionItems *PetscOptionsObject,TS ts)
295: {
296: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
298: PetscOptionsHead(PetscOptionsObject,"Basic symplectic integrator options");
299: {
300: BasicSymplecticSchemeLink link;
301: PetscInt count,choice;
302: PetscBool flg;
303: const char **namelist;
305: for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) ;
306: PetscMalloc1(count,(char***)&namelist);
307: for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) namelist[count] = link->sch.name;
308: PetscOptionsEList("-ts_basicsymplectic_type","Family of basic symplectic integration method","TSBasicSymplecticSetType",(const char*const*)namelist,count,bsymp->scheme->name,&choice,&flg);
309: if (flg) TSBasicSymplecticSetType(ts,namelist[choice]);
310: PetscFree(namelist);
311: }
312: PetscOptionsTail();
313: return 0;
314: }
316: static PetscErrorCode TSView_BasicSymplectic(TS ts,PetscViewer viewer)
317: {
318: return 0;
319: }
321: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts,PetscReal t,Vec X)
322: {
323: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
324: Vec update = bsymp->update;
325: PetscReal alpha = (ts->ptime - t)/ts->time_step;
327: VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
328: VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
329: return 0;
330: }
332: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
333: {
334: *yr = 1.0 + xr;
335: *yi = xi;
336: return 0;
337: }
339: /*@C
340: TSBasicSymplecticSetType - Set the type of the basic symplectic method
342: Logically Collective on TS
344: Input Parameters:
345: + ts - timestepping context
346: - bsymptype - type of the symplectic scheme
348: Options Database:
349: . -ts_basicsymplectic_type <scheme> - select the scheme
351: Notes:
352: 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
353: that is intended to store the user-provided RHS function.
355: Level: intermediate
356: @*/
357: PetscErrorCode TSBasicSymplecticSetType(TS ts,TSBasicSymplecticType bsymptype)
358: {
360: PetscTryMethod(ts,"TSBasicSymplecticSetType_C",(TS,TSBasicSymplecticType),(ts,bsymptype));
361: return 0;
362: }
364: /*@C
365: TSBasicSymplecticGetType - Get the type of the basic symplectic method
367: Logically Collective on TS
369: Input Parameters:
370: + ts - timestepping context
371: - bsymptype - type of the basic symplectic scheme
373: Level: intermediate
374: @*/
375: PetscErrorCode TSBasicSymplecticGetType(TS ts,TSBasicSymplecticType *bsymptype)
376: {
378: PetscUseMethod(ts,"TSBasicSymplecticGetType_C",(TS,TSBasicSymplecticType*),(ts,bsymptype));
379: return 0;
380: }
382: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts,TSBasicSymplecticType bsymptype)
383: {
384: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
385: BasicSymplecticSchemeLink link;
386: PetscBool match;
388: if (bsymp->scheme) {
389: PetscStrcmp(bsymp->scheme->name,bsymptype,&match);
390: if (match) return 0;
391: }
392: for (link = BasicSymplecticSchemeList; link; link=link->next) {
393: PetscStrcmp(link->sch.name,bsymptype,&match);
394: if (match) {
395: bsymp->scheme = &link->sch;
396: return 0;
397: }
398: }
399: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",bsymptype);
400: }
402: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts,TSBasicSymplecticType *bsymptype)
403: {
404: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
406: *bsymptype = bsymp->scheme->name;
407: return 0;
408: }
410: /*MC
411: TSBasicSymplectic - ODE solver using basic symplectic integration schemes
413: These methods are intened for separable Hamiltonian systems
415: $ qdot = dH(q,p,t)/dp
416: $ pdot = -dH(q,p,t)/dq
418: where the Hamiltonian can be split into the sum of kinetic energy and potential energy
420: $ H(q,p,t) = T(p,t) + V(q,t).
422: As a result, the system can be genearlly represented by
424: $ qdot = f(p,t) = dT(p,t)/dp
425: $ pdot = g(q,t) = -dV(q,t)/dq
427: and solved iteratively with
429: $ q_new = q_old + d_i*h*f(p_old,t_old)
430: $ t_new = t_old + d_i*h
431: $ p_new = p_old + c_i*h*g(p_new,t_new)
432: $ i=0,1,...,n.
434: 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.
435: 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.
437: Reference: wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)
439: Level: beginner
441: .seealso: TSCreate(), TSSetType(), TSRHSSplitSetIS(), TSRHSSplitSetRHSFunction()
443: M*/
444: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
445: {
446: TS_BasicSymplectic *bsymp;
448: TSBasicSymplecticInitializePackage();
449: PetscNewLog(ts,&bsymp);
450: ts->data = (void*)bsymp;
452: ts->ops->setup = TSSetUp_BasicSymplectic;
453: ts->ops->step = TSStep_BasicSymplectic;
454: ts->ops->reset = TSReset_BasicSymplectic;
455: ts->ops->destroy = TSDestroy_BasicSymplectic;
456: ts->ops->setfromoptions = TSSetFromOptions_BasicSymplectic;
457: ts->ops->view = TSView_BasicSymplectic;
458: ts->ops->interpolate = TSInterpolate_BasicSymplectic;
459: ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;
461: PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticSetType_C",TSBasicSymplecticSetType_BasicSymplectic);
462: PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticGetType_C",TSBasicSymplecticGetType_BasicSymplectic);
464: TSBasicSymplecticSetType(ts,TSBasicSymplecticDefault);
465: return 0;
466: }