1 /*
2 Code for Timestepping with basic symplectic integrators for separable Hamiltonian systems
3 */
4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
5 #include <petscdm.h>
6
7 static TSBasicSymplecticType TSBasicSymplecticDefault = TSBASICSYMPLECTICSIEULER;
8 static PetscBool TSBasicSymplecticRegisterAllCalled;
9 static PetscBool TSBasicSymplecticPackageInitialized;
10
11 typedef struct _BasicSymplecticScheme *BasicSymplecticScheme;
12 typedef struct _BasicSymplecticSchemeLink *BasicSymplecticSchemeLink;
13
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;
31
32 /*MC
33 TSBASICSYMPLECTICSIEULER - first order semi-implicit Euler method
34
35 Level: intermediate
36
37 .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
38 M*/
39
40 /*MC
41 TSBASICSYMPLECTICVELVERLET - second order Velocity Verlet method (leapfrog method with starting process and determining velocity and position at the same time)
42
43 Level: intermediate
44
45 .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
46 M*/
47
48 /*@C
49 TSBasicSymplecticRegisterAll - Registers all of the basic symplectic integration methods in `TSBASICSYMPLECTIC`
50
51 Not Collective, but should be called by all processes which will need the schemes to be registered
52
53 Level: advanced
54
55 .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticRegisterDestroy()`
56 @*/
TSBasicSymplecticRegisterAll(void)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 }
81
82 /*@C
83 TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by `TSBasicSymplecticRegister()`.
84
85 Not Collective
86
87 Level: advanced
88
89 .seealso: [](ch_ts), `TSBasicSymplecticRegister()`, `TSBasicSymplecticRegisterAll()`, `TSBASICSYMPLECTIC`
90 @*/
TSBasicSymplecticRegisterDestroy(void)91 PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
92 {
93 BasicSymplecticSchemeLink link;
94
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 }
106
107 /*@C
108 TSBasicSymplecticInitializePackage - This function initializes everything in the `TSBASICSYMPLECTIC` package. It is called
109 from `TSInitializePackage()`.
110
111 Level: developer
112
113 .seealso: [](ch_ts), `PetscInitialize()`, `TSBASICSYMPLECTIC`
114 @*/
TSBasicSymplecticInitializePackage(void)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 }
124
125 /*@C
126 TSBasicSymplecticFinalizePackage - This function destroys everything in the `TSBASICSYMPLECTIC` package. It is
127 called from `PetscFinalize()`.
128
129 Level: developer
130
131 .seealso: [](ch_ts), `PetscFinalize()`, `TSBASICSYMPLECTIC`
132 @*/
TSBasicSymplecticFinalizePackage(void)133 PetscErrorCode TSBasicSymplecticFinalizePackage(void)
134 {
135 PetscFunctionBegin;
136 TSBasicSymplecticPackageInitialized = PETSC_FALSE;
137 PetscCall(TSBasicSymplecticRegisterDestroy());
138 PetscFunctionReturn(PETSC_SUCCESS);
139 }
140
141 /*@C
142 TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.
143
144 Not Collective, but the same schemes should be registered on all processes on which they will be used
145
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)
152
153 Level: advanced
154
155 Note:
156 Several symplectic methods are provided, this function is only needed to create new methods.
157
158 .seealso: [](ch_ts), `TSBASICSYMPLECTIC`
159 @*/
TSBasicSymplecticRegister(TSRosWType name,PetscInt order,PetscInt s,PetscReal c[],PetscReal d[])160 PetscErrorCode TSBasicSymplecticRegister(TSRosWType name, PetscInt order, PetscInt s, PetscReal c[], PetscReal d[])
161 {
162 BasicSymplecticSchemeLink link;
163 BasicSymplecticScheme scheme;
164
165 PetscFunctionBegin;
166 PetscAssertPointer(name, 1);
167 PetscAssertPointer(c, 4);
168 PetscAssertPointer(d, 5);
169
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 }
183
184 /*
185 The simplified form of the equations are:
186
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
191
192 Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.
193
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
201
202 */
TSStep_BasicSymplectic(TS ts)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;
213
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 }
246
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 }
254
DMCoarsenHook_BasicSymplectic(DM fine,DM coarse,PetscCtx ctx)255 static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine, DM coarse, PetscCtx ctx)
256 {
257 PetscFunctionBegin;
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260
DMRestrictHook_BasicSymplectic(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)261 static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
262 {
263 PetscFunctionBegin;
264 PetscFunctionReturn(PETSC_SUCCESS);
265 }
266
DMSubDomainHook_BasicSymplectic(DM dm,DM subdm,PetscCtx ctx)267 static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm, DM subdm, PetscCtx ctx)
268 {
269 PetscFunctionBegin;
270 PetscFunctionReturn(PETSC_SUCCESS);
271 }
272
DMSubDomainRestrictHook_BasicSymplectic(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,PetscCtx ctx)273 static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
274 {
275 PetscFunctionBegin;
276 PetscFunctionReturn(PETSC_SUCCESS);
277 }
278
TSSetUp_BasicSymplectic(TS ts)279 static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
280 {
281 TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
282 DM dm;
283
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");
291
292 PetscCall(VecDuplicate(ts->vec_sol, &bsymp->update));
293
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 }
303
TSReset_BasicSymplectic(TS ts)304 static PetscErrorCode TSReset_BasicSymplectic(TS ts)
305 {
306 TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
307
308 PetscFunctionBegin;
309 PetscCall(VecDestroy(&bsymp->update));
310 PetscFunctionReturn(PETSC_SUCCESS);
311 }
312
TSDestroy_BasicSymplectic(TS ts)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 }
322
TSSetFromOptions_BasicSymplectic(TS ts,PetscOptionItems PetscOptionsObject)323 static PetscErrorCode TSSetFromOptions_BasicSymplectic(TS ts, PetscOptionItems PetscOptionsObject)
324 {
325 TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
326
327 PetscFunctionBegin;
328 PetscOptionsHeadBegin(PetscOptionsObject, "Basic symplectic integrator options");
329 {
330 BasicSymplecticSchemeLink link;
331 PetscInt count, choice;
332 PetscBool flg;
333 const char **namelist;
334
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 }
345
TSInterpolate_BasicSymplectic(TS ts,PetscReal t,Vec X)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;
351
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 }
357
TSComputeLinearStability_BasicSymplectic(TS ts,PetscReal xr,PetscReal xi,PetscReal * yr,PetscReal * yi)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 }
365
366 /*@
367 TSBasicSymplecticSetType - Set the type of the basic symplectic method
368
369 Logically Collective
370
371 Input Parameters:
372 + ts - timestepping context
373 - bsymptype - type of the symplectic scheme
374
375 Options Database Key:
376 . -ts_basicsymplectic_type <scheme> - select the scheme
377
378 Level: intermediate
379
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.
384
385 .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`
386 @*/
TSBasicSymplecticSetType(TS ts,TSBasicSymplecticType bsymptype)387 PetscErrorCode TSBasicSymplecticSetType(TS ts, TSBasicSymplecticType bsymptype)
388 {
389 PetscFunctionBegin;
390 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
391 PetscTryMethod(ts, "TSBasicSymplecticSetType_C", (TS, TSBasicSymplecticType), (ts, bsymptype));
392 PetscFunctionReturn(PETSC_SUCCESS);
393 }
394
395 /*@
396 TSBasicSymplecticGetType - Get the type of the basic symplectic method
397
398 Logically Collective
399
400 Input Parameters:
401 + ts - timestepping context
402 - bsymptype - type of the basic symplectic scheme
403
404 Level: intermediate
405
406 .seealso: [](ch_ts), `TSBASICSYMPLECTIC`, `TSBasicSymplecticType`, `TSBasicSymplecticSetType()`
407 @*/
TSBasicSymplecticGetType(TS ts,TSBasicSymplecticType * bsymptype)408 PetscErrorCode TSBasicSymplecticGetType(TS ts, TSBasicSymplecticType *bsymptype)
409 {
410 PetscFunctionBegin;
411 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
412 PetscUseMethod(ts, "TSBasicSymplecticGetType_C", (TS, TSBasicSymplecticType *), (ts, bsymptype));
413 PetscFunctionReturn(PETSC_SUCCESS);
414 }
415
TSBasicSymplecticSetType_BasicSymplectic(TS ts,TSBasicSymplecticType bsymptype)416 static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts, TSBasicSymplecticType bsymptype)
417 {
418 TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
419 BasicSymplecticSchemeLink link;
420 PetscBool match;
421
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 }
436
TSBasicSymplecticGetType_BasicSymplectic(TS ts,TSBasicSymplecticType * bsymptype)437 static PetscErrorCode TSBasicSymplecticGetType_BasicSymplectic(TS ts, TSBasicSymplecticType *bsymptype)
438 {
439 TS_BasicSymplectic *bsymp = (TS_BasicSymplectic *)ts->data;
440
441 PetscFunctionBegin;
442 *bsymptype = bsymp->scheme->name;
443 PetscFunctionReturn(PETSC_SUCCESS);
444 }
445
446 /*MC
447 TSBASICSYMPLECTIC - ODE solver using basic symplectic integration schemes <https://en.wikipedia.org/wiki/Symplectic_integrator>
448
449 These methods are intended for separable Hamiltonian systems
450
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 $$
457
458 where the Hamiltonian can be split into the sum of kinetic energy and potential energy
459
460 $$
461 H(q,p,t) = T(p,t) + V(q,t).
462 $$
463
464 As a result, the system can be generally represented by
465
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 $$
472
473 and solved iteratively with $i \in [0, n]$
474
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 $$
482
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.
486
487 Level: beginner
488
489 .seealso: [](ch_ts), `TSCreate()`, `TSSetType()`, `TSRHSSplitSetIS()`, `TSRHSSplitSetRHSFunction()`, `TSType`
490 M*/
TSCreate_BasicSymplectic(TS ts)491 PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
492 {
493 TS_BasicSymplectic *bsymp;
494
495 PetscFunctionBegin;
496 PetscCall(TSBasicSymplecticInitializePackage());
497 PetscCall(PetscNew(&bsymp));
498 ts->data = (void *)bsymp;
499
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;
507
508 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticSetType_C", TSBasicSymplecticSetType_BasicSymplectic));
509 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBasicSymplecticGetType_C", TSBasicSymplecticGetType_BasicSymplectic));
510
511 PetscCall(TSBasicSymplecticSetType(ts, TSBasicSymplecticDefault));
512 PetscFunctionReturn(PETSC_SUCCESS);
513 }
514