Actual source code: basicsymplectic.c
petsc-3.10.5 2019-03-28
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: .keywords: TS, TSBASICSYMPLECTIC, register, all
53: .seealso: TSBasicSymplecticRegisterDestroy()
54: @*/
55: PetscErrorCode TSBasicSymplecticRegisterAll(void)
56: {
60: if (TSBasicSymplecticRegisterAllCalled) return(0);
61: TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
62: {
63: PetscReal c[1] = {1.0},d[1] = {1.0};
64: TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER,1,1,c,d);
65: }
66: {
67: PetscReal c[2] = {0,1.0},d[2] = {0.5,0.5};
68: 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: TSBasicSymplecticRegister(TSBASICSYMPLECTIC3,3,3,c,d);
73: }
74: {
75: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
76: 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};
77: TSBasicSymplecticRegister(TSBASICSYMPLECTIC4,4,4,c,d);
78: }
79: return(0);
80: }
82: /*@C
83: TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by TSBasicSymplecticRegister().
85: Not Collective
87: Level: advanced
89: .keywords: TSBasicSymplectic, register, destroy
90: .seealso: TSBasicSymplecticRegister(), TSBasicSymplecticRegisterAll()
91: @*/
92: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
93: {
94: PetscErrorCode ierr;
95: BasicSymplecticSchemeLink link;
98: while ((link = BasicSymplecticSchemeList)) {
99: BasicSymplecticScheme scheme = &link->sch;
100: BasicSymplecticSchemeList = link->next;
101: PetscFree2(scheme->c,scheme->d);
102: PetscFree(scheme->name);
103: PetscFree(link);
104: }
105: TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
106: return(0);
107: }
109: /*@C
110: TSBasicSymplecticInitializePackage - This function initializes everything in the TSBasicSymplectic package. It is called
111: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_BasicSymplectic()
112: when using static libraries.
114: Level: developer
116: .keywords: TS, TSBasicSymplectic, initialize, package
117: .seealso: PetscInitialize()
118: @*/
119: PetscErrorCode TSBasicSymplecticInitializePackage(void)
120: {
124: if (TSBasicSymplecticPackageInitialized) return(0);
125: TSBasicSymplecticPackageInitialized = PETSC_TRUE;
126: TSBasicSymplecticRegisterAll();
127: PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
128: return(0);
129: }
131: /*@C
132: TSBasicSymplecticFinalizePackage - This function destroys everything in the TSBasicSymplectic package. It is
133: called from PetscFinalize().
135: Level: developer
137: .keywords: Petsc, destroy, package
138: .seealso: PetscFinalize()
139: @*/
140: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
141: {
145: TSBasicSymplecticPackageInitialized = PETSC_FALSE;
146: TSBasicSymplecticRegisterDestroy();
147: return(0);
148: }
150: /*@C
151: TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.
153: Not Collective, but the same schemes should be registered on all processes on which they will be used
155: Input Parameters:
156: + name - identifier for method
157: . order - approximation order of method
158: . s - number of stages, this is the dimension of the matrices below
159: . c - coefficients for updating generalized position (dimension s)
160: - d - coefficients for updating generalized momentum (dimension s)
162: Notes:
163: Several symplectic methods are provided, this function is only needed to create new methods.
165: Level: advanced
167: .keywords: TS, register
169: .seealso: TSBasicSymplectic
170: @*/
171: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name,PetscInt order,PetscInt s,PetscReal c[],PetscReal d[])
172: {
173: BasicSymplecticSchemeLink link;
174: BasicSymplecticScheme scheme;
175: PetscErrorCode ierr;
182: TSBasicSymplecticInitializePackage();
183: PetscCalloc1(1,&link);
184: scheme = &link->sch;
185: PetscStrallocpy(name,&scheme->name);
186: scheme->order = order;
187: scheme->s = s;
188: PetscMalloc2(s,&scheme->c,s,&scheme->d);
189: PetscMemcpy(scheme->c,c,s*sizeof(c[0]));
190: PetscMemcpy(scheme->d,d,s*sizeof(d[0]));
191: link->next = BasicSymplecticSchemeList;
192: BasicSymplecticSchemeList = link;
193: return(0);
194: }
196: /*
197: The simplified form of the equations are:
199: $ p_{i+1} = p_i + c_i*g(q_i)*h
200: $ q_{i+1} = q_i + d_i*f(p_{i+1},t_{i+1})*h
202: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.
204: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:
206: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
207: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
208: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
209: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2
211: */
212: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
213: {
214: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
215: BasicSymplecticScheme scheme = bsymp->scheme;
216: Vec solution = ts->vec_sol,update = bsymp->update,q,p,q_update,p_update;
217: IS is_q = bsymp->is_q,is_p = bsymp->is_p;
218: TS subts_q = bsymp->subts_q,subts_p = bsymp->subts_p;
219: PetscBool stageok;
220: PetscReal next_time_step = ts->time_step;
221: PetscInt iter;
222: PetscErrorCode ierr;
225: VecGetSubVector(solution,is_q,&q);
226: VecGetSubVector(solution,is_p,&p);
227: VecGetSubVector(update,is_q,&q_update);
228: VecGetSubVector(update,is_p,&p_update);
230: for (iter = 0;iter<scheme->s;iter++) {
231: TSPreStage(ts,ts->ptime);
232: /* update velocity p */
233: if (scheme->c[iter]) {
234: TSComputeRHSFunction(subts_p,ts->ptime,q,p_update);
235: VecAXPY(p,scheme->c[iter]*ts->time_step,p_update);
236: }
237: /* update position q */
238: if (scheme->d[iter]) {
239: TSComputeRHSFunction(subts_q,ts->ptime,p,q_update);
240: VecAXPY(q,scheme->d[iter]*ts->time_step,q_update);
241: ts->ptime = ts->ptime+scheme->d[iter]*ts->time_step;
242: }
243: TSPostStage(ts,ts->ptime,0,&solution);
244: TSAdaptCheckStage(ts->adapt,ts,ts->ptime,solution,&stageok);
245: if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
246: TSFunctionDomainError(ts,ts->ptime+ts->time_step,update,&stageok);
247: if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
248: }
250: ts->time_step = next_time_step;
251: VecRestoreSubVector(solution,is_q,&q);
252: VecRestoreSubVector(solution,is_p,&p);
253: VecRestoreSubVector(update,is_q,&q_update);
254: VecRestoreSubVector(update,is_p,&p_update);
255: return(0);
256: }
258: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine,DM coarse,void *ctx)
259: {
261: return(0);
262: }
264: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
265: {
267: return(0);
268: }
270: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm,DM subdm,void *ctx)
271: {
273: return(0);
274: }
276: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
277: {
280: return(0);
281: }
283: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
284: {
285: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
286: DM dm;
287: PetscErrorCode ierr;
290: TSRHSSplitGetIS(ts,"position",&bsymp->is_q);
291: TSRHSSplitGetIS(ts,"momentum",&bsymp->is_p);
292: 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");
293: TSRHSSplitGetSubTS(ts,"position",&bsymp->subts_q);
294: TSRHSSplitGetSubTS(ts,"momentum",&bsymp->subts_p);
295: 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");
297: VecDuplicate(ts->vec_sol,&bsymp->update);
299: TSGetAdapt(ts,&ts->adapt);
300: TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
301: TSGetDM(ts,&dm);
302: if (dm) {
303: DMCoarsenHookAdd(dm,DMCoarsenHook_BasicSymplectic,DMRestrictHook_BasicSymplectic,ts);
304: DMSubDomainHookAdd(dm,DMSubDomainHook_BasicSymplectic,DMSubDomainRestrictHook_BasicSymplectic,ts);
305: }
306: return(0);
307: }
309: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
310: {
311: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
312: PetscErrorCode ierr;
315: VecDestroy(&bsymp->update);
316: return(0);
317: }
319: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
320: {
324: TSReset_BasicSymplectic(ts);
325: PetscFree(ts->data);
326: return(0);
327: }
329: static PetscErrorCode TSSetFromOptions_BasicSymplectic(PetscOptionItems *PetscOptionsObject,TS ts)
330: {
331: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
332: PetscErrorCode ierr;
335: PetscOptionsHead(PetscOptionsObject,"Basic symplectic integrator options");
336: {
337: BasicSymplecticSchemeLink link;
338: PetscInt count,choice;
339: PetscBool flg;
340: const char **namelist;
342: for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) ;
343: PetscMalloc1(count,(char***)&namelist);
344: for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) namelist[count] = link->sch.name;
345: PetscOptionsEList("-ts_basicsymplectic_type","Family of basic symplectic integration method","TSBasicSymplecticSetType",(const char*const*)namelist,count,bsymp->scheme->name,&choice,&flg);
346: if (flg) {TSBasicSymplecticSetType(ts,namelist[choice]);}
347: PetscFree(namelist);
348: }
349: PetscOptionsTail();
350: return(0);
351: }
353: static PetscErrorCode TSView_BasicSymplectic(TS ts,PetscViewer viewer)
354: {
356: return(0);
357: }
359: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts,PetscReal t,Vec X)
360: {
361: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
362: Vec update = bsymp->update;
363: PetscReal alpha = (ts->ptime - t)/ts->time_step;
364: PetscErrorCode ierr;
367: VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
368: VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
369: return(0);
370: }
372: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
373: {
375: *yr = 1.0 + xr;
376: *yi = xi;
377: return(0);
378: }
380: /*@C
381: TSBasicSymplecticSetType - Set the type of the basic symplectic method
383: Logically Collective on TS
385: Input Parameter:
386: + ts - timestepping context
387: - bsymptype - type of the symplectic scheme
389: Options Database:
390: . -ts_basicsymplectic_type <scheme>
392: Notes:
393: 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.
395: Level: intermediate
396: @*/
397: PetscErrorCode TSBasicSymplecticSetType(TS ts,TSBasicSymplecticType bsymptype)
398: {
403: PetscTryMethod(ts,"TSBasicSymplecticSetType_C",(TS,TSBasicSymplecticType),(ts,bsymptype));
404: return(0);
405: }
407: /*@C
408: TSBasicSymplecticGetType - Get the type of the basic symplectic method
410: Logically Collective on TS
412: Input Parameter:
413: + ts - timestepping context
414: - bsymptype - type of the basic symplectic scheme
416: Level: intermediate
417: @*/
418: PetscErrorCode TSBasicSymplecticGetType(TS ts,TSBasicSymplecticType *bsymptype)
419: {
424: PetscUseMethod(ts,"TSBasicSymplecticGetType_C",(TS,TSBasicSymplecticType*),(ts,bsymptype));
425: return(0);
426: }
428: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts,TSBasicSymplecticType bsymptype)
429: {
430: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
431: BasicSymplecticSchemeLink link;
432: PetscBool match;
433: PetscErrorCode ierr;
436: if (bsymp->scheme) {
437: PetscStrcmp(bsymp->scheme->name,bsymptype,&match);
438: if (match) return(0);
439: }
440: for (link = BasicSymplecticSchemeList; link; link=link->next) {
441: PetscStrcmp(link->sch.name,bsymptype,&match);
442: if (match) {
443: bsymp->scheme = &link->sch;
444: return(0);
445: }
446: }
447: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",bsymptype);
448: return(0);
449: }
451: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts,TSBasicSymplecticType *bsymptype)
452: {
453: TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
456: *bsymptype = bsymp->scheme->name;
457: return(0);
458: }
460: /*MC
461: TSBasicSymplectic - ODE solver using basic symplectic integration schemes
463: These methods are intened for separable Hamiltonian systems
465: $ qdot = dH(q,p,t)/dp
466: $ pdot = -dH(q,p,t)/dq
468: where the Hamiltonian can be split into the sum of kinetic energy and potential energy
470: $ H(q,p,t) = T(p,t) + V(q,t).
472: As a result, the system can be genearlly represented by
474: $ qdot = f(p,t) = dT(p,t)/dp
475: $ pdot = g(q,t) = -dV(q,t)/dq
477: and solved iteratively with
479: $ q_new = q_old + d_i*h*f(p_old,t_old)
480: $ t_new = t_old + d_i*h
481: $ p_new = p_old + c_i*h*g(p_new,t_new)
482: $ i=0,1,...,n.
484: 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.
485: 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.
487: Reference: wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)
489: Level: beginner
491: .seealso: TSCreate(), TSSetType(), TSRHSSplitSetIS(), TSRHSSplitSetRHSFunction()
493: M*/
494: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
495: {
496: TS_BasicSymplectic *bsymp;
497: PetscErrorCode ierr;
500: TSBasicSymplecticInitializePackage();
501: PetscNewLog(ts,&bsymp);
502: ts->data = (void*)bsymp;
504: ts->ops->setup = TSSetUp_BasicSymplectic;
505: ts->ops->step = TSStep_BasicSymplectic;
506: ts->ops->reset = TSReset_BasicSymplectic;
507: ts->ops->destroy = TSDestroy_BasicSymplectic;
508: ts->ops->setfromoptions = TSSetFromOptions_BasicSymplectic;
509: ts->ops->view = TSView_BasicSymplectic;
510: ts->ops->interpolate = TSInterpolate_BasicSymplectic;
511: ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;
513: PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticSetType_C",TSBasicSymplecticSetType_BasicSymplectic);
514: PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticGetType_C",TSBasicSymplecticGetType_BasicSymplectic);
516: TSBasicSymplecticSetType(ts,TSBasicSymplecticDefault);
517: return(0);
518: }