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: [](ch_ts), `TSBASICSYMPLECTIC`
 38: M*/

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

 43: Level: intermediate

 45: .seealso: [](ch_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: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticRegisterDestroy()`
 56: @*/
 57: PetscErrorCode TSBasicSymplecticRegisterAll(void)
 58: {
 59:   PetscFunctionBegin;
 60:   if (TSBasicSymplecticRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
 61:   TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
 62:   {
 63:     PetscReal c[1] = {1.0}, d[1] = {1.0};
 64:     PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER, 1, 1, c, d));
 65:   }
 66:   {
 67:     PetscReal c[2] = {0, 1.0}, d[2] = {0.5, 0.5};
 68:     PetscCall(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:     PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTIC3, 3, 3, c, d));
 73:   }
 74:   {
 75: #define CUBEROOTOFTWO 1.2599210498948731647672106
 76:     PetscReal c[4] = {1.0 / 2.0 / (2.0 - CUBEROOTOFTWO), (1.0 - CUBEROOTOFTWO) / 2.0 / (2.0 - CUBEROOTOFTWO), (1.0 - CUBEROOTOFTWO) / 2.0 / (2.0 - CUBEROOTOFTWO), 1.0 / 2.0 / (2.0 - CUBEROOTOFTWO)}, d[4] = {1.0 / (2.0 - CUBEROOTOFTWO), -CUBEROOTOFTWO / (2.0 - CUBEROOTOFTWO), 1.0 / (2.0 - CUBEROOTOFTWO), 0};
 77:     PetscCall(TSBasicSymplecticRegister(TSBASICSYMPLECTIC4, 4, 4, c, d));
 78:   }
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

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

 85:   Not Collective

 87:   Level: advanced

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

 95:   PetscFunctionBegin;
 96:   while ((link = BasicSymplecticSchemeList)) {
 97:     BasicSymplecticScheme scheme = &link->sch;
 98:     BasicSymplecticSchemeList    = link->next;
 99:     PetscCall(PetscFree2(scheme->c, scheme->d));
100:     PetscCall(PetscFree(scheme->name));
101:     PetscCall(PetscFree(link));
102:   }
103:   TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

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

111:   Level: developer

113: .seealso: [](ch_ts), `PetscInitialize()`, `TSBASICSYMPLECTIC`
114: @*/
115: PetscErrorCode TSBasicSymplecticInitializePackage(void)
116: {
117:   PetscFunctionBegin;
118:   if (TSBasicSymplecticPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
119:   TSBasicSymplecticPackageInitialized = PETSC_TRUE;
120:   PetscCall(TSBasicSymplecticRegisterAll());
121:   PetscCall(PetscRegisterFinalize(TSBasicSymplecticFinalizePackage));
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

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

129:   Level: developer

131: .seealso: [](ch_ts), `PetscFinalize()`, `TSBASICSYMPLECTIC`
132: @*/
133: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
134: {
135:   PetscFunctionBegin;
136:   TSBasicSymplecticPackageInitialized = PETSC_FALSE;
137:   PetscCall(TSBasicSymplecticRegisterDestroy());
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

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

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

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

153:   Level: advanced

155:   Note:
156:   Several symplectic methods are provided, this function is only needed to create new methods.

158: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
159: @*/
160: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name, PetscInt order, PetscInt s, PetscReal c[], PetscReal d[])
161: {
162:   BasicSymplecticSchemeLink link;
163:   BasicSymplecticScheme     scheme;

165:   PetscFunctionBegin;
166:   PetscAssertPointer(name, 1);
167:   PetscAssertPointer(c, 4);
168:   PetscAssertPointer(d, 5);

170:   PetscCall(TSBasicSymplecticInitializePackage());
171:   PetscCall(PetscNew(&link));
172:   scheme = &link->sch;
173:   PetscCall(PetscStrallocpy(name, &scheme->name));
174:   scheme->order = order;
175:   scheme->s     = s;
176:   PetscCall(PetscMalloc2(s, &scheme->c, s, &scheme->d));
177:   PetscCall(PetscArraycpy(scheme->c, c, s));
178:   PetscCall(PetscArraycpy(scheme->d, d, s));
179:   link->next                = BasicSymplecticSchemeList;
180:   BasicSymplecticSchemeList = link;
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: /*
185: The simplified form of the equations are:

187: .vb
188:  q_{i+1} = q_i + c_i*g(p_i)*h
189:  p_{i+1} = p_i + d_i*f(q_{i+1})*h
190: .ve

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

194: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:
195: .vb
196: - Update the position of the particle by adding to it its velocity multiplied by c_1
197: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by d_1
198: - Update the position of the particle by adding to it its (updated) velocity multiplied by c_2
199: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by d_2
200: .ve

202: */
203: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
204: {
205:   TS_BasicSymplectic   *bsymp    = (TS_BasicSymplectic *)ts->data;
206:   BasicSymplecticScheme scheme   = bsymp->scheme;
207:   Vec                   solution = ts->vec_sol, update = bsymp->update, q, p, q_update, p_update;
208:   IS                    is_q = bsymp->is_q, is_p = bsymp->is_p;
209:   TS                    subts_q = bsymp->subts_q, subts_p = bsymp->subts_p;
210:   PetscBool             stageok = PETSC_TRUE;
211:   PetscReal             ptime = ts->ptime, next_time_step = ts->time_step;
212:   PetscInt              iter;

214:   PetscFunctionBegin;
215:   PetscCall(VecGetSubVector(update, is_q, &q_update));
216:   PetscCall(VecGetSubVector(update, is_p, &p_update));
217:   for (iter = 0; iter < scheme->s; iter++) {
218:     PetscCall(TSPreStage(ts, ptime));
219:     PetscCall(VecGetSubVector(solution, is_q, &q));
220:     PetscCall(VecGetSubVector(solution, is_p, &p));
221:     /* update position q */
222:     if (scheme->c[iter]) {
223:       PetscCall(TSComputeRHSFunction(subts_q, ptime, p, q_update));
224:       PetscCall(VecAXPY(q, scheme->c[iter] * ts->time_step, q_update));
225:     }
226:     /* update velocity p */
227:     if (scheme->d[iter]) {
228:       ptime = ptime + scheme->d[iter] * ts->time_step;
229:       PetscCall(TSComputeRHSFunction(subts_p, ptime, q, p_update));
230:       PetscCall(VecAXPY(p, scheme->d[iter] * ts->time_step, p_update));
231:     }
232:     PetscCall(VecRestoreSubVector(solution, is_q, &q));
233:     PetscCall(VecRestoreSubVector(solution, is_p, &p));
234:     PetscCall(TSPostStage(ts, ptime, 0, &solution));
235:     PetscCall(TSAdaptCheckStage(ts->adapt, ts, ptime, solution, &stageok));
236:     if (!stageok) goto finally;
237:     PetscCall(TSFunctionDomainError(ts, ptime, solution, &stageok));
238:     if (!stageok) goto finally;
239:   }

241: finally:
242:   if (!stageok) ts->reason = TS_DIVERGED_STEP_REJECTED;
243:   else ts->ptime += next_time_step;
244:   PetscCall(VecRestoreSubVector(update, is_q, &q_update));
245:   PetscCall(VecRestoreSubVector(update, is_p, &p_update));
246:   PetscFunctionReturn(PETSC_SUCCESS);
247: }

249: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, void *ctx)
250: {
251:   PetscFunctionBegin;
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
256: {
257:   PetscFunctionBegin;
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm, DM subdm, void *ctx)
262: {
263:   PetscFunctionBegin;
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
268: {
269:   PetscFunctionBegin;
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
274: {
275:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
276:   DM                  dm;

278:   PetscFunctionBegin;
279:   PetscCall(TSRHSSplitGetIS(ts, "position", &bsymp->is_q));
280:   PetscCall(TSRHSSplitGetIS(ts, "momentum", &bsymp->is_p));
281:   PetscCheck(bsymp->is_q && bsymp->is_p, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up RHSSplits with TSRHSSplitSetIS() using split names position and momentum respectively in order to use -ts_type basicsymplectic");
282:   PetscCall(TSRHSSplitGetSubTS(ts, "position", &bsymp->subts_q));
283:   PetscCall(TSRHSSplitGetSubTS(ts, "momentum", &bsymp->subts_p));
284:   PetscCheck(bsymp->subts_q && bsymp->subts_p, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up the RHSFunctions for position and momentum using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS");

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

288:   PetscCall(TSGetAdapt(ts, &ts->adapt));
289:   PetscCall(TSAdaptCandidatesClear(ts->adapt)); /* make sure to use fixed time stepping */
290:   PetscCall(TSGetDM(ts, &dm));
291:   if (dm) {
292:     PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_BasicSymplectic, DMRestrictHook_BasicSymplectic, ts));
293:     PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_BasicSymplectic, DMSubDomainRestrictHook_BasicSymplectic, ts));
294:   }
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
299: {
300:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

302:   PetscFunctionBegin;
303:   PetscCall(VecDestroy(&bsymp->update));
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
308: {
309:   PetscFunctionBegin;
310:   PetscCall(TSReset_BasicSymplectic(ts));
311:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", NULL));
312:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", NULL));
313:   PetscCall(PetscFree(ts->data));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems *PetscOptionsObject)
318: {
319:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

321:   PetscFunctionBegin;
322:   PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
323:   {
324:     BasicSymplecticSchemeLink link;
325:     PetscInt                  count, choice;
326:     PetscBool                 flg;
327:     const char              **namelist;

329:     for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++);
330:     PetscCall(PetscMalloc1(count, (char ***)&namelist));
331:     for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++) namelist[count] = link->sch.name;
332:     PetscCall(PetscOptionsEList("-ts_basicsymplectic_type", "Family of basic symplectic integration method", "TSBasicSymplecticSetType", (const char *const *)namelist, count, bsymp->scheme->name, &choice, &flg));
333:     if (flg) PetscCall(TSBasicSymplecticSetType(ts, namelist[choice]));
334:     PetscCall(PetscFree(namelist));
335:   }
336:   PetscOptionsHeadEnd();
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: static PetscErrorCode TSView_BasicSymplectic(TS ts, PetscViewer viewer)
341: {
342:   PetscFunctionBegin;
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }

346: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts, PetscReal t, Vec X)
347: {
348:   TS_BasicSymplectic *bsymp  = (TS_BasicSymplectic *)ts->data;
349:   Vec                 update = bsymp->update;
350:   PetscReal           alpha  = (ts->ptime - t) / ts->time_step;

352:   PetscFunctionBegin;
353:   PetscCall(VecWAXPY(X, -ts->time_step, update, ts->vec_sol));
354:   PetscCall(VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol));
355:   PetscFunctionReturn(PETSC_SUCCESS);
356: }

358: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
359: {
360:   PetscFunctionBegin;
361:   *yr = 1.0 + xr;
362:   *yi = xi;
363:   PetscFunctionReturn(PETSC_SUCCESS);
364: }

366: /*@C
367:   TSBasicSymplecticSetType - Set the type of the basic symplectic method

369:   Logically Collective

371:   Input Parameters:
372: + ts        - timestepping context
373: - bsymptype - type of the symplectic scheme

375:   Options Database Key:
376: . -ts_basicsymplectic_type <scheme> - select the scheme

378:   Level: intermediate

380:   Note:
381:   The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum".
382:   Each split is associated with an `IS` object and a sub-`TS`
383:   that is intended to store the user-provided RHS function.

385: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`
386: @*/
387: PetscErrorCode TSBasicSymplecticSetType(TS ts, TSBasicSymplecticType bsymptype)
388: {
389:   PetscFunctionBegin;
391:   PetscTryMethod(ts, "TSBasicSymplecticSetType_C", (TS, TSBasicSymplecticType), (ts, bsymptype));
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: /*@C
396:   TSBasicSymplecticGetType - Get the type of the basic symplectic method

398:   Logically Collective

400:   Input Parameters:
401: + ts        - timestepping context
402: - bsymptype - type of the basic symplectic scheme

404:   Level: intermediate

406: .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
407: @*/
408: PetscErrorCode TSBasicSymplecticGetType(TS ts, TSBasicSymplecticType *bsymptype)
409: {
410:   PetscFunctionBegin;
412:   PetscUseMethod(ts, "TSBasicSymplecticGetType_C", (TS, TSBasicSymplecticType *), (ts, bsymptype));
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts, TSBasicSymplecticType bsymptype)
417: {
418:   TS_BasicSymplectic       *bsymp = (TS_BasicSymplectic *)ts->data;
419:   BasicSymplecticSchemeLink link;
420:   PetscBool                 match;

422:   PetscFunctionBegin;
423:   if (bsymp->scheme) {
424:     PetscCall(PetscStrcmp(bsymp->scheme->name, bsymptype, &match));
425:     if (match) PetscFunctionReturn(PETSC_SUCCESS);
426:   }
427:   for (link = BasicSymplecticSchemeList; link; link = link->next) {
428:     PetscCall(PetscStrcmp(link->sch.name, bsymptype, &match));
429:     if (match) {
430:       bsymp->scheme = &link->sch;
431:       PetscFunctionReturn(PETSC_SUCCESS);
432:     }
433:   }
434:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", bsymptype);
435: }

437: static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts, TSBasicSymplecticType *bsymptype)
438: {
439:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

441:   PetscFunctionBegin;
442:   *bsymptype = bsymp->scheme->name;
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /*MC
447:   TSBASICSYMPLECTIC - ODE solver using basic symplectic integration schemes <https://en.wikipedia.org/wiki/Symplectic_integrator>

449:   These methods are intended for separable Hamiltonian systems

451:   $$
452:   \begin{align*}
453:   qdot = dH(q,p,t)/dp   \\
454:   pdot = -dH(q,p,t)/dq
455:   \end{align*}
456:   $$

458:   where the Hamiltonian can be split into the sum of kinetic energy and potential energy

460:   $$
461:   H(q,p,t) = T(p,t) + V(q,t).
462:   $$

464:   As a result, the system can be generally represented by

466:   $$
467:   \begin{align*}
468:   qdot = f(p,t) = dT(p,t)/dp \\
469:   pdot = g(q,t) = -dV(q,t)/dq
470:   \end{align*}
471:   $$

473:   and solved iteratively with

475:   $$
476:   \begin{align*}
477:   q_new = q_old + d_i*h*f(p_old,t_old) \\
478:   t_new = t_old + d_i*h \\
479:   p_new = p_old + c_i*h*g(p_new,t_new) \\
480:   i     = 0,1,...,n.
481:   \end{align*}
482:   $$

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

488:   Level: beginner

490: .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TSRHSSplitSetIS()`, `TSRHSSplitSetRHSFunction()`, `TSType`
491: M*/
492: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
493: {
494:   TS_BasicSymplectic *bsymp;

496:   PetscFunctionBegin;
497:   PetscCall(TSBasicSymplecticInitializePackage());
498:   PetscCall(PetscNew(&bsymp));
499:   ts->data = (void *)bsymp;

501:   ts->ops->setup           = TSSetUp_BasicSymplectic;
502:   ts->ops->step            = TSStep_BasicSymplectic;
503:   ts->ops->reset           = TSReset_BasicSymplectic;
504:   ts->ops->destroy         = TSDestroy_BasicSymplectic;
505:   ts->ops->setfromoptions  = TSSetFromOptions_BasicSymplectic;
506:   ts->ops->view            = TSView_BasicSymplectic;
507:   ts->ops->interpolate     = TSInterpolate_BasicSymplectic;
508:   ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;

510:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic));
511:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic));

513:   PetscCall(TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault));
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }