Actual source code: basicsymplectic.c
petsc-3.13.6 2020-09-29
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: {
58: if (TSBasicSymplecticRegisterAllCalled) return(0);
59: TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
60: {
61: PetscReal c[1] = {1.0},d[1] = {1.0};
62: TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER,1,1,c,d);
63: }
64: {
65: PetscReal c[2] = {0,1.0},d[2] = {0.5,0.5};
66: TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET,2,2,c,d);
67: }
68: {
69: 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};
70: TSBasicSymplecticRegister(TSBASICSYMPLECTIC3,3,3,c,d);
71: }
72: {
73: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
74: 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};
75: TSBasicSymplecticRegister(TSBASICSYMPLECTIC4,4,4,c,d);
76: }
77: return(0);
78: }
80: /*@C
81: TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by TSBasicSymplecticRegister().
83: Not Collective
85: Level: advanced
87: .seealso: TSBasicSymplecticRegister(), TSBasicSymplecticRegisterAll()
88: @*/
89: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
90: {
91: PetscErrorCode ierr;
92: BasicSymplecticSchemeLink link;
95: while ((link = BasicSymplecticSchemeList)) {
96: BasicSymplecticScheme scheme = &link->sch;
97: BasicSymplecticSchemeList = link->next;
98: PetscFree2(scheme->c,scheme->d);
99: PetscFree(scheme->name);
100: PetscFree(link);
101: }
102: TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
103: return(0);
104: }
106: /*@C
107: TSBasicSymplecticInitializePackage - This function initializes everything in the TSBasicSymplectic package. It is called
108: from TSInitializePackage().
110: Level: developer
112: .seealso: PetscInitialize()
113: @*/
114: PetscErrorCode TSBasicSymplecticInitializePackage(void)
115: {
119: if (TSBasicSymplecticPackageInitialized) return(0);
120: TSBasicSymplecticPackageInitialized = PETSC_TRUE;
121: TSBasicSymplecticRegisterAll();
122: PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
123: return(0);
124: }
126: /*@C
127: TSBasicSymplecticFinalizePackage - This function destroys everything in the TSBasicSymplectic package. It is
128: called from PetscFinalize().
130: Level: developer
132: .seealso: PetscFinalize()
133: @*/
134: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
135: {
139: TSBasicSymplecticPackageInitialized = PETSC_FALSE;
140: TSBasicSymplecticRegisterDestroy();
141: return(0);
142: }
144: /*@C
145: TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.
147: Not Collective, but the same schemes should be registered on all processes on which they will be used
149: Input Parameters:
150: + name - identifier for method
151: . order - approximation order of method
152: . s - number of stages, this is the dimension of the matrices below
153: . c - coefficients for updating generalized position (dimension s)
154: - d - coefficients for updating generalized momentum (dimension s)
156: Notes:
157: Several symplectic methods are provided, this function is only needed to create new methods.
159: Level: advanced
161: .seealso: TSBasicSymplectic
162: @*/
163: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name,PetscInt order,PetscInt s,PetscReal c[],PetscReal d[])
164: {
165: BasicSymplecticSchemeLink link;
166: BasicSymplecticScheme scheme;
167: PetscErrorCode ierr;
174: TSBasicSymplecticInitializePackage();
175: PetscNew(&link);
176: scheme = &link->sch;
177: PetscStrallocpy(name,&scheme->name);
178: scheme->order = order;
179: scheme->s = s;
180: PetscMalloc2(s,&scheme->c,s,&scheme->d);
181: PetscArraycpy(scheme->c,c,s);
182: PetscArraycpy(scheme->d,d,s);
183: link->next = BasicSymplecticSchemeList;
184: BasicSymplecticSchemeList = link;
185: return(0);
186: }
188: /*
189: The simplified form of the equations are:
191: $ p_{i+1} = p_i + c_i*g(q_i)*h
192: $ q_{i+1} = q_i + d_i*f(p_{i+1},t_{i+1})*h
194: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.
196: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:
198: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
199: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
200: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
201: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2
203: */
204: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
205: {
206: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
207: BasicSymplecticScheme scheme = bsymp->scheme;
208: Vec solution = ts->vec_sol,update = bsymp->update,q,p,q_update,p_update;
209: IS is_q = bsymp->is_q,is_p = bsymp->is_p;
210: TS subts_q = bsymp->subts_q,subts_p = bsymp->subts_p;
211: PetscBool stageok;
212: PetscReal next_time_step = ts->time_step;
213: PetscInt iter;
214: PetscErrorCode ierr;
217: VecGetSubVector(solution,is_q,&q);
218: VecGetSubVector(solution,is_p,&p);
219: VecGetSubVector(update,is_q,&q_update);
220: VecGetSubVector(update,is_p,&p_update);
222: for (iter = 0;iter<scheme->s;iter++) {
223: TSPreStage(ts,ts->ptime);
224: /* update velocity p */
225: if (scheme->c[iter]) {
226: TSComputeRHSFunction(subts_p,ts->ptime,q,p_update);
227: VecAXPY(p,scheme->c[iter]*ts->time_step,p_update);
228: }
229: /* update position q */
230: if (scheme->d[iter]) {
231: TSComputeRHSFunction(subts_q,ts->ptime,p,q_update);
232: VecAXPY(q,scheme->d[iter]*ts->time_step,q_update);
233: ts->ptime = ts->ptime+scheme->d[iter]*ts->time_step;
234: }
235: TSPostStage(ts,ts->ptime,0,&solution);
236: TSAdaptCheckStage(ts->adapt,ts,ts->ptime,solution,&stageok);
237: if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
238: TSFunctionDomainError(ts,ts->ptime+ts->time_step,update,&stageok);
239: if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
240: }
242: ts->time_step = next_time_step;
243: VecRestoreSubVector(solution,is_q,&q);
244: VecRestoreSubVector(solution,is_p,&p);
245: VecRestoreSubVector(update,is_q,&q_update);
246: VecRestoreSubVector(update,is_p,&p_update);
247: return(0);
248: }
250: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine,DM coarse,void *ctx)
251: {
253: return(0);
254: }
256: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
257: {
259: return(0);
260: }
262: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm,DM subdm,void *ctx)
263: {
265: return(0);
266: }
268: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
269: {
272: return(0);
273: }
275: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
276: {
277: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
278: DM dm;
279: PetscErrorCode ierr;
282: TSRHSSplitGetIS(ts,"position",&bsymp->is_q);
283: TSRHSSplitGetIS(ts,"momentum",&bsymp->is_p);
284: if (!bsymp->is_q || !bsymp->is_p) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names positon and momentum respectively in order to use -ts_type basicsymplectic");
285: TSRHSSplitGetSubTS(ts,"position",&bsymp->subts_q);
286: TSRHSSplitGetSubTS(ts,"momentum",&bsymp->subts_p);
287: if (!bsymp->subts_q || !bsymp->subts_p) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up the RHSFunctions for position and momentum using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS");
289: VecDuplicate(ts->vec_sol,&bsymp->update);
291: TSGetAdapt(ts,&ts->adapt);
292: TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
293: TSGetDM(ts,&dm);
294: if (dm) {
295: DMCoarsenHookAdd(dm,DMCoarsenHook_BasicSymplectic,DMRestrictHook_BasicSymplectic,ts);
296: DMSubDomainHookAdd(dm,DMSubDomainHook_BasicSymplectic,DMSubDomainRestrictHook_BasicSymplectic,ts);
297: }
298: return(0);
299: }
301: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
302: {
303: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
304: PetscErrorCode ierr;
307: VecDestroy(&bsymp->update);
308: return(0);
309: }
311: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
312: {
316: TSReset_BasicSymplectic(ts);
317: PetscFree(ts->data);
318: return(0);
319: }
321: static PetscErrorCode TSSetFromOptions_BasicSymplectic(PetscOptionItems *PetscOptionsObject,TS ts)
322: {
323: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
324: PetscErrorCode ierr;
327: PetscOptionsHead(PetscOptionsObject,"Basic symplectic integrator options");
328: {
329: BasicSymplecticSchemeLink link;
330: PetscInt count,choice;
331: PetscBool flg;
332: const char **namelist;
334: for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) ;
335: PetscMalloc1(count,(char***)&namelist);
336: for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) namelist[count] = link->sch.name;
337: PetscOptionsEList("-ts_basicsymplectic_type","Family of basic symplectic integration method","TSBasicSymplecticSetType",(const char*const*)namelist,count,bsymp->scheme->name,&choice,&flg);
338: if (flg) {TSBasicSymplecticSetType(ts,namelist[choice]);}
339: PetscFree(namelist);
340: }
341: PetscOptionsTail();
342: return(0);
343: }
345: static PetscErrorCode TSView_BasicSymplectic(TS ts,PetscViewer viewer)
346: {
348: return(0);
349: }
351: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts,PetscReal t,Vec X)
352: {
353: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
354: Vec update = bsymp->update;
355: PetscReal alpha = (ts->ptime - t)/ts->time_step;
356: PetscErrorCode ierr;
359: VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
360: VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
361: return(0);
362: }
364: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
365: {
367: *yr = 1.0 + xr;
368: *yi = xi;
369: return(0);
370: }
372: /*@C
373: TSBasicSymplecticSetType - Set the type of the basic symplectic method
375: Logically Collective on TS
377: Input Parameter:
378: + ts - timestepping context
379: - bsymptype - type of the symplectic scheme
381: Options Database:
382: . -ts_basicsymplectic_type <scheme>
384: Notes:
385: 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.
387: Level: intermediate
388: @*/
389: PetscErrorCode TSBasicSymplecticSetType(TS ts,TSBasicSymplecticType bsymptype)
390: {
395: PetscTryMethod(ts,"TSBasicSymplecticSetType_C",(TS,TSBasicSymplecticType),(ts,bsymptype));
396: return(0);
397: }
399: /*@C
400: TSBasicSymplecticGetType - Get the type of the basic symplectic method
402: Logically Collective on TS
404: Input Parameter:
405: + ts - timestepping context
406: - bsymptype - type of the basic symplectic scheme
408: Level: intermediate
409: @*/
410: PetscErrorCode TSBasicSymplecticGetType(TS ts,TSBasicSymplecticType *bsymptype)
411: {
416: PetscUseMethod(ts,"TSBasicSymplecticGetType_C",(TS,TSBasicSymplecticType*),(ts,bsymptype));
417: return(0);
418: }
420: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts,TSBasicSymplecticType bsymptype)
421: {
422: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
423: BasicSymplecticSchemeLink link;
424: PetscBool match;
425: PetscErrorCode ierr;
428: if (bsymp->scheme) {
429: PetscStrcmp(bsymp->scheme->name,bsymptype,&match);
430: if (match) return(0);
431: }
432: for (link = BasicSymplecticSchemeList; link; link=link->next) {
433: PetscStrcmp(link->sch.name,bsymptype,&match);
434: if (match) {
435: bsymp->scheme = &link->sch;
436: return(0);
437: }
438: }
439: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",bsymptype);
440: return(0);
441: }
443: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts,TSBasicSymplecticType *bsymptype)
444: {
445: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
448: *bsymptype = bsymp->scheme->name;
449: return(0);
450: }
452: /*MC
453: TSBasicSymplectic - ODE solver using basic symplectic integration schemes
455: These methods are intened for separable Hamiltonian systems
457: $ qdot = dH(q,p,t)/dp
458: $ pdot = -dH(q,p,t)/dq
460: where the Hamiltonian can be split into the sum of kinetic energy and potential energy
462: $ H(q,p,t) = T(p,t) + V(q,t).
464: As a result, the system can be genearlly represented by
466: $ qdot = f(p,t) = dT(p,t)/dp
467: $ pdot = g(q,t) = -dV(q,t)/dq
469: and solved iteratively with
471: $ q_new = q_old + d_i*h*f(p_old,t_old)
472: $ t_new = t_old + d_i*h
473: $ p_new = p_old + c_i*h*g(p_new,t_new)
474: $ i=0,1,...,n.
476: 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.
477: 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.
479: Reference: wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)
481: Level: beginner
483: .seealso: TSCreate(), TSSetType(), TSRHSSplitSetIS(), TSRHSSplitSetRHSFunction()
485: M*/
486: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
487: {
488: TS_BasicSymplectic *bsymp;
489: PetscErrorCode ierr;
492: TSBasicSymplecticInitializePackage();
493: PetscNewLog(ts,&bsymp);
494: ts->data = (void*)bsymp;
496: ts->ops->setup = TSSetUp_BasicSymplectic;
497: ts->ops->step = TSStep_BasicSymplectic;
498: ts->ops->reset = TSReset_BasicSymplectic;
499: ts->ops->destroy = TSDestroy_BasicSymplectic;
500: ts->ops->setfromoptions = TSSetFromOptions_BasicSymplectic;
501: ts->ops->view = TSView_BasicSymplectic;
502: ts->ops->interpolate = TSInterpolate_BasicSymplectic;
503: ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;
505: PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticSetType_C",TSBasicSymplecticSetType_BasicSymplectic);
506: PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticGetType_C",TSBasicSymplecticGetType_BasicSymplectic);
508: TSBasicSymplecticSetType(ts,TSBasicSymplecticDefault);
509: return(0);
510: }