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: }