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, next_time_step = ts->time_step;
    212:   PetscInt              n;

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

    247: finally:
    248:   if (!stageok) ts->reason = TS_DIVERGED_STEP_REJECTED;
    249:   else ts->ptime += next_time_step;
    250:   PetscCall(VecRestoreSubVector(update, is_q, &q_update));
    251:   PetscCall(VecRestoreSubVector(update, is_p, &p_update));
    252:   PetscFunctionReturn(PETSC_SUCCESS);
    253: }

    255: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, void *ctx)
    256: {
    257:   PetscFunctionBegin;
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

    261: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
    262: {
    263:   PetscFunctionBegin;
    264:   PetscFunctionReturn(PETSC_SUCCESS);
    265: }

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

    273: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
    274: {
    275:   PetscFunctionBegin;
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
    280: {
    281:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
    282:   DM                  dm;

    284:   PetscFunctionBegin;
    285:   PetscCall(TSRHSSplitGetIS(ts, "position", &bsymp->is_q));
    286:   PetscCall(TSRHSSplitGetIS(ts, "momentum", &bsymp->is_p));
    287:   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");
    288:   PetscCall(TSRHSSplitGetSubTS(ts, "position", &bsymp->subts_q));
    289:   PetscCall(TSRHSSplitGetSubTS(ts, "momentum", &bsymp->subts_p));
    290:   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");

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

    294:   PetscCall(TSGetAdapt(ts, &ts->adapt));
    295:   PetscCall(TSAdaptCandidatesClear(ts->adapt)); /* make sure to use fixed time stepping */
    296:   PetscCall(TSGetDM(ts, &dm));
    297:   if (dm) {
    298:     PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_BasicSymplectic, DMRestrictHook_BasicSymplectic, ts));
    299:     PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_BasicSymplectic, DMSubDomainRestrictHook_BasicSymplectic, ts));
    300:   }
    301:   PetscFunctionReturn(PETSC_SUCCESS);
    302: }

    304: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
    305: {
    306:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

    308:   PetscFunctionBegin;
    309:   PetscCall(VecDestroy(&bsymp->update));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
    314: {
    315:   PetscFunctionBegin;
    316:   PetscCall(TSReset_BasicSymplectic(ts));
    317:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", NULL));
    318:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", NULL));
    319:   PetscCall(PetscFree(ts->data));
    320:   PetscFunctionReturn(PETSC_SUCCESS);
    321: }

    323: static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems PetscOptionsObject)
    324: {
    325:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;

    327:   PetscFunctionBegin;
    328:   PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
    329:   {
    330:     BasicSymplecticSchemeLink link;
    331:     PetscInt                  count, choice;
    332:     PetscBool                 flg;
    333:     const char              **namelist;

    335:     for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++);
    336:     PetscCall(PetscMalloc1(count, (char ***)&namelist));
    337:     for (link = BasicSymplecticSchemeList, count = 0; link; link = link->next, count++) namelist[count] = link->sch.name;
    338:     PetscCall(PetscOptionsEList("-ts_basicsymplectic_type", "Family of basic symplectic integration method", "TSBasicSymplecticSetType", (const char *const *)namelist, count, bsymp->scheme->name, &choice, &flg));
    339:     if (flg) PetscCall(TSBasicSymplecticSetType(ts, namelist[choice]));
    340:     PetscCall(PetscFree(namelist));
    341:   }
    342:   PetscOptionsHeadEnd();
    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: /*@
    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: /*@
    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:   \dot q &= \frac{dH(q,p,t)}{dp}   \\
    454:   \dot p &= -\frac{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:   \dot q &= f(p,t) = \frac{dT(p,t)}{dp} \\
    469:   \dot p &= g(q,t) = -\frac{dV(q,t)}{dq}
    470:   \end{align*}
    471:   $$

    473:   and solved iteratively with $i \in [0, n]$

    475:   $$
    476:   \begin{align*}
    477:   q_{new} &= q_{old} + h d_i f(p_{old}, t_{old}) \\
    478:   t_{new} &= t_{old} + h d_i \\
    479:   p_{new} &= p_{old} + h c_i g(q_{new}, t_{new})
    480:   \end{align*}
    481:   $$

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

    487:   Level: beginner

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

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

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

    508:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic));
    509:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic));

    511:   PetscCall(TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault));
    512:   PetscFunctionReturn(PETSC_SUCCESS);
    513: }