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

 35:   Level: intermediate

 37: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`
 38: M*/

 40: /*MC
 41:   TSBASICSYMPLECTICVELVERLET - second order Velocity Verlet method (leapfrog method with starting process and determing velocity and position at the same time)

 43: Level: intermediate

 45: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`
 46: M*/

 48: /*@C
 49:   TSBasicSymplecticRegisterAll - Registers all of the basic symplectic integration methods in `TSBASICSYMPLECTIC`

 51:   Not Collective, but should be called by all processes which will need the schemes to be registered

 53:   Level: advanced

 55: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticRegisterDestroy()`
 56: @*/
 57: PetscErrorCode TSBasicSymplecticRegisterAll(void)
 58: {
 59:   if (TSBasicSymplecticRegisterAllCalled) return 0;
 60:   TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
 61:   {
 62:     PetscReal c[1] = {1.0}, d[1] = {1.0};
 63:     TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER, 1, 1, c, d);
 64:   }
 65:   {
 66:     PetscReal c[2] = {0, 1.0}, d[2] = {0.5, 0.5};
 67:     TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET, 2, 2, c, d);
 68:   }
 69:   {
 70:     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};
 71:     TSBasicSymplecticRegister(TSBASICSYMPLECTIC3, 3, 3, c, d);
 72:   }
 73:   {
 74: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
 75:     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};
 76:     TSBasicSymplecticRegister(TSBASICSYMPLECTIC4, 4, 4, c, d);
 77:   }
 78:   return 0;
 79: }

 81: /*@C
 82:    TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by `TSBasicSymplecticRegister()`.

 84:    Not Collective

 86:    Level: advanced

 88: .seealso: [](chapter_ts), `TSBasicSymplecticRegister()`, `TSBasicSymplecticRegisterAll()`, `TSBASICSYMPLECTIC`
 89: @*/
 90: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
 91: {
 92:   BasicSymplecticSchemeLink link;

 94:   while ((link = BasicSymplecticSchemeList)) {
 95:     BasicSymplecticScheme scheme = &link->sch;
 96:     BasicSymplecticSchemeList    = link->next;
 97:     PetscFree2(scheme->c, scheme->d);
 98:     PetscFree(scheme->name);
 99:     PetscFree(link);
100:   }
101:   TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
102:   return 0;
103: }

105: /*@C
106:   TSBasicSymplecticInitializePackage - This function initializes everything in the `TSBASICSYMPLECTIC` package. It is called
107:   from `TSInitializePackage()`.

109:   Level: developer

111: .seealso: [](chapter_ts), `PetscInitialize()`, `TSBASICSYMPLECTIC`
112: @*/
113: PetscErrorCode TSBasicSymplecticInitializePackage(void)
114: {
115:   if (TSBasicSymplecticPackageInitialized) return 0;
116:   TSBasicSymplecticPackageInitialized = PETSC_TRUE;
117:   TSBasicSymplecticRegisterAll();
118:   PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
119:   return 0;
120: }

122: /*@C
123:   TSBasicSymplecticFinalizePackage - This function destroys everything in the `TSBASICSYMPLECTIC` package. It is
124:   called from `PetscFinalize()`.

126:   Level: developer

128: .seealso: [](chapter_ts), `PetscFinalize()`, `TSBASICSYMPLECTIC`
129: @*/
130: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
131: {
132:   TSBasicSymplecticPackageInitialized = PETSC_FALSE;
133:   TSBasicSymplecticRegisterDestroy();
134:   return 0;
135: }

137: /*@C
138:    TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.

140:    Not Collective, but the same schemes should be registered on all processes on which they will be used

142:    Input Parameters:
143: +  name - identifier for method
144: .  order - approximation order of method
145: .  s - number of stages, this is the dimension of the matrices below
146: .  c - coefficients for updating generalized position (dimension s)
147: -  d - coefficients for updating generalized momentum (dimension s)

149:    Level: advanced

151:    Notes:
152:    Several symplectic methods are provided, this function is only needed to create new methods.

154: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`
155: @*/
156: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name, PetscInt order, PetscInt s, PetscReal c[], PetscReal d[])
157: {
158:   BasicSymplecticSchemeLink link;
159:   BasicSymplecticScheme     scheme;


165:   TSBasicSymplecticInitializePackage();
166:   PetscNew(&link);
167:   scheme = &link->sch;
168:   PetscStrallocpy(name, &scheme->name);
169:   scheme->order = order;
170:   scheme->s     = s;
171:   PetscMalloc2(s, &scheme->c, s, &scheme->d);
172:   PetscArraycpy(scheme->c, c, s);
173:   PetscArraycpy(scheme->d, d, s);
174:   link->next                = BasicSymplecticSchemeList;
175:   BasicSymplecticSchemeList = link;
176:   return 0;
177: }

179: /*
180: The simplified form of the equations are:

182: $ p_{i+1} = p_i + c_i*g(q_i)*h
183: $ q_{i+1} = q_i + d_i*f(p_{i+1},t_{i+1})*h

185: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.

187: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:

189: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
190: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
191: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
192: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2

194: */
195: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
196: {
197:   TS_BasicSymplectic   *bsymp    = (TS_BasicSymplectic *)ts->data;
198:   BasicSymplecticScheme scheme   = bsymp->scheme;
199:   Vec                   solution = ts->vec_sol, update = bsymp->update, q, p, q_update, p_update;
200:   IS                    is_q = bsymp->is_q, is_p = bsymp->is_p;
201:   TS                    subts_q = bsymp->subts_q, subts_p = bsymp->subts_p;
202:   PetscBool             stageok;
203:   PetscReal             next_time_step = ts->time_step;
204:   PetscInt              iter;

206:   VecGetSubVector(solution, is_q, &q);
207:   VecGetSubVector(solution, is_p, &p);
208:   VecGetSubVector(update, is_q, &q_update);
209:   VecGetSubVector(update, is_p, &p_update);

211:   for (iter = 0; iter < scheme->s; iter++) {
212:     TSPreStage(ts, ts->ptime);
213:     /* update velocity p */
214:     if (scheme->c[iter]) {
215:       TSComputeRHSFunction(subts_p, ts->ptime, q, p_update);
216:       VecAXPY(p, scheme->c[iter] * ts->time_step, p_update);
217:     }
218:     /* update position q */
219:     if (scheme->d[iter]) {
220:       TSComputeRHSFunction(subts_q, ts->ptime, p, q_update);
221:       VecAXPY(q, scheme->d[iter] * ts->time_step, q_update);
222:       ts->ptime = ts->ptime + scheme->d[iter] * ts->time_step;
223:     }
224:     TSPostStage(ts, ts->ptime, 0, &solution);
225:     TSAdaptCheckStage(ts->adapt, ts, ts->ptime, solution, &stageok);
226:     if (!stageok) {
227:       ts->reason = TS_DIVERGED_STEP_REJECTED;
228:       return 0;
229:     }
230:     TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok);
231:     if (!stageok) {
232:       ts->reason = TS_DIVERGED_STEP_REJECTED;
233:       return 0;
234:     }
235:   }

237:   ts->time_step = next_time_step;
238:   VecRestoreSubVector(solution, is_q, &q);
239:   VecRestoreSubVector(solution, is_p, &p);
240:   VecRestoreSubVector(update, is_q, &q_update);
241:   VecRestoreSubVector(update, is_p, &p_update);
242:   return 0;
243: }

245: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, void *ctx)
246: {
247:   return 0;
248: }

250: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
251: {
252:   return 0;
253: }

255: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm, DM subdm, void *ctx)
256: {
257:   return 0;
258: }

260: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
261: {
262:   return 0;
263: }

265: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
266: {
267:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
268:   DM                  dm;

270:   TSRHSSplitGetIS(ts, "position", &bsymp->is_q);
271:   TSRHSSplitGetIS(ts, "momentum", &bsymp->is_p);
273:   TSRHSSplitGetSubTS(ts, "position", &bsymp->subts_q);
274:   TSRHSSplitGetSubTS(ts, "momentum", &bsymp->subts_p);

277:   VecDuplicate(ts->vec_sol, &bsymp->update);

279:   TSGetAdapt(ts, &ts->adapt);
280:   TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
281:   TSGetDM(ts, &dm);
282:   if (dm) {
283:     DMCoarsenHookAdd(dm, DMCoarsenHook_BasicSymplectic, DMRestrictHook_BasicSymplectic, ts);
284:     DMSubDomainHookAdd(dm, DMSubDomainHook_BasicSymplectic, DMSubDomainRestrictHook_BasicSymplectic, ts);
285:   }
286:   return 0;
287: }

289: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
290: {
291:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

293:   VecDestroy(&bsymp->update);
294:   return 0;
295: }

297: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
298: {
299:   TSReset_BasicSymplectic(ts);
300:   PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", NULL);
301:   PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", NULL);
302:   PetscFree(ts->data);
303:   return 0;
304: }

306: static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems *PetscOptionsObject)
307: {
308:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

310:   PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
311:   {
312:     BasicSymplecticSchemeLink link;
313:     PetscInt                  count, choice;
314:     PetscBool                 flg;
315:     const char              **namelist;

317:     for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++)
318:       ;
319:     PetscMalloc1(count, (char ***)&namelist);
320:     for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++) namelist[count] = link->sch.name;
321:     PetscOptionsEList("-ts_basicsymplectic_type", "Family of basic symplectic integration method", "TSBasicSymplecticSetType", (const char *const *)namelist, count, bsymp->scheme->name, &choice, &flg);
322:     if (flg) TSBasicSymplecticSetType(ts, namelist[choice]);
323:     PetscFree(namelist);
324:   }
325:   PetscOptionsHeadEnd();
326:   return 0;
327: }

329: static PetscErrorCode TSView_BasicSymplectic(TS ts, PetscViewer viewer)
330: {
331:   return 0;
332: }

334: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts, PetscReal t, Vec X)
335: {
336:   TS_BasicSymplectic *bsymp  = (TS_BasicSymplectic *)ts->data;
337:   Vec                 update = bsymp->update;
338:   PetscReal           alpha  = (ts->ptime - t) / ts->time_step;

340:   VecWAXPY(X, -ts->time_step, update, ts->vec_sol);
341:   VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol);
342:   return 0;
343: }

345: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
346: {
347:   *yr = 1.0 + xr;
348:   *yi = xi;
349:   return 0;
350: }

352: /*@C
353:   TSBasicSymplecticSetType - Set the type of the basic symplectic method

355:   Logically Collective

357:   Input Parameters:
358: +  ts - timestepping context
359: -  bsymptype - type of the symplectic scheme

361:   Options Database Key:
362: .  -ts_basicsymplectic_type <scheme> - select the scheme

364:   Level: intermediate

366:   Note:
367:     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`
368:     that is intended to store the user-provided RHS function.

370: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
371: @*/
372: PetscErrorCode TSBasicSymplecticSetType(TS ts, TSBasicSymplecticType bsymptype)
373: {
375:   PetscTryMethod(ts, "TSBasicSymplecticSetType_C", (TS, TSBasicSymplecticType), (ts, bsymptype));
376:   return 0;
377: }

379: /*@C
380:   TSBasicSymplecticGetType - Get the type of the basic symplectic method

382:   Logically Collective

384:   Input Parameters:
385: +  ts - timestepping context
386: -  bsymptype - type of the basic symplectic scheme

388:   Level: intermediate

390: .seealso: [](chapter_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
391: @*/
392: PetscErrorCode TSBasicSymplecticGetType(TS ts, TSBasicSymplecticType *bsymptype)
393: {
395:   PetscUseMethod(ts, "TSBasicSymplecticGetType_C", (TS, TSBasicSymplecticType *), (ts, bsymptype));
396:   return 0;
397: }

399: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts, TSBasicSymplecticType bsymptype)
400: {
401:   TS_BasicSymplectic       *bsymp = (TS_BasicSymplectic *)ts->data;
402:   BasicSymplecticSchemeLink link;
403:   PetscBool                 match;

405:   if (bsymp->scheme) {
406:     PetscStrcmp(bsymp->scheme->name, bsymptype, &match);
407:     if (match) return 0;
408:   }
409:   for (link = BasicSymplecticSchemeList; link; link = link->next) {
410:     PetscStrcmp(link->sch.name, bsymptype, &match);
411:     if (match) {
412:       bsymp->scheme = &link->sch;
413:       return 0;
414:     }
415:   }
416:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", bsymptype);
417: }

419: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts, TSBasicSymplecticType *bsymptype)
420: {
421:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

423:   *bsymptype = bsymp->scheme->name;
424:   return 0;
425: }

427: /*MC
428:   TSBASICSYMPLECTIC - ODE solver using basic symplectic integration schemes

430:   These methods are intended for separable Hamiltonian systems
431: .vb
432:   qdot = dH(q,p,t)/dp
433:   pdot = -dH(q,p,t)/dq
434: .ve

436:   where the Hamiltonian can be split into the sum of kinetic energy and potential energy
437: .vb
438:   H(q,p,t) = T(p,t) + V(q,t).
439: .ve

441:   As a result, the system can be genearlly represented by
442: .vb
443:   qdot = f(p,t) = dT(p,t)/dp
444:   pdot = g(q,t) = -dV(q,t)/dq
445: .ve

447:   and solved iteratively with
448: .vb
449:   q_new = q_old + d_i*h*f(p_old,t_old)
450:   t_new = t_old + d_i*h
451:   p_new = p_old + c_i*h*g(p_new,t_new)
452:   i=0,1,...,n.
453: .ve

455:   The solution vector should contain both q and p, which correspond to (generalized) position and momentum respectively. Note that the momentum component
456:   could simply be velocity in some representations. The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum".
457:   Each split is associated with an `IS` object and a sub-`TS` that is intended to store the user-provided RHS function.

459:   Level: beginner

461:   Reference:
462: . * -  wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)

464: .seealso: [](chapter_ts), `TSCreate()`, `TSSetType()`, `TSRHSSplitSetIS()`, `TSRHSSplitSetRHSFunction()`, `TSType`
465: M*/
466: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
467: {
468:   TS_BasicSymplectic *bsymp;

470:   TSBasicSymplecticInitializePackage();
471:   PetscNew(&bsymp);
472:   ts->data = (void *)bsymp;

474:   ts->ops->setup           = TSSetUp_BasicSymplectic;
475:   ts->ops->step            = TSStep_BasicSymplectic;
476:   ts->ops->reset           = TSReset_BasicSymplectic;
477:   ts->ops->destroy         = TSDestroy_BasicSymplectic;
478:   ts->ops->setfromoptions  = TSSetFromOptions_BasicSymplectic;
479:   ts->ops->view            = TSView_BasicSymplectic;
480:   ts->ops->interpolate     = TSInterpolate_BasicSymplectic;
481:   ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;

483:   PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic);
484:   PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic);

486:   TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault);
487:   return 0;
488: }