1 /*
2 Code for timestepping with additive Runge-Kutta IMEX method or Diagonally Implicit Runge-Kutta methods.
3
4 Notes:
5 For ARK, the general system is written as
6
7 F(t,U,Udot) = G(t,U)
8
9 where F represents the stiff part of the physics and G represents the non-stiff part.
10
11 */
12 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
13 #include <petscdm.h>
14 #include <../src/ts/impls/arkimex/arkimex.h>
15 #include <../src/ts/impls/arkimex/fsarkimex.h>
16
17 static ARKTableauLink ARKTableauList;
18 static TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3;
19 static TSDIRKType TSDIRKDefault = TSDIRKES213SAL;
20 static PetscBool TSARKIMEXRegisterAllCalled;
21 static PetscBool TSARKIMEXPackageInitialized;
22 static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec);
23
24 /*MC
25 TSARKIMEXARS122 - Second order ARK IMEX scheme, {cite}`ascher_1997`
26
27 This method has one explicit stage and one implicit stage.
28
29 Options Database Key:
30 . -ts_arkimex_type ars122 - set arkimex type to ars122
31
32 Level: advanced
33
34 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
35 M*/
36
37 /*MC
38 TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part.
39
40 This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu.
41
42 Options Database Key:
43 . -ts_arkimex_type a2 - set arkimex type to a2
44
45 Level: advanced
46
47 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
48 M*/
49
50 /*MC
51 TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part, {cite}`pareschi_2005`
52
53 This method has two implicit stages, and L-stable implicit scheme.
54
55 Options Database Key:
56 . -ts_arkimex_type l2 - set arkimex type to l2
57
58 Level: advanced
59
60 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
61 M*/
62
63 /*MC
64 TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method.
65
66 This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used.
67
68 Options Database Key:
69 . -ts_arkimex_type 1bee - set arkimex type to 1bee
70
71 Level: advanced
72
73 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
74 M*/
75
76 /*MC
77 TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part.
78
79 This method has one explicit stage and two implicit stages. The implicit part is the same as in TSARKIMEX2D and TSARKIMEX2E, but the explicit part has a larger stability region on the negative real axis. This method was provided by Emil Constantinescu.
80
81 Options Database Key:
82 . -ts_arkimex_type 2c - set arkimex type to 2c
83
84 Level: advanced
85
86 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
87 M*/
88
89 /*MC
90 TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part.
91
92 This method has one explicit stage and two implicit stages. The stability function is independent of the explicit part in the infinity limit of the implicit component. This method was provided by Emil Constantinescu.
93
94 Options Database Key:
95 . -ts_arkimex_type 2d - set arkimex type to 2d
96
97 Level: advanced
98
99 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
100 M*/
101
102 /*MC
103 TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part.
104
105 This method has one explicit stage and two implicit stages. It is an optimal method developed by Emil Constantinescu.
106
107 Options Database Key:
108 . -ts_arkimex_type 2e - set arkimex type to 2e
109
110 Level: advanced
111
112 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
113 M*/
114
115 /*MC
116 TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme, {cite}`pareschi_2005`
117
118 This method has three implicit stages.
119
120 This method is referred to as SSP2-(3,3,2) in <https://arxiv.org/abs/1110.4375>
121
122 Options Database Key:
123 . -ts_arkimex_type prssp2 - set arkimex type to prssp2
124
125 Level: advanced
126
127 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
128 M*/
129
130 /*MC
131 TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`
132
133 This method has one explicit stage and three implicit stages.
134
135 Options Database Key:
136 . -ts_arkimex_type 3 - set arkimex type to 3
137
138 Level: advanced
139
140 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
141 M*/
142
143 /*MC
144 TSARKIMEXARS443 - Third order ARK IMEX scheme, {cite}`ascher_1997`
145
146 This method has one explicit stage and four implicit stages.
147
148 Options Database Key:
149 . -ts_arkimex_type ars443 - set arkimex type to ars443
150
151 Level: advanced
152
153 Notes:
154 This method is referred to as ARS(4,4,3) in <https://arxiv.org/abs/1110.4375>
155
156 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
157 M*/
158
159 /*MC
160 TSARKIMEXBPR3 - Third order ARK IMEX scheme. Referred to as ARK3 in <https://arxiv.org/abs/1110.4375>
161
162 This method has one explicit stage and four implicit stages.
163
164 Options Database Key:
165 . -ts_arkimex_type bpr3 - set arkimex type to bpr3
166
167 Level: advanced
168
169 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
170 M*/
171
172 /*MC
173 TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`.
174
175 This method has one explicit stage and four implicit stages.
176
177 Options Database Key:
178 . -ts_arkimex_type 4 - set arkimex type to4
179
180 Level: advanced
181
182 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
183 M*/
184
185 /*MC
186 TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`.
187
188 This method has one explicit stage and five implicit stages.
189
190 Options Database Key:
191 . -ts_arkimex_type 5 - set arkimex type to 5
192
193 Level: advanced
194
195 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
196 M*/
197
198 /*MC
199 TSDIRKS212 - Second order DIRK scheme.
200
201 This method has two implicit stages with an embedded method of other 1.
202 See `TSDIRK` for additional details.
203
204 Options Database Key:
205 . -ts_dirk_type s212 - select this method.
206
207 Level: advanced
208
209 Note:
210 This is the default DIRK scheme in SUNDIALS.
211
212 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
213 M*/
214
215 /*MC
216 TSDIRKES122SAL - First order DIRK scheme <https://arxiv.org/abs/1803.01613>
217
218 Uses backward Euler as advancing method and trapezoidal rule as embedded method. See `TSDIRK` for additional details.
219
220 Options Database Key:
221 . -ts_dirk_type es122sal - select this method.
222
223 Level: advanced
224
225 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
226 M*/
227
228 /*MC
229 TSDIRKES213SAL - Second order DIRK scheme {cite}`kennedy2019diagonally`. Also known as TR-BDF2, see{cite}`hosea1996analysis`
230
231 See `TSDIRK` for additional details.
232
233 Options Database Key:
234 . -ts_dirk_type es213sal - select this method.
235
236 Level: advanced
237
238 Note:
239 This is the default DIRK scheme used in PETSc.
240
241 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
242 M*/
243
244 /*MC
245 TSDIRKES324SAL - Third order DIRK scheme, {cite}`kennedy2019diagonally`
246
247 See `TSDIRK` for additional details.
248
249 Options Database Key:
250 . -ts_dirk_type es324sal - select this method.
251
252 Level: advanced
253
254 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
255 M*/
256
257 /*MC
258 TSDIRKES325SAL - Third order DIRK scheme {cite}`kennedy2019diagonally`.
259
260 See `TSDIRK` for additional details.
261
262 Options Database Key:
263 . -ts_dirk_type es325sal - select this method.
264
265 Level: advanced
266
267 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
268 M*/
269
270 /*MC
271 TSDIRK657A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
272
273 See `TSDIRK` for additional details.
274
275 Options Database Key:
276 . -ts_dirk_type 657a - select this method.
277
278 Level: advanced
279
280 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
281 M*/
282
283 /*MC
284 TSDIRKES648SA - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
285
286 See `TSDIRK` for additional details.
287
288 Options Database Key:
289 . -ts_dirk_type es648sa - select this method.
290
291 Level: advanced
292
293 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
294 M*/
295
296 /*MC
297 TSDIRK658A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
298
299 See `TSDIRK` for additional details.
300
301 Options Database Key:
302 . -ts_dirk_type 658a - select this method.
303
304 Level: advanced
305
306 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
307 M*/
308
309 /*MC
310 TSDIRKS659A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
311
312 See `TSDIRK` for additional details.
313
314 Options Database Key:
315 . -ts_dirk_type s659a - select this method.
316
317 Level: advanced
318
319 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
320 M*/
321
322 /*MC
323 TSDIRK7510SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
324
325 See `TSDIRK` for additional details.
326
327 Options Database Key:
328 . -ts_dirk_type 7510sal - select this method.
329
330 Level: advanced
331
332 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
333 M*/
334
335 /*MC
336 TSDIRKES7510SA - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
337
338 See `TSDIRK` for additional details.
339
340 Options Database Key:
341 . -ts_dirk_type es7510sa - select this method.
342
343 Level: advanced
344
345 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
346 M*/
347
348 /*MC
349 TSDIRK759A - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
350
351 See `TSDIRK` for additional details.
352
353 Options Database Key:
354 . -ts_dirk_type 759a - select this method.
355
356 Level: advanced
357
358 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
359 M*/
360
361 /*MC
362 TSDIRKS7511SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
363
364 See `TSDIRK` for additional details.
365
366 Options Database Key:
367 . -ts_dirk_type s7511sal - select this method.
368
369 Level: advanced
370
371 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
372 M*/
373
374 /*MC
375 TSDIRK8614A - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
376
377 See `TSDIRK` for additional details.
378
379 Options Database Key:
380 . -ts_dirk_type 8614a - select this method.
381
382 Level: advanced
383
384 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
385 M*/
386
387 /*MC
388 TSDIRK8616SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
389
390 See `TSDIRK` for additional details.
391
392 Options Database Key:
393 . -ts_dirk_type 8616sal - select this method.
394
395 Level: advanced
396
397 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
398 M*/
399
400 /*MC
401 TSDIRKES8516SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
402
403 See `TSDIRK` for additional details.
404
405 Options Database Key:
406 . -ts_dirk_type es8516sal - select this method.
407
408 Level: advanced
409
410 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
411 M*/
412
TSHasRHSFunction(TS ts,PetscBool * has)413 PetscErrorCode TSHasRHSFunction(TS ts, PetscBool *has)
414 {
415 TSRHSFunctionFn *func;
416
417 PetscFunctionBegin;
418 *has = PETSC_FALSE;
419 PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL));
420 if (func) *has = PETSC_TRUE;
421 PetscFunctionReturn(PETSC_SUCCESS);
422 }
423
424 /*@C
425 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX`
426
427 Not Collective, but should be called by all processes which will need the schemes to be registered
428
429 Level: advanced
430
431 .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()`
432 @*/
TSARKIMEXRegisterAll(void)433 PetscErrorCode TSARKIMEXRegisterAll(void)
434 {
435 PetscFunctionBegin;
436 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
437 TSARKIMEXRegisterAllCalled = PETSC_TRUE;
438
439 #define RC PetscRealConstant
440 #define s2 RC(1.414213562373095048802) /* PetscSqrtReal((PetscReal)2.0) */
441 #define us2 RC(0.2928932188134524755992) /* 1.0-1.0/PetscSqrtReal((PetscReal)2.0); */
442
443 /* Diagonally implicit methods */
444 {
445 /* DIRK212, default of SUNDIALS */
446 const PetscReal A[2][2] = {
447 {RC(1.0), RC(0.0)},
448 {RC(-1.0), RC(1.0)}
449 };
450 const PetscReal b[2] = {RC(0.5), RC(0.5)};
451 const PetscReal bembed[2] = {RC(1.0), RC(0.0)};
452 PetscCall(TSDIRKRegister(TSDIRKS212, 2, 2, &A[0][0], b, NULL, bembed, 1, b));
453 }
454
455 {
456 /* ESDIRK12 from https://arxiv.org/pdf/1803.01613.pdf */
457 const PetscReal A[2][2] = {
458 {RC(0.0), RC(0.0)},
459 {RC(0.0), RC(1.0)}
460 };
461 const PetscReal b[2] = {RC(0.0), RC(1.0)};
462 const PetscReal bembed[2] = {RC(0.5), RC(0.5)};
463 PetscCall(TSDIRKRegister(TSDIRKES122SAL, 1, 2, &A[0][0], b, NULL, bembed, 1, b));
464 }
465
466 {
467 /* ESDIRK213L[2]SA from KC16.
468 TR-BDF2 from Hosea and Shampine
469 ESDIRK23 in https://arxiv.org/pdf/1803.01613.pdf */
470 const PetscReal A[3][3] = {
471 {RC(0.0), RC(0.0), RC(0.0)},
472 {us2, us2, RC(0.0)},
473 {s2 / RC(4.0), s2 / RC(4.0), us2 },
474 };
475 const PetscReal b[3] = {s2 / RC(4.0), s2 / RC(4.0), us2};
476 const PetscReal bembed[3] = {(RC(1.0) - s2 / RC(4.0)) / RC(3.0), (RC(3.0) * s2 / RC(4.0) + RC(1.0)) / RC(3.0), us2 / RC(3.0)};
477 PetscCall(TSDIRKRegister(TSDIRKES213SAL, 2, 3, &A[0][0], b, NULL, bembed, 1, b));
478 }
479
480 {
481 #define g RC(0.43586652150845899941601945)
482 #define g2 PetscSqr(g)
483 #define g3 g *g2
484 #define g4 PetscSqr(g2)
485 #define g5 g *g4
486 #define c3 RC(1.0)
487 #define a32 c3 *(c3 - RC(2.0) * g) / (RC(4.0) * g)
488 #define b2 (-RC(2.0) + RC(3.0) * c3 + RC(6.0) * g * (RC(1.0) - c3)) / (RC(12.0) * g * (c3 - RC(2.0) * g))
489 #define b3 (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) / (RC(3.0) * c3 * (c3 - RC(2.0) * g))
490 #if 0
491 /* This is for c3 = 3/5 */
492 #define bh2 \
493 c3 * (-RC(1.0) + RC(6.0) * g - RC(23.0) * g3 + RC(12.0) * g4 - RC(6.0) * g5) / (RC(4.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) + (RC(3.0) - RC(27.0) * g + RC(68.0) * g2 - RC(55.0) * g3 + RC(21.0) * g4 - RC(6.0) * g5) / (RC(2.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2))
494 #define bh3 -g * (-RC(2.0) + RC(21.0) * g - RC(68.0) * g2 + RC(79.0) * g3 - RC(33.0) * g4 + RC(12.0) * g5) / (RC(2.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2))
495 #define bh4 -RC(3.0) * g2 * (-RC(1.0) + RC(4.0) * g - RC(2.0) * g2 + g3) / (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)
496 #else
497 /* This is for c3 = 1.0 */
498 #define bh2 a32
499 #define bh3 g
500 #define bh4 RC(0.0)
501 #endif
502 /* ESDIRK3(2I)4L[2]SA from KC16 with c3 = 1.0 */
503 /* Given by Kvaerno https://link.springer.com/article/10.1023/b:bitn.0000046811.70614.38 */
504 const PetscReal A[4][4] = {
505 {RC(0.0), RC(0.0), RC(0.0), RC(0.0)},
506 {g, g, RC(0.0), RC(0.0)},
507 {c3 - a32 - g, a32, g, RC(0.0)},
508 {RC(1.0) - b2 - b3 - g, b2, b3, g },
509 };
510 const PetscReal b[4] = {RC(1.0) - b2 - b3 - g, b2, b3, g};
511 const PetscReal bembed[4] = {RC(1.0) - bh2 - bh3 - bh4, bh2, bh3, bh4};
512 PetscCall(TSDIRKRegister(TSDIRKES324SAL, 3, 4, &A[0][0], b, NULL, bembed, 1, b));
513 #undef g
514 #undef g2
515 #undef g3
516 #undef c3
517 #undef a32
518 #undef b2
519 #undef b3
520 #undef bh2
521 #undef bh3
522 #undef bh4
523 }
524
525 {
526 /* ESDIRK3(2I)5L[2]SA from KC16 */
527 const PetscReal A[5][5] = {
528 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
529 {RC(9.0) / RC(40.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0), RC(0.0) },
530 {RC(19.0) / RC(72.0), RC(14.0) / RC(45.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0) },
531 {RC(3337.0) / RC(11520.0), RC(233.0) / RC(720.0), RC(207.0) / RC(1280.0), RC(9.0) / RC(40.0), RC(0.0) },
532 {RC(7415.0) / RC(34776.0), RC(9920.0) / RC(30429.0), RC(4845.0) / RC(9016.0), -RC(5827.0) / RC(19320.0), RC(9.0) / RC(40.0)},
533 };
534 const PetscReal b[5] = {RC(7415.0) / RC(34776.0), RC(9920.0) / RC(30429.0), RC(4845.0) / RC(9016.0), -RC(5827.0) / RC(19320.0), RC(9.0) / RC(40.0)};
535 const PetscReal bembed[5] = {RC(23705.0) / RC(104328.0), RC(29720.0) / RC(91287.0), RC(4225.0) / RC(9016.0), -RC(69304987.0) / RC(337732920.0), RC(42843.0) / RC(233080.0)};
536 PetscCall(TSDIRKRegister(TSDIRKES325SAL, 3, 5, &A[0][0], b, NULL, bembed, 1, b));
537 }
538
539 {
540 // DIRK(6,6)[1]A[(7,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
541 const PetscReal A[7][7] = {
542 {RC(0.303487844706747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
543 {RC(-0.279756492709814), RC(0.500032236020747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
544 {RC(0.280583215743895), RC(-0.438560061586751), RC(0.217250734515736), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
545 {RC(-0.0677678738539846), RC(0.984312781232293), RC(-0.266720192540149), RC(0.2476680834526), RC(0.0), RC(0.0), RC(0.0) },
546 {RC(0.125671616147993), RC(-0.995401751002415), RC(0.761333109549059), RC(-0.210281837202208), RC(0.866743712636936), RC(0.0), RC(0.0) },
547 {RC(-0.368056238801488), RC(-0.999928082701516), RC(0.534734253232519), RC(-0.174856916279082), RC(0.615007160285509), RC(0.696549912132029), RC(0.0) },
548 {RC(-0.00570546839653984), RC(-0.113110431835656), RC(-0.000965563207671587), RC(-0.000130490084629567), RC(0.00111737736895673), RC(-0.279385587378871), RC(0.618455906845342)}
549 };
550 const PetscReal b[7] = {RC(0.257561510484877), RC(0.234281287047716), RC(0.126658904241469), RC(0.252363215441784), RC(0.396701083526306), RC(-0.267566000742152), RC(0.0)};
551 const PetscReal bembed[7] = {RC(0.257561510484945), RC(0.387312822934391), RC(0.126658904241468), RC(0.252363215441784), RC(0.396701083526306), RC(-0.267566000742225), RC(-0.153031535886669)};
552 PetscCall(TSDIRKRegister(TSDIRK657A, 6, 7, &A[0][0], b, NULL, bembed, 1, b));
553 }
554 {
555 // ESDIRK(8,6)[2]SA[(8,4)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
556 const PetscReal A[8][8] = {
557 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
558 {RC(0.333222149217725), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
559 {RC(0.0639743773182214), RC(-0.0830330224410214), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
560 {RC(-0.728522201369326), RC(-0.210414479522485), RC(0.532519916559342), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
561 {RC(-0.175135269272067), RC(0.666675582067552), RC(-0.304400907370867), RC(0.656797712445756), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0) },
562 {RC(0.222695802705462), RC(-0.0948971794681061), RC(-0.0234336346686545), RC(-0.45385925012042), RC(0.0283910313826958), RC(0.333222149217725), RC(0.0), RC(0.0) },
563 {RC(-0.132534078051299), RC(0.702597935004879), RC(-0.433316453128078), RC(0.893717488547587), RC(0.057381454791406), RC(-0.207798411552402), RC(0.333222149217725), RC(0.0) },
564 {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)}
565 };
566 const PetscReal b[8] = {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)};
567 const PetscReal bembed[8] = {RC(0.0), RC(0.292331064554014), RC(0.409676102283681), RC(-0.002094718084982), RC(-0.282771520835975), RC(0.113862336644901), RC(0.181973572260693), RC(0.287023163177669)};
568 PetscCall(TSDIRKRegister(TSDIRKES648SA, 6, 8, &A[0][0], b, NULL, bembed, 1, b));
569 }
570 {
571 // DIRK(8,6)[1]SAL[(8,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
572 const PetscReal A[8][8] = {
573 {RC(0.477264457385826), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
574 {RC(-0.197052588415002), RC(0.476363428459584), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
575 {RC(-0.0347674430372966), RC(0.633051807335483), RC(0.193634310075028), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
576 {RC(0.0967797668578702), RC(-0.193533526466535), RC(-0.000207622945800473), RC(0.159572204849431), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
577 {RC(0.162527231819875), RC(-0.249672513547382), RC(-0.0459079972041795), RC(0.36579476400859), RC(0.255752838307699), RC(0.0), RC(0.0), RC(0.0) },
578 {RC(-0.00707603197171262), RC(0.846299854860295), RC(0.344020016925018), RC(-0.0720926054548865), RC(-0.215492331980875), RC(0.104341097622161), RC(0.0), RC(0.0) },
579 {RC(0.00176857935179744), RC(0.0779960013127515), RC(0.303333277564557), RC(0.213160806732836), RC(0.351769320319038), RC(-0.381545894386538), RC(0.433517909105558), RC(0.0) },
580 {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)}
581 };
582 const PetscReal b[8] = {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)};
583 const PetscReal bembed[8] = {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.103483943222765), RC(-0.0103721771642262), RC(0.177302191576001)};
584 PetscCall(TSDIRKRegister(TSDIRK658A, 6, 8, &A[0][0], b, NULL, bembed, 1, b));
585 }
586 {
587 // SDIRK(9,6)[1]SAL[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
588 const PetscReal A[9][9] = {
589 {RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
590 {RC(-0.0903514856119419), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
591 {RC(0.172952039138937), RC(-0.35365501036282), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
592 {RC(0.511999875919193), RC(0.0289640332201925), RC(-0.0144030945657094), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
593 {RC(0.00465303495506782), RC(-0.075635818766597), RC(0.217273030786712), RC(-0.0206519428725472), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
594 {RC(0.896145501762472), RC(0.139267327700498), RC(-0.186920979752805), RC(0.0672971012371723), RC(-0.350891963442176), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0) },
595 {RC(0.552959701885751), RC(-0.439360579793662), RC(0.333704002325091), RC(-0.0339426520778416), RC(-0.151947445912595), RC(0.0213825661026943), RC(0.218127781944908), RC(0.0), RC(0.0) },
596 {RC(0.631360374036476), RC(0.724733619641466), RC(-0.432170625425258), RC(0.598611382182477), RC(-0.709087197034345), RC(-0.483986685696934), RC(0.378391562905131), RC(0.218127781944908), RC(0.0) },
597 {RC(0.0), RC(-0.15504452530869), RC(0.194518478660789), RC(0.63515640279203), RC(0.81172278664173), RC(0.110736108691585), RC(-0.495304692414479), RC(-0.319912341007872), RC(0.218127781944908)}
598 };
599 const PetscReal b[9] = {RC(0.0), RC(-0.15504452530869), RC(0.194518478660789), RC(0.63515640279203), RC(0.81172278664173), RC(0.110736108691585), RC(-0.495304692414479), RC(-0.319912341007872), RC(0.218127781944908)};
600 const PetscReal bembed[9] = {RC(3.62671059311602e-16), RC(0.0736615558278942), RC(0.103527397262229), RC(1.00247481935499), RC(0.361377289250057), RC(-0.785425929961365), RC(-0.0170499047960784), RC(0.296321252214769), RC(-0.0348864791524953)};
601 PetscCall(TSDIRKRegister(TSDIRKS659A, 6, 9, &A[0][0], b, NULL, bembed, 1, b));
602 }
603 {
604 // DIRK(10,7)[1]SAL[(10,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
605 const PetscReal A[10][10] = {
606 {RC(0.233704632125264), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
607 {RC(-0.0739324813149407), RC(0.200056838146104), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
608 {RC(0.0943790344044812), RC(0.264056067701605), RC(0.133245202456465), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
609 {RC(0.269084810601201), RC(-0.503479002548384), RC(-0.00486736469695022), RC(0.251518716213569), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
610 {RC(0.145665801918423), RC(0.204983170463176), RC(0.407154634069484), RC(-0.0121039135200389), RC(0.190243622486334), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
611 {RC(0.985450198547345), RC(0.806942652811456), RC(-0.808130934167263), RC(-0.669035819439391), RC(0.0269384406756128), RC(0.462144080607327), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
612 {RC(0.163902957809563), RC(0.228315094960095), RC(0.0745971021260249), RC(0.000509793400156559), RC(0.0166533681378294), RC(-0.0229383879045797), RC(0.103505486637336), RC(0.0), RC(0.0), RC(0.0) },
613 {RC(-0.162694156858437), RC(0.0453478837428434), RC(0.997443481211424), RC(0.200251514941093), RC(-0.000161755458839048), RC(-0.0848134335980281), RC(-0.36438666566666), RC(0.158604420136055), RC(0.0), RC(0.0) },
614 {RC(0.200733156477425), RC(0.239686443444433), RC(0.303837014418929), RC(-0.0534390596279896), RC(0.0314067599640569), RC(-0.00764032790448536), RC(0.0609191260198661), RC(-0.0736319201590642), RC(0.204602530607021), RC(0.0) },
615 {RC(0.0), RC(0.235563761744267), RC(0.658651488684319), RC(0.0308877804992098), RC(-0.906514945595336), RC(-0.0248488551739974), RC(-0.309967582365257), RC(0.191663316925525), RC(0.923933712199542), RC(0.200631323081727)}
616 };
617 const PetscReal b[10] = {RC(0.0), RC(0.235563761744267), RC(0.658651488684319), RC(0.0308877804992098), RC(-0.906514945595336), RC(-0.0248488551739974), RC(-0.309967582365257), RC(0.191663316925525), RC(0.923933712199542), RC(0.200631323081727)};
618 const PetscReal bembed[10] =
619 {RC(0.0), RC(0.222929376486581), RC(0.950668440138169), RC(0.0342694607044032), RC(0.362875840545746), RC(0.223572979288581), RC(-0.764361723526727), RC(0.563476909230026), RC(-0.690896961894185), RC(0.0974656790270323)};
620 PetscCall(TSDIRKRegister(TSDIRK7510SAL, 7, 10, &A[0][0], b, NULL, bembed, 1, b));
621 }
622 {
623 // ESDIRK(10,7)[2]SA[(10,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
624 const PetscReal A[10][10] = {
625 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
626 {RC(0.210055790203419), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
627 {RC(0.255781739921086), RC(0.239850916980976), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
628 {RC(0.286789624880437), RC(0.230494748834778), RC(0.263925149885491), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
629 {RC(-0.0219118128774335), RC(0.897684380345907), RC(-0.657954605498907), RC(0.124962304722633), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
630 {RC(-0.065614879584776), RC(-0.0565630711859497), RC(0.0254881105065311), RC(-0.00368981790650006), RC(-0.0115178258446329), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
631 {RC(0.399860851232098), RC(0.915588469718705), RC(-0.0758429094934412), RC(-0.263369154872759), RC(0.719687583564526), RC(-0.787410407015369), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0) },
632 {RC(0.51693616104628), RC(1.00000540846973), RC(-0.0485110663289207), RC(-0.315208041581942), RC(0.749742806451587), RC(-0.990975090921248), RC(0.0159279583407308), RC(0.210055790203419), RC(0.0), RC(0.0) },
633 {RC(-0.0303062129076945), RC(-0.297035174659034), RC(0.184724697462164), RC(-0.0351876079516183), RC(-0.00324668230690761), RC(0.216151004053531), RC(-0.126676252098317), RC(0.114040254365262), RC(0.210055790203419), RC(0.0) },
634 {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)}
635 };
636 const PetscReal b[10] = {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742),
637 RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)};
638 const PetscReal bembed[10] = {RC(-0.015494246543626), RC(0.167657963820093), RC(0.269858958144236), RC(-0.0443258997755156), RC(0.150049236875266),
639 RC(0.259452082755846), RC(0.244624573502521), RC(-0.215528446920284), RC(0.0487601760292619), RC(0.134945602112201)};
640 PetscCall(TSDIRKRegister(TSDIRKES7510SA, 7, 10, &A[0][0], b, NULL, bembed, 1, b));
641 }
642 {
643 // DIRK(9,7)[1]A[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
644 const PetscReal A[9][9] = {
645 {RC(0.179877789855839), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
646 {RC(-0.100405844885157), RC(0.214948590644819), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
647 {RC(0.112251360198995), RC(-0.206162139150298), RC(0.125159642941958), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
648 {RC(-0.0335164000768257), RC(0.999942349946143), RC(-0.491470853833294), RC(0.19820086325566), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
649 {RC(-0.0417345265478321), RC(0.187864510308215), RC(0.0533789224305102), RC(-0.00822060284862916), RC(0.127670843671646), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
650 {RC(-0.0278257925239257), RC(0.600979340683382), RC(-0.242632273241134), RC(-0.11318753652081), RC(0.164326917632931), RC(0.284116597781395), RC(0.0), RC(0.0), RC(0.0) },
651 {RC(0.041465583858922), RC(0.429657872601836), RC(-0.381323410582524), RC(0.391934277498434), RC(-0.245918275501241), RC(-0.35960669741231), RC(0.184000022289158), RC(0.0), RC(0.0) },
652 {RC(-0.105565651574538), RC(-0.0557833155018609), RC(0.358967568942643), RC(-0.13489263413921), RC(0.129553247260677), RC(0.0992493795371489), RC(-0.15716610563461), RC(0.17918862279814), RC(0.0) },
653 {RC(0.00439696079965225), RC(0.960250486570491), RC(0.143558372286706), RC(0.0819015241056593), RC(0.999562318563625), RC(0.325203439314358), RC(-0.679013149331228), RC(-0.990589559837246), RC(0.0773648037639896)}
654 };
655
656 const PetscReal b[9] = {RC(0.0), RC(0.179291520437966), RC(0.115310295273026), RC(-0.857943261453138), RC(0.654911318641998), RC(1.18713633508094), RC(-0.0949482361570542), RC(-0.37661430946407), RC(0.19285633764033)};
657 const PetscReal bembed[9] = {RC(0.0), RC(0.1897135479408), RC(0.127461414808862), RC(-0.835810807663404), RC(0.665114177777166), RC(1.16481046518346), RC(-0.11661858889792), RC(-0.387303251022099), RC(0.192633041873135)};
658 PetscCall(TSDIRKRegister(TSDIRK759A, 7, 9, &A[0][0], b, NULL, bembed, 1, b));
659 }
660 {
661 // SDIRK(11,7)[1]SAL[(11,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
662 const PetscReal A[11][11] = {
663 {RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
664 {RC(-0.082947368165267), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
665 {RC(0.483452690540751), RC(0.0), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
666 {RC(0.771076453481321), RC(-0.22936926341842), RC(0.289733373208823), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
667 {RC(0.0329683054968892), RC(-0.162397421903366), RC(0.000951777538562805), RC(0.0), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
668 {RC(0.265888743485945), RC(0.606743151103931), RC(0.173443800537369), RC(-0.0433968261546912), RC(-0.385211017224481), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
669 {RC(0.220662294551146), RC(-0.0465078507657608), RC(-0.0333111995282464), RC(0.011801580836998), RC(0.169480801030105), RC(-0.0167974432139385), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
670 {RC(0.323099728365267), RC(0.0288371831672575), RC(-0.0543404318773196), RC(0.0137765831431662), RC(0.0516799019060702), RC(-0.0421359763835713), RC(0.181297932037826), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0) },
671 {RC(-0.164226696476538), RC(0.187552004946792), RC(0.0628674420973025), RC(-0.0108886582703428), RC(-0.0117628641717889), RC(0.0432176880867965), RC(-0.0315206836275473), RC(-0.0846007021638797), RC(0.200252661187742), RC(0.0), RC(0.0) },
672 {RC(0.651428598623771), RC(-0.10208078475356), RC(0.198305701801888), RC(-0.0117354096673789), RC(-0.0440385966743686), RC(-0.0358364455795087), RC(-0.0075408087654097), RC(0.160320941654639), RC(0.017940248694499), RC(0.200252661187742), RC(0.0) },
673 {RC(0.0), RC(-0.266259448580236), RC(-0.615982357748271), RC(0.561474126687165), RC(0.266911112787025), RC(0.219775952207137), RC(0.387847665451514), RC(0.612483137773236), RC(0.330027015806089), RC(-0.6965298655714), RC(0.200252661187742)}
674 };
675 const PetscReal b[11] =
676 {RC(0.0), RC(-0.266259448580236), RC(-0.615982357748271), RC(0.561474126687165), RC(0.266911112787025), RC(0.219775952207137), RC(0.387847665451514), RC(0.612483137773236), RC(0.330027015806089), RC(-0.6965298655714), RC(0.200252661187742)};
677 const PetscReal bembed[11] =
678 {RC(0.0), RC(0.180185524442613), RC(-0.628869710835338), RC(0.186185675988647), RC(0.0484716652630425), RC(0.203927720607141), RC(0.44041662512573), RC(0.615710527731245), RC(0.0689648839032607), RC(-0.253599870605903), RC(0.138606958379488)};
679 PetscCall(TSDIRKRegister(TSDIRKS7511SAL, 7, 11, &A[0][0], b, NULL, bembed, 1, b));
680 }
681 {
682 // DIRK(13,8)[1]A[(14,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
683 const PetscReal A[14][14] = {
684 {RC(0.421050745442291), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
685 {RC(-0.0761079419591268), RC(0.264353986580857), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
686 {RC(0.0727106904170694), RC(-0.204265976977285), RC(0.181608196544136), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
687 {RC(0.55763054816611), RC(-0.409773579543499), RC(0.510926516886944), RC(0.259892204518476), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
688 {RC(0.0228083864844437), RC(-0.445569051836454), RC(-0.0915242778636248), RC(0.00450055909321655), RC(0.6397807199983), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
689 {RC(-0.135945849505152), RC(0.0946509646963754), RC(-0.236110197279175), RC(0.00318944206456517), RC(0.255453021028118), RC(0.174805219173446), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
690 {RC(-0.147960260670772), RC(-0.402188192230535), RC(-0.703014530043888), RC(0.00941974677418186), RC(0.885747111289207), RC(0.261314066449028), RC(0.16307697503668), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
691 {RC(0.165597241042244), RC(0.824182962188923), RC(-0.0280136160783609), RC(0.282372386631758), RC(-0.957721354131182), RC(0.489439550159977), RC(0.170094415598103), RC(0.0522519785718563), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
692 {RC(0.0335292011495618), RC(0.575750388029166), RC(0.223289855356637), RC(-0.00317458833242804), RC(-0.112890382135193), RC(-0.419809267954284), RC(0.0466136902102104), RC(-0.00115413813041085), RC(0.109685363692383), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
693 {RC(-0.0512616878252355), RC(0.699261265830807), RC(-0.117939611738769), RC(0.0021745241931243), RC(-0.00932826702640947), RC(-0.267575057469428), RC(0.126949139814065), RC(0.00330353204502163), RC(0.185949445053766), RC(0.0938215615963721), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
694 {RC(-0.106521517960343), RC(0.41835889096168), RC(0.353585905881916), RC(-0.0746474161579599), RC(-0.015450626460289), RC(-0.46224659192275), RC(-0.0576406327329181), RC(-0.00712066942504018), RC(0.377776558014452), RC(0.36890054338294), RC(0.0618488746331837), RC(0.0), RC(0.0), RC(0.0) },
695 {RC(-0.163079104890997), RC(0.644561721693806), RC(0.636968661639572), RC(-0.122346720085377), RC(-0.333062564990312), RC(-0.3054226490478), RC(-0.357820712828352), RC(-0.0125510510334706), RC(0.371263681186311), RC(0.371979640363694), RC(0.0531090658708968), RC(0.0518279459132049), RC(0.0), RC(0.0) },
696 {RC(0.579993784455521), RC(-0.188833728676494), RC(0.999975696843775), RC(0.0572810855901161), RC(-0.264374735003671), RC(0.165091739976854), RC(-0.546675809010452), RC(-0.0283821822291982), RC(-0.102639860418374), RC(-0.0343251040446405), RC(0.4762598462591), RC(-0.304153104931261), RC(0.0953911855943621), RC(0.0) },
697 {RC(0.0848552694007844), RC(0.287193912340074), RC(0.543683503004232), RC(-0.081311059300692), RC(-0.0328661289388557), RC(-0.323456834372922), RC(-0.240378871658975), RC(-0.0189913019930369), RC(0.220663114082036), RC(0.253029984360864), RC(0.252011799370563), RC(-0.154882222605423), RC(0.0315202264687415), RC(0.0514095812104714)}
698 };
699 const PetscReal b[14] = {RC(0.0), RC(0.516650324205117), RC(0.0773227217357826), RC(-0.12474204666975), RC(-0.0241052115180679), RC(-0.325821145180359), RC(0.0907237460123951), RC(0.0459271880596652), RC(0.221012259404702), RC(0.235510906761942), RC(0.491109674204385), RC(-0.323506525837343), RC(0.119918108821531), RC(0.0)};
700 const PetscReal bembed[14] = {RC(2.32345691433618e-16), RC(0.499150900944401), RC(0.080991997189243), RC(-0.0359440417166322), RC(-0.0258910397441454), RC(-0.304540350278636), RC(0.0836627473632563),
701 RC(0.0417664613347638), RC(0.223636394275293), RC(0.231569156867596), RC(0.240526201277663), RC(-0.222933582911926), RC(-0.0111479879597561), RC(0.19915314335888)};
702 PetscCall(TSDIRKRegister(TSDIRK8614A, 8, 14, &A[0][0], b, NULL, bembed, 1, b));
703 }
704 {
705 // DIRK(15,8)[1]SAL[(16,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
706 const PetscReal A[16][16] = {
707 {RC(0.498904981271193), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
708 {RC(-0.303806037341816), RC(0.886299445992379), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
709 {RC(-0.581440223471476), RC(0.371003719460259), RC(0.43844717752802), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
710 {RC(0.531852638870051), RC(-0.339363014907108), RC(0.422373239795441), RC(0.223854203543397), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
711 {RC(0.118517891868867), RC(-0.0756235584174296), RC(-0.0864284870668712), RC(0.000536692838658312), RC(0.10101418329932), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
712 {RC(0.218733626116401), RC(-0.139568928299635), RC(0.30473612813488), RC(0.00354038623073564), RC(0.0932085751160559), RC(0.140161806097591), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
713 {RC(0.0692944686081835), RC(-0.0442152168939502), RC(-0.0903375348855603), RC(0.00259030241156141), RC(0.204514233679515), RC(-0.0245383758960002), RC(0.199289437094059), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
714 {RC(0.990640016505571), RC(-0.632104756315967), RC(0.856971425234221), RC(0.174494099232246), RC(-0.113715829680145), RC(-0.151494045307366), RC(-0.438268629569005), RC(0.120578398912139), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
715 {RC(-0.099415677713136), RC(0.211832014309207), RC(-0.245998265866888), RC(-0.182249672235861), RC(0.167897635713799), RC(0.212850335030069), RC(-0.391739299440123), RC(-0.0118718506876767), RC(0.526293701659093), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
716 {RC(0.383983914845461), RC(-0.245011361219604), RC(0.46717278554955), RC(-0.0361272447593202), RC(0.0742234660511333), RC(-0.0474816271948766), RC(-0.229859978525756), RC(0.0516283729206322), RC(0.0), RC(0.193823890777594), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
717 {RC(0.0967855003180134), RC(-0.0481037037916184), RC(0.191268138832434), RC(0.234977164564126), RC(0.0620265921753097), RC(0.403432826534738), RC(0.152403846687238), RC(-0.118420429237746), RC(0.0582141598685892), RC(-0.13924540906863), RC(0.106661313117545), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
718 {RC(0.133941307432055), RC(-0.0722076602896254), RC(0.217086297689275), RC(0.00495499602192887), RC(0.0306090174933995), RC(0.26483526755746), RC(0.204442440745605), RC(0.196883395136708), RC(0.056527012583996), RC(-0.150216381356784), RC(-0.217209415757333), RC(0.330353722743315), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
719 {RC(0.157014274561299), RC(-0.0883810256381874), RC(0.117193033885034), RC(-0.0362304243769466), RC(0.0169030211466111), RC(-0.169835753576141), RC(0.399749979234113), RC(0.31806704093008), RC(0.050340008347693), RC(0.120284837472214), RC(-0.235313193645423), RC(0.232488522208926), RC(0.117719679450729), RC(0.0), RC(0.0), RC(0.0) },
720 {RC(0.00276453816875833), RC(-0.00366028255231782), RC(-0.331078914515559), RC(0.623377549031949), RC(0.167618142989491), RC(0.0748467945312516), RC(0.797629286699677), RC(-0.390714256799583), RC(-0.00808553925131555), RC(0.014840324980952), RC(-0.0856180410248133), RC(0.602943304937827), RC(-0.5771359338496), RC(0.112273026653282), RC(0.0), RC(0.0) },
721 {RC(0.0), RC(0.0), RC(0.085283971980307), RC(0.51334393454179), RC(0.144355978013514), RC(0.255379109487853), RC(0.225075750790524), RC(-0.343241323394982), RC(0.0), RC(0.0798250392218852), RC(0.0528824734082655), RC(-0.0830350888900362), RC(0.022567388707279), RC(-0.0592631119040204), RC(0.106825878037621), RC(0.0) },
722 {RC(0.173784481207652), RC(-0.110887906116241), RC(0.190052513365204), RC(-0.0688345422674029), RC(0.10326505079603), RC(0.267127097115219), RC(0.141703423176897), RC(0.0117966866651728), RC(-6.65650091812762e-15), RC(-0.0213725083662519), RC(-0.00931148598712566), RC(-0.10007679077114), RC(0.123471797451553), RC(0.00203684241073055), RC(-0.0294320891781173), RC(0.195746619921528)}
723 };
724 const PetscReal b[16] = {RC(0.0), RC(0.0), RC(0.085283971980307), RC(0.51334393454179), RC(0.144355978013514), RC(0.255379109487853), RC(0.225075750790524), RC(-0.343241323394982), RC(0.0), RC(0.0798250392218852), RC(0.0528824734082655), RC(-0.0830350888900362), RC(0.022567388707279), RC(-0.0592631119040204), RC(0.106825878037621), RC(0.0)};
725 const PetscReal bembed[16] = {RC(-1.31988512519898e-15), RC(7.53606601764004e-16), RC(0.0886789133915965), RC(0.0968726531622137), RC(0.143815375874267), RC(0.335214773313601), RC(0.221862366978063), RC(-0.147408947987273),
726 RC(4.16297599203445e-16), RC(0.000727276166520566), RC(-0.00284892677941246), RC(0.00512492274297611), RC(-0.000275595071215218), RC(0.0136014719350733), RC(0.0165190013607726), RC(0.228116714912817)};
727 PetscCall(TSDIRKRegister(TSDIRK8616SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b));
728 }
729 {
730 // ESDIRK(16,8)[2]SAL[(16,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs
731 const PetscReal A[16][16] = {
732 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
733 {RC(0.117318819358521), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
734 {RC(0.0557014605974616), RC(0.385525646638742), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
735 {RC(0.063493276428895), RC(0.373556126263681), RC(0.0082994166438953), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
736 {RC(0.0961351856230088), RC(0.335558324517178), RC(0.207077765910132), RC(-0.0581917140797146), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
737 {RC(0.0497669214238319), RC(0.384288616546039), RC(0.0821728117583936), RC(0.120337007107103), RC(0.202262782645888), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
738 {RC(0.00626710666809847), RC(0.496491452640725), RC(-0.111303249827358), RC(0.170478821683603), RC(0.166517073971103), RC(-0.0328669811542241), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
739 {RC(0.0463439767281591), RC(0.00306724391019652), RC(-0.00816305222386205), RC(-0.0353302599538294), RC(0.0139313601702569), RC(-0.00992014507967429), RC(0.0210087909090165), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
740 {RC(0.111574049232048), RC(0.467639166482209), RC(0.237773114804619), RC(0.0798895699267508), RC(0.109580615914593), RC(0.0307353103825936), RC(-0.0404391509541147), RC(-0.16942110744293), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
741 {RC(-0.0107072484863877), RC(-0.231376703354252), RC(0.017541113036611), RC(0.144871527682418), RC(-0.041855459769806), RC(0.0841832168332261), RC(-0.0850020937282192), RC(0.486170343825899), RC(-0.0526717116822739), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
742 {RC(-0.0142238262314935), RC(0.14752923682514), RC(0.238235830732566), RC(0.037950291904103), RC(0.252075123381518), RC(0.0474266904224567), RC(-0.00363139069342027), RC(0.274081442388563), RC(-0.0599166970745255), RC(-0.0527138812389185), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
743 {RC(-0.11837020183211), RC(-0.635712481821264), RC(0.239738832602538), RC(0.330058936651707), RC(-0.325784087988237), RC(-0.0506514314589253), RC(-0.281914404487009), RC(0.852596345144291), RC(0.651444614298805), RC(-0.103476387303591), RC(-0.354835880209975), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0) },
744 {RC(-0.00458164025442349), RC(0.296219694015248), RC(0.322146049419995), RC(0.15917778285238), RC(0.284864871688843), RC(0.185509526463076), RC(-0.0784621067883274), RC(0.166312223692047), RC(-0.284152486083397), RC(-0.357125104338944), RC(0.078437074055306), RC(0.0884129667114481), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0) },
745 {RC(-0.0545561913848106), RC(0.675785423442753), RC(0.423066443201941), RC(-0.000165300126841193), RC(0.104252994793763), RC(-0.105763019303021), RC(-0.15988308809318), RC(0.0515050001032011), RC(0.56013979290924), RC(-0.45781539708603), RC(-0.255870699752664), RC(0.026960254296416), RC(-0.0721245985053681), RC(0.117318819358521), RC(0.0), RC(0.0) },
746 {RC(0.0649253995775223), RC(-0.0216056457922249), RC(-0.073738139377975), RC(0.0931033310077225), RC(-0.0194339577299149), RC(-0.0879623837313009), RC(0.057125517179467), RC(0.205120850488097), RC(0.132576503537441), RC(0.489416890627328), RC(-0.1106765720501), RC(-0.081038793996096), RC(0.0606031613503788), RC(-0.00241467937442272), RC(0.117318819358521), RC(0.0) },
747 {RC(0.0459979286336779), RC(0.0780075394482806), RC(0.015021874148058), RC(0.195180277284195), RC(-0.00246643310153235), RC(0.0473977117068314), RC(-0.0682773558610363), RC(0.19568019123878), RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)}
748 };
749 const PetscReal b[16] = {RC(0.0459979286336779), RC(0.0780075394482806), RC(0.015021874148058), RC(0.195180277284195), RC(-0.00246643310153235), RC(0.0473977117068314), RC(-0.0682773558610363), RC(0.19568019123878),
750 RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)};
751 const PetscReal bembed[16] = {RC(0.0603373529853206), RC(0.175453809423998), RC(0.0537707777611352), RC(0.195309248607308), RC(0.0135893741970232), RC(-0.0221160259296707), RC(-0.00726526156430691), RC(0.102961059369124),
752 RC(0.000900215457460583), RC(0.0547959465692338), RC(-0.334995726863153), RC(0.0464409662093384), RC(0.301388101652194), RC(0.00524851570622031), RC(0.229538601845236), RC(0.124643044573514)};
753 PetscCall(TSDIRKRegister(TSDIRKES8516SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b));
754 }
755
756 /* Additive methods */
757 {
758 const PetscReal A[3][3] = {
759 {0.0, 0.0, 0.0},
760 {0.0, 0.0, 0.0},
761 {0.0, 0.5, 0.0}
762 };
763 const PetscReal At[3][3] = {
764 {1.0, 0.0, 0.0},
765 {0.0, 0.5, 0.0},
766 {0.0, 0.5, 0.5}
767 };
768 const PetscReal b[3] = {0.0, 0.5, 0.5};
769 const PetscReal bembedt[3] = {1.0, 0.0, 0.0};
770 PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL));
771 }
772 {
773 const PetscReal A[2][2] = {
774 {0.0, 0.0},
775 {0.5, 0.0}
776 };
777 const PetscReal At[2][2] = {
778 {0.0, 0.0},
779 {0.0, 0.5}
780 };
781 const PetscReal b[2] = {0.0, 1.0};
782 const PetscReal bembedt[2] = {0.5, 0.5};
783 /* binterpt[2][2] = {{1.0,-1.0},{0.0,1.0}}; second order dense output has poor stability properties and hence it is not currently in use */
784 PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL));
785 }
786 {
787 const PetscReal A[2][2] = {
788 {0.0, 0.0},
789 {1.0, 0.0}
790 };
791 const PetscReal At[2][2] = {
792 {0.0, 0.0},
793 {0.5, 0.5}
794 };
795 const PetscReal b[2] = {0.5, 0.5};
796 const PetscReal bembedt[2] = {0.0, 1.0};
797 /* binterpt[2][2] = {{1.0,-0.5},{0.0,0.5}} second order dense output has poor stability properties and hence it is not currently in use */
798 PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL));
799 }
800 {
801 const PetscReal A[2][2] = {
802 {0.0, 0.0},
803 {1.0, 0.0}
804 };
805 const PetscReal At[2][2] = {
806 {us2, 0.0},
807 {1.0 - 2.0 * us2, us2}
808 };
809 const PetscReal b[2] = {0.5, 0.5};
810 const PetscReal bembedt[2] = {0.0, 1.0};
811 const PetscReal binterpt[2][2] = {
812 {(us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))},
813 {1 - (us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))}
814 };
815 const PetscReal binterp[2][2] = {
816 {1.0, -0.5},
817 {0.0, 0.5 }
818 };
819 PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0]));
820 }
821 {
822 const PetscReal A[3][3] = {
823 {0, 0, 0},
824 {2 - s2, 0, 0},
825 {0.5, 0.5, 0}
826 };
827 const PetscReal At[3][3] = {
828 {0, 0, 0 },
829 {1 - 1 / s2, 1 - 1 / s2, 0 },
830 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2}
831 };
832 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)};
833 const PetscReal binterpt[3][2] = {
834 {1.0 / s2, -1.0 / (2.0 * s2)},
835 {1.0 / s2, -1.0 / (2.0 * s2)},
836 {1.0 - s2, 1.0 / s2 }
837 };
838 PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
839 }
840 {
841 const PetscReal A[3][3] = {
842 {0, 0, 0},
843 {2 - s2, 0, 0},
844 {0.75, 0.25, 0}
845 };
846 const PetscReal At[3][3] = {
847 {0, 0, 0 },
848 {1 - 1 / s2, 1 - 1 / s2, 0 },
849 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2}
850 };
851 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)};
852 const PetscReal binterpt[3][2] = {
853 {1.0 / s2, -1.0 / (2.0 * s2)},
854 {1.0 / s2, -1.0 / (2.0 * s2)},
855 {1.0 - s2, 1.0 / s2 }
856 };
857 PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
858 }
859 { /* Optimal for linear implicit part */
860 const PetscReal A[3][3] = {
861 {0, 0, 0},
862 {2 - s2, 0, 0},
863 {(3 - 2 * s2) / 6, (3 + 2 * s2) / 6, 0}
864 };
865 const PetscReal At[3][3] = {
866 {0, 0, 0 },
867 {1 - 1 / s2, 1 - 1 / s2, 0 },
868 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2}
869 };
870 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)};
871 const PetscReal binterpt[3][2] = {
872 {1.0 / s2, -1.0 / (2.0 * s2)},
873 {1.0 / s2, -1.0 / (2.0 * s2)},
874 {1.0 - s2, 1.0 / s2 }
875 };
876 PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
877 }
878 { /* Optimal for linear implicit part */
879 const PetscReal A[3][3] = {
880 {0, 0, 0},
881 {0.5, 0, 0},
882 {0.5, 0.5, 0}
883 };
884 const PetscReal At[3][3] = {
885 {0.25, 0, 0 },
886 {0, 0.25, 0 },
887 {1. / 3, 1. / 3, 1. / 3}
888 };
889 PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL));
890 }
891 {
892 const PetscReal A[4][4] = {
893 {0, 0, 0, 0},
894 {1767732205903. / 2027836641118., 0, 0, 0},
895 {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0},
896 {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0}
897 };
898 const PetscReal At[4][4] = {
899 {0, 0, 0, 0 },
900 {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0 },
901 {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0 },
902 {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.}
903 };
904 const PetscReal bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.};
905 const PetscReal binterpt[4][2] = {
906 {4655552711362. / 22874653954995., -215264564351. / 13552729205753. },
907 {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119. },
908 {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.},
909 {584795268549. / 6622622206610., 2508943948391. / 7218656332882. }
910 };
911 PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
912 }
913 {
914 const PetscReal A[5][5] = {
915 {0, 0, 0, 0, 0},
916 {1. / 2, 0, 0, 0, 0},
917 {11. / 18, 1. / 18, 0, 0, 0},
918 {5. / 6, -5. / 6, .5, 0, 0},
919 {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0}
920 };
921 const PetscReal At[5][5] = {
922 {0, 0, 0, 0, 0 },
923 {0, 1. / 2, 0, 0, 0 },
924 {0, 1. / 6, 1. / 2, 0, 0 },
925 {0, -1. / 2, 1. / 2, 1. / 2, 0 },
926 {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2}
927 };
928 PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL));
929 }
930 {
931 const PetscReal A[5][5] = {
932 {0, 0, 0, 0, 0},
933 {1, 0, 0, 0, 0},
934 {4. / 9, 2. / 9, 0, 0, 0},
935 {1. / 4, 0, 3. / 4, 0, 0},
936 {1. / 4, 0, 3. / 5, 0, 0}
937 };
938 const PetscReal At[5][5] = {
939 {0, 0, 0, 0, 0 },
940 {.5, .5, 0, 0, 0 },
941 {5. / 18, -1. / 9, .5, 0, 0 },
942 {.5, 0, 0, .5, 0 },
943 {.25, 0, .75, -.5, .5}
944 };
945 PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL));
946 }
947 {
948 const PetscReal A[6][6] = {
949 {0, 0, 0, 0, 0, 0},
950 {1. / 2, 0, 0, 0, 0, 0},
951 {13861. / 62500., 6889. / 62500., 0, 0, 0, 0},
952 {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0},
953 {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0},
954 {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0}
955 };
956 const PetscReal At[6][6] = {
957 {0, 0, 0, 0, 0, 0 },
958 {1. / 4, 1. / 4, 0, 0, 0, 0 },
959 {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0 },
960 {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0 },
961 {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0 },
962 {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4}
963 };
964 const PetscReal bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.};
965 const PetscReal binterpt[6][3] = {
966 {6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025. },
967 {0, 0, 0 },
968 {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035. },
969 {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.},
970 {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469. },
971 {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.}
972 };
973 PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL));
974 }
975 {
976 const PetscReal A[8][8] = {
977 {0, 0, 0, 0, 0, 0, 0, 0},
978 {41. / 100, 0, 0, 0, 0, 0, 0, 0},
979 {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0},
980 {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0},
981 {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0},
982 {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0},
983 {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0},
984 {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0}
985 };
986 const PetscReal At[8][8] = {
987 {0, 0, 0, 0, 0, 0, 0, 0 },
988 {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0 },
989 {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0 },
990 {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0 },
991 {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0 },
992 {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0 },
993 {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0 },
994 {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.}
995 };
996 const PetscReal bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.};
997 const PetscReal binterpt[8][3] = {
998 {-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439. },
999 {0, 0, 0 },
1000 {0, 0, 0 },
1001 {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.},
1002 {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.},
1003 {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.},
1004 {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.},
1005 {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460. }
1006 };
1007 PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL));
1008 }
1009 #undef RC
1010 #undef us2
1011 #undef s2
1012 PetscFunctionReturn(PETSC_SUCCESS);
1013 }
1014
1015 /*@C
1016 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`.
1017
1018 Not Collective
1019
1020 Level: advanced
1021
1022 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()`
1023 @*/
TSARKIMEXRegisterDestroy(void)1024 PetscErrorCode TSARKIMEXRegisterDestroy(void)
1025 {
1026 ARKTableauLink link;
1027
1028 PetscFunctionBegin;
1029 while ((link = ARKTableauList)) {
1030 ARKTableau t = &link->tab;
1031 ARKTableauList = link->next;
1032 PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c));
1033 PetscCall(PetscFree2(t->bembedt, t->bembed));
1034 PetscCall(PetscFree2(t->binterpt, t->binterp));
1035 PetscCall(PetscFree(t->name));
1036 PetscCall(PetscFree(link));
1037 }
1038 TSARKIMEXRegisterAllCalled = PETSC_FALSE;
1039 PetscFunctionReturn(PETSC_SUCCESS);
1040 }
1041
1042 /*@C
1043 TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called
1044 from `TSInitializePackage()`.
1045
1046 Level: developer
1047
1048 .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()`
1049 @*/
TSARKIMEXInitializePackage(void)1050 PetscErrorCode TSARKIMEXInitializePackage(void)
1051 {
1052 PetscFunctionBegin;
1053 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
1054 TSARKIMEXPackageInitialized = PETSC_TRUE;
1055 PetscCall(TSARKIMEXRegisterAll());
1056 PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage));
1057 PetscFunctionReturn(PETSC_SUCCESS);
1058 }
1059
1060 /*@C
1061 TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is
1062 called from `PetscFinalize()`.
1063
1064 Level: developer
1065
1066 .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()`
1067 @*/
TSARKIMEXFinalizePackage(void)1068 PetscErrorCode TSARKIMEXFinalizePackage(void)
1069 {
1070 PetscFunctionBegin;
1071 TSARKIMEXPackageInitialized = PETSC_FALSE;
1072 PetscCall(TSARKIMEXRegisterDestroy());
1073 PetscFunctionReturn(PETSC_SUCCESS);
1074 }
1075
1076 /*@C
1077 TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
1078
1079 Logically Collective
1080
1081 Input Parameters:
1082 + name - identifier for method
1083 . order - approximation order of method
1084 . s - number of stages, this is the dimension of the matrices below
1085 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major)
1086 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At)
1087 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At)
1088 . A - Non-stiff stage coefficients (dimension s*s, row-major)
1089 . b - Non-stiff step completion table (dimension s; NULL to use last row of At)
1090 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A)
1091 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available)
1092 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided)
1093 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp
1094 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp)
1095 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt)
1096
1097 Level: advanced
1098
1099 Note:
1100 Several `TSARKIMEX` methods are provided, this function is only needed to create new methods.
1101
1102 .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS`
1103 @*/
TSARKIMEXRegister(TSARKIMEXType name,PetscInt order,PetscInt s,const PetscReal At[],const PetscReal bt[],const PetscReal ct[],const PetscReal A[],const PetscReal b[],const PetscReal c[],const PetscReal bembedt[],const PetscReal bembed[],PetscInt pinterp,const PetscReal binterpt[],const PetscReal binterp[])1104 PetscErrorCode TSARKIMEXRegister(TSARKIMEXType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal A[], const PetscReal b[], const PetscReal c[], const PetscReal bembedt[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[], const PetscReal binterp[])
1105 {
1106 ARKTableauLink link;
1107 ARKTableau t;
1108 PetscInt i, j;
1109
1110 PetscFunctionBegin;
1111 PetscCheck(s > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected number of stages s %" PetscInt_FMT " > 0", s);
1112 PetscCall(TSARKIMEXInitializePackage());
1113 for (link = ARKTableauList; link; link = link->next) {
1114 PetscBool match;
1115
1116 PetscCall(PetscStrcmp(link->tab.name, name, &match));
1117 PetscCheck(!match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Method with name \"%s\" already registered", name);
1118 }
1119 PetscCall(PetscNew(&link));
1120 t = &link->tab;
1121 PetscCall(PetscStrallocpy(name, &t->name));
1122 t->order = order;
1123 t->s = s;
1124 PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c));
1125 PetscCall(PetscArraycpy(t->At, At, s * s));
1126 if (A) {
1127 PetscCall(PetscArraycpy(t->A, A, s * s));
1128 t->additive = PETSC_TRUE;
1129 }
1130
1131 if (bt) PetscCall(PetscArraycpy(t->bt, bt, s));
1132 else
1133 for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i];
1134
1135 if (t->additive) {
1136 if (b) PetscCall(PetscArraycpy(t->b, b, s));
1137 else
1138 for (i = 0; i < s; i++) t->b[i] = t->bt[i];
1139 } else PetscCall(PetscArrayzero(t->b, s));
1140
1141 if (ct) PetscCall(PetscArraycpy(t->ct, ct, s));
1142 else
1143 for (i = 0; i < s; i++)
1144 for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j];
1145
1146 if (t->additive) {
1147 if (c) PetscCall(PetscArraycpy(t->c, c, s));
1148 else
1149 for (i = 0; i < s; i++)
1150 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j];
1151 } else PetscCall(PetscArrayzero(t->c, s));
1152
1153 t->stiffly_accurate = PETSC_TRUE;
1154 for (i = 0; i < s; i++)
1155 if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE;
1156
1157 t->explicit_first_stage = PETSC_TRUE;
1158 for (i = 0; i < s; i++)
1159 if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE;
1160
1161 /* def of FSAL can be made more precise */
1162 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate);
1163
1164 if (bembedt) {
1165 PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed));
1166 PetscCall(PetscArraycpy(t->bembedt, bembedt, s));
1167 PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s));
1168 }
1169
1170 t->pinterp = pinterp;
1171 PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp));
1172 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp));
1173 PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp));
1174
1175 link->next = ARKTableauList;
1176 ARKTableauList = link;
1177 PetscFunctionReturn(PETSC_SUCCESS);
1178 }
1179
1180 /*@C
1181 TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation
1182
1183 Logically Collective.
1184
1185 Input Parameters:
1186 + name - identifier for method
1187 . order - approximation order of method
1188 . s - number of stages, this is the dimension of the matrices below
1189 . At - Butcher table of stage coefficients (dimension `s`*`s`, row-major order)
1190 . bt - Butcher table for completing the step (dimension `s`; pass `NULL` to use the last row of `At`)
1191 . ct - Abscissa of each stage (dimension s, NULL to use row sums of At)
1192 . bembedt - Stiff part of completion table for embedded method (dimension s; `NULL` if not available)
1193 . pinterp - Order of the interpolation scheme, equal to the number of columns of `binterpt` and `binterp`
1194 - binterpt - Coefficients of the interpolation formula (dimension s*pinterp)
1195
1196 Level: advanced
1197
1198 Note:
1199 Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods.
1200
1201 .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS`
1202 @*/
TSDIRKRegister(TSDIRKType name,PetscInt order,PetscInt s,const PetscReal At[],const PetscReal bt[],const PetscReal ct[],const PetscReal bembedt[],PetscInt pinterp,const PetscReal binterpt[])1203 PetscErrorCode TSDIRKRegister(TSDIRKType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal bembedt[], PetscInt pinterp, const PetscReal binterpt[])
1204 {
1205 PetscFunctionBegin;
1206 PetscCall(TSARKIMEXRegister(name, order, s, At, bt, ct, NULL, NULL, NULL, bembedt, NULL, pinterp, binterpt, NULL));
1207 PetscFunctionReturn(PETSC_SUCCESS);
1208 }
1209
1210 /*
1211 The step completion formula is
1212
1213 x1 = x0 - h bt^T YdotI + h b^T YdotRHS
1214
1215 This function can be called before or after ts->vec_sol has been updated.
1216 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order.
1217 We can write
1218
1219 x1e = x0 - h bet^T YdotI + h be^T YdotRHS
1220 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS
1221 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS
1222
1223 so we can evaluate the method with different order even after the step has been optimistically completed.
1224 */
TSEvaluateStep_ARKIMEX(TS ts,PetscInt order,Vec X,PetscBool * done)1225 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done)
1226 {
1227 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1228 ARKTableau tab = ark->tableau;
1229 PetscScalar *w = ark->work;
1230 PetscReal h;
1231 PetscInt s = tab->s, j;
1232 PetscBool hasE;
1233
1234 PetscFunctionBegin;
1235 switch (ark->status) {
1236 case TS_STEP_INCOMPLETE:
1237 case TS_STEP_PENDING:
1238 h = ts->time_step;
1239 break;
1240 case TS_STEP_COMPLETE:
1241 h = ts->ptime - ts->ptime_prev;
1242 break;
1243 default:
1244 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1245 }
1246 if (order == tab->order) {
1247 if (ark->status == TS_STEP_INCOMPLETE) {
1248 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */
1249 PetscCall(VecCopy(ark->Y[s - 1], X));
1250 } else { /* Use the standard completion formula (bt,b) */
1251 PetscCall(VecCopy(ts->vec_sol, X));
1252 for (j = 0; j < s; j++) w[j] = h * tab->bt[j];
1253 PetscCall(VecMAXPY(X, s, w, ark->YdotI));
1254 if (tab->additive && ark->imex) { /* Method is IMEX, complete the explicit formula */
1255 PetscCall(TSHasRHSFunction(ts, &hasE));
1256 if (hasE) {
1257 for (j = 0; j < s; j++) w[j] = h * tab->b[j];
1258 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS));
1259 }
1260 }
1261 }
1262 } else PetscCall(VecCopy(ts->vec_sol, X));
1263 if (done) *done = PETSC_TRUE;
1264 PetscFunctionReturn(PETSC_SUCCESS);
1265 } else if (order == tab->order - 1) {
1266 if (!tab->bembedt) goto unavailable;
1267 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */
1268 PetscCall(VecCopy(ts->vec_sol, X));
1269 for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j];
1270 PetscCall(VecMAXPY(X, s, w, ark->YdotI));
1271 if (tab->additive) {
1272 PetscCall(TSHasRHSFunction(ts, &hasE));
1273 if (hasE) {
1274 for (j = 0; j < s; j++) w[j] = h * tab->bembed[j];
1275 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS));
1276 }
1277 }
1278 } else { /* Rollback and re-complete using (bet-be,be-b) */
1279 PetscCall(VecCopy(ts->vec_sol, X));
1280 for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]);
1281 PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI));
1282 if (tab->additive) {
1283 PetscCall(TSHasRHSFunction(ts, &hasE));
1284 if (hasE) {
1285 for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]);
1286 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS));
1287 }
1288 }
1289 }
1290 if (done) *done = PETSC_TRUE;
1291 PetscFunctionReturn(PETSC_SUCCESS);
1292 }
1293 unavailable:
1294 PetscCheck(done, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "ARKIMEX '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.",
1295 tab->name, tab->order, order);
1296 *done = PETSC_FALSE;
1297 PetscFunctionReturn(PETSC_SUCCESS);
1298 }
1299
TSARKIMEXTestMassIdentity(TS ts,PetscBool * id)1300 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id)
1301 {
1302 Vec Udot, Y1, Y2;
1303 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1304 PetscReal norm;
1305
1306 PetscFunctionBegin;
1307 PetscCall(VecDuplicate(ts->vec_sol, &Udot));
1308 PetscCall(VecDuplicate(ts->vec_sol, &Y1));
1309 PetscCall(VecDuplicate(ts->vec_sol, &Y2));
1310 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex));
1311 PetscCall(VecSetRandom(Udot, NULL));
1312 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex));
1313 PetscCall(VecAXPY(Y2, -1.0, Y1));
1314 PetscCall(VecAXPY(Y2, -1.0, Udot));
1315 PetscCall(VecNorm(Y2, NORM_2, &norm));
1316 if (norm < 100.0 * PETSC_MACHINE_EPSILON) {
1317 *id = PETSC_TRUE;
1318 } else {
1319 *id = PETSC_FALSE;
1320 PetscCall(PetscInfo(ts, "IFunction(Udot = random) - IFunction(Udot = 0) is not near Udot, %g, suspect mass matrix implied in IFunction() is not the identity as required\n", (double)norm));
1321 }
1322 PetscCall(VecDestroy(&Udot));
1323 PetscCall(VecDestroy(&Y1));
1324 PetscCall(VecDestroy(&Y2));
1325 PetscFunctionReturn(PETSC_SUCCESS);
1326 }
1327
1328 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS, PetscReal, Vec, IS *);
1329
TSStep_ARKIMEX(TS ts)1330 static PetscErrorCode TSStep_ARKIMEX(TS ts)
1331 {
1332 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1333 ARKTableau tab = ark->tableau;
1334 const PetscInt s = tab->s;
1335 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c;
1336 PetscScalar *w = ark->work;
1337 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z;
1338 PetscBool extrapolate = ark->extrapolate;
1339 TSAdapt adapt;
1340 SNES snes;
1341 PetscInt i, j, its, lits;
1342 PetscInt rejections = 0;
1343 PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE;
1344 PetscReal next_time_step = ts->time_step;
1345
1346 PetscFunctionBegin;
1347 if (ark->extrapolate && !ark->Y_prev) {
1348 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev));
1349 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev));
1350 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev));
1351 }
1352
1353 if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE));
1354 if (!hasE) dirk = PETSC_TRUE;
1355
1356 if (!ts->steprollback) {
1357 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */
1358 PetscCall(VecCopy(YdotI[s - 1], Ydot0));
1359 }
1360 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */
1361 for (i = 0; i < s; i++) {
1362 PetscCall(VecCopy(Y[i], ark->Y_prev[i]));
1363 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i]));
1364 if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i]));
1365 }
1366 }
1367 }
1368
1369 /*
1370 For fully implicit formulations we must solve the equations
1371
1372 F(t_n,x_n,xdot) = 0
1373
1374 for the explicit first stage.
1375 Here we call SNESSolve using PETSC_MAX_REAL as shift to flag it.
1376 Special handling is inside SNESTSFormFunction_ARKIMEX and SNESTSFormJacobian_ARKIMEX
1377 We compute Ydot0 if we restart the step or if we resized the problem after remeshing
1378 */
1379 if (dirk && tab->explicit_first_stage && (ts->steprestart || ts->stepresize)) {
1380 ark->scoeff = PETSC_MAX_REAL;
1381 PetscCall(VecCopy(ts->vec_sol, Z));
1382 if (!ark->alg_is) {
1383 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ts->ptime, Z, &ark->alg_is));
1384 PetscCall(ISViewFromOptions(ark->alg_is, (PetscObject)ts, "-ts_arkimex_algebraic_is_view"));
1385 }
1386 PetscCall(TSGetSNES(ts, &snes));
1387 PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, 1));
1388 PetscCall(SNESSolve(snes, NULL, Ydot0));
1389 if (ark->alg_is) PetscCall(VecISSet(Ydot0, ark->alg_is, 0.0));
1390 PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, -1));
1391 }
1392
1393 /* For IMEX we compute a step */
1394 if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) {
1395 TS ts_start;
1396 if (PetscDefined(USE_DEBUG) && hasE) {
1397 PetscBool id = PETSC_FALSE;
1398 PetscCall(TSARKIMEXTestMassIdentity(ts, &id));
1399 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "This scheme requires an identity mass matrix, however the TSIFunctionFn you provided does not utilize an identity mass matrix");
1400 }
1401 PetscCall(TSClone(ts, &ts_start));
1402 PetscCall(TSSetSolution(ts_start, ts->vec_sol));
1403 PetscCall(TSSetTime(ts_start, ts->ptime));
1404 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1));
1405 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step));
1406 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER));
1407 PetscCall(TSSetTimeStep(ts_start, ts->time_step));
1408 PetscCall(TSSetType(ts_start, TSARKIMEX));
1409 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE));
1410 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE));
1411
1412 PetscCall(TSRestartStep(ts_start));
1413 PetscCall(TSSolve(ts_start, ts->vec_sol));
1414 PetscCall(TSGetTime(ts_start, &ts->ptime));
1415 PetscCall(TSGetTimeStep(ts_start, &ts->time_step));
1416
1417 { /* Save the initial slope for the next step */
1418 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data;
1419 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0));
1420 }
1421 ts->steps++;
1422
1423 /* Set the correct TS in SNES */
1424 /* We'll try to bypass this by changing the method on the fly */
1425 {
1426 PetscCall(TSGetSNES(ts, &snes));
1427 PetscCall(TSSetSNES(ts, snes));
1428 }
1429 PetscCall(TSDestroy(&ts_start));
1430 }
1431
1432 ark->status = TS_STEP_INCOMPLETE;
1433 while (!ts->reason && ark->status != TS_STEP_COMPLETE) {
1434 PetscReal t = ts->ptime;
1435 PetscReal h = ts->time_step;
1436 for (i = 0; i < s; i++) {
1437 ark->stage_time = t + h * ct[i];
1438 PetscCall(TSPreStage(ts, ark->stage_time));
1439 if (At[i * s + i] == 0) { /* This stage is explicit */
1440 PetscCheck(i == 0 || ts->equation_type < TS_EQ_IMPLICIT, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Explicit stages other than the first one are not supported for implicit problems");
1441 PetscCall(VecCopy(ts->vec_sol, Y[i]));
1442 for (j = 0; j < i; j++) w[j] = h * At[i * s + j];
1443 PetscCall(VecMAXPY(Y[i], i, w, YdotI));
1444 if (tab->additive && hasE) {
1445 for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
1446 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS));
1447 }
1448 PetscCall(TSGetSNES(ts, &snes));
1449 PetscCall(SNESResetCounters(snes));
1450 } else {
1451 ark->scoeff = 1. / At[i * s + i];
1452 /* Ydot = shift*(Y-Z) */
1453 PetscCall(VecCopy(ts->vec_sol, Z));
1454 for (j = 0; j < i; j++) w[j] = h * At[i * s + j];
1455 PetscCall(VecMAXPY(Z, i, w, YdotI));
1456 if (tab->additive && hasE) {
1457 for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
1458 PetscCall(VecMAXPY(Z, i, w, YdotRHS));
1459 }
1460 if (extrapolate && !ts->steprestart) {
1461 /* Initial guess extrapolated from previous time step stage values */
1462 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i]));
1463 } else {
1464 /* Initial guess taken from last stage */
1465 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i]));
1466 }
1467 PetscCall(TSGetSNES(ts, &snes));
1468 PetscCall(SNESSolve(snes, NULL, Y[i]));
1469 PetscCall(SNESGetIterationNumber(snes, &its));
1470 PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1471 ts->snes_its += its;
1472 ts->ksp_its += lits;
1473 PetscCall(TSGetAdapt(ts, &adapt));
1474 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok));
1475 if (!stageok) {
1476 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to
1477 * use extrapolation to initialize the solves on the next attempt. */
1478 extrapolate = PETSC_FALSE;
1479 goto reject_step;
1480 }
1481 }
1482 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) {
1483 if (i == 0 && tab->explicit_first_stage) {
1484 PetscCheck(tab->stiffly_accurate, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s is not stiffly accurate and therefore explicit-first stage methods cannot be used if the equation is implicit because the slope cannot be evaluated",
1485 ((PetscObject)ts)->type_name, ark->tableau->name);
1486 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */
1487 } else {
1488 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */
1489 }
1490 } else {
1491 if (i == 0 && tab->explicit_first_stage) {
1492 PetscCall(VecZeroEntries(Ydot));
1493 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */
1494 PetscCall(VecScale(YdotI[i], -1.0));
1495 } else {
1496 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */
1497 }
1498 if (hasE) {
1499 if (ark->imex) {
1500 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i]));
1501 } else {
1502 PetscCall(VecZeroEntries(YdotRHS[i]));
1503 }
1504 }
1505 }
1506 PetscCall(TSPostStage(ts, ark->stage_time, i, Y));
1507 }
1508
1509 ark->status = TS_STEP_INCOMPLETE;
1510 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL));
1511 ark->status = TS_STEP_PENDING;
1512 PetscCall(TSGetAdapt(ts, &adapt));
1513 PetscCall(TSAdaptCandidatesClear(adapt));
1514 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
1515 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
1516 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1517 if (!accept) { /* Roll back the current step */
1518 PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
1519 ts->time_step = next_time_step;
1520 goto reject_step;
1521 }
1522
1523 ts->ptime += ts->time_step;
1524 ts->time_step = next_time_step;
1525 break;
1526
1527 reject_step:
1528 ts->reject++;
1529 accept = PETSC_FALSE;
1530 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1531 ts->reason = TS_DIVERGED_STEP_REJECTED;
1532 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
1533 }
1534 }
1535 PetscFunctionReturn(PETSC_SUCCESS);
1536 }
1537
1538 /*
1539 This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form
1540 Udot = H(t,U) + G(t,U)
1541 This corresponds to F(t,U,Udot) = Udot-H(t,U).
1542
1543 The complete adjoint equations are
1544 (shift*I - dHdu) lambda_s[i] = 1/at[i][i] (
1545 dGdU (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j])
1546 + dHdU (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0
1547 lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j]
1548 mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i]
1549 + dGdP (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j])
1550 + dHdP (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0
1551 mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j]
1552 where shift = 1/(at[i][i]*h)
1553
1554 If at[i][i] is 0, the first equation falls back to
1555 lambda_s[i] = h (
1556 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j]
1557 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j])
1558
1559 */
TSAdjointStep_ARKIMEX(TS ts)1560 static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts)
1561 {
1562 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1563 ARKTableau tab = ark->tableau;
1564 const PetscInt s = tab->s;
1565 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt;
1566 PetscScalar *w = ark->work;
1567 Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp;
1568 Mat Jex, Jim, Jimpre;
1569 PetscInt i, j, nadj;
1570 PetscReal t = ts->ptime, stage_time_ex;
1571 PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */
1572 KSP ksp;
1573
1574 PetscFunctionBegin;
1575 ark->status = TS_STEP_INCOMPLETE;
1576 PetscCall(SNESGetKSP(ts->snes, &ksp));
1577 PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL));
1578 PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL));
1579
1580 for (i = s - 1; i >= 0; i--) {
1581 ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]);
1582 stage_time_ex = t - adjoint_time_step * (1.0 - c[i]);
1583 if (At[i * s + i] == 0) { // This stage is explicit
1584 ark->scoeff = 0.;
1585 } else {
1586 ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative
1587 }
1588 PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre));
1589 PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex));
1590 if (ts->vecs_sensip) {
1591 PetscCall(TSComputeIJacobianP(ts, ark->stage_time, Y[i], Ydot, ark->scoeff / adjoint_time_step, ts->Jacp, PETSC_TRUE)); // get dFdP (-dHdP), Ydot not really used since mass matrix is identity
1592 PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP
1593 }
1594 /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */
1595 for (nadj = 0; nadj < ts->numcost; nadj++) {
1596 /* build implicit part */
1597 PetscCall(VecSet(VecsSensiTemp[nadj], 0));
1598 if (s - i - 1 > 0) {
1599 /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */
1600 for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i];
1601 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
1602 }
1603 /* Temp = Temp - bt[i] lambda_{n+1} */
1604 PetscCall(VecAXPY(VecsSensiTemp[nadj], -bt[i], ts->vecs_sensi[nadj]));
1605 if (bt[i] || s - i - 1 > 0) {
1606 /* (shift I - dHdU) Temp */
1607 PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i]));
1608 /* cancel out shift Temp where shift=-scoeff/h */
1609 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj]));
1610 if (ts->vecs_sensip) {
1611 /* - dHdP Temp */
1612 PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj]));
1613 /* mu_n += -h dHdP Temp */
1614 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj]));
1615 }
1616 } else {
1617 PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized
1618 }
1619 /* build explicit part */
1620 PetscCall(VecSet(VecsSensiTemp[nadj], 0));
1621 if (s - i - 1 > 0) {
1622 /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */
1623 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i];
1624 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
1625 }
1626 /* Temp = Temp + b[i] lambda_{n+1} */
1627 PetscCall(VecAXPY(VecsSensiTemp[nadj], b[i], ts->vecs_sensi[nadj]));
1628 if (b[i] || s - i - 1 > 0) {
1629 /* dGdU Temp */
1630 PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i]));
1631 if (ts->vecs_sensip) {
1632 /* dGdP Temp */
1633 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj]));
1634 /* mu_n += h dGdP Temp */
1635 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj]));
1636 }
1637 }
1638 /* Build LHS for first-order adjoint */
1639 if (At[i * s + i] == 0) { // This stage is explicit
1640 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step));
1641 } else {
1642 KSP ksp;
1643 KSPConvergedReason kspreason;
1644 PetscCall(SNESGetKSP(ts->snes, &ksp));
1645 PetscCall(KSPSetOperators(ksp, Jim, Jimpre));
1646 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i]));
1647 PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i]));
1648 PetscCall(KSPGetConvergedReason(ksp, &kspreason));
1649 if (kspreason < 0) {
1650 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
1651 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj));
1652 }
1653 if (ts->vecs_sensip) {
1654 /* -dHdP lambda_s[i] */
1655 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj]));
1656 /* mu_n += h at[i][i] dHdP lambda_s[i] */
1657 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj]));
1658 }
1659 }
1660 }
1661 }
1662 for (j = 0; j < s; j++) w[j] = 1.0;
1663 for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's
1664 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s]));
1665 ark->status = TS_STEP_COMPLETE;
1666 PetscFunctionReturn(PETSC_SUCCESS);
1667 }
1668
TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X)1669 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X)
1670 {
1671 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1672 ARKTableau tab = ark->tableau;
1673 PetscInt s = tab->s, pinterp = tab->pinterp, i, j;
1674 PetscReal h;
1675 PetscReal tt, t;
1676 PetscScalar *bt = ark->work, *b = ark->work + s;
1677 const PetscReal *Bt = tab->binterpt, *B = tab->binterp;
1678
1679 PetscFunctionBegin;
1680 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name);
1681 switch (ark->status) {
1682 case TS_STEP_INCOMPLETE:
1683 case TS_STEP_PENDING:
1684 h = ts->time_step;
1685 t = (itime - ts->ptime) / h;
1686 break;
1687 case TS_STEP_COMPLETE:
1688 h = ts->ptime - ts->ptime_prev;
1689 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1690 break;
1691 default:
1692 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1693 }
1694 for (i = 0; i < s; i++) bt[i] = b[i] = 0;
1695 for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1696 for (i = 0; i < s; i++) {
1697 bt[i] += h * Bt[i * pinterp + j] * tt;
1698 b[i] += h * B[i * pinterp + j] * tt;
1699 }
1700 }
1701 PetscCall(VecCopy(ark->Y[0], X));
1702 PetscCall(VecMAXPY(X, s, bt, ark->YdotI));
1703 if (tab->additive) {
1704 PetscBool hasE;
1705 PetscCall(TSHasRHSFunction(ts, &hasE));
1706 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS));
1707 }
1708 PetscFunctionReturn(PETSC_SUCCESS);
1709 }
1710
TSExtrapolate_ARKIMEX(TS ts,PetscReal c,Vec X)1711 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X)
1712 {
1713 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1714 ARKTableau tab = ark->tableau;
1715 PetscInt s = tab->s, pinterp = tab->pinterp, i, j;
1716 PetscReal h, h_prev, t, tt;
1717 PetscScalar *bt = ark->work, *b = ark->work + s;
1718 const PetscReal *Bt = tab->binterpt, *B = tab->binterp;
1719
1720 PetscFunctionBegin;
1721 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name);
1722 h = ts->time_step;
1723 h_prev = ts->ptime - ts->ptime_prev;
1724 t = 1 + h / h_prev * c;
1725 for (i = 0; i < s; i++) bt[i] = b[i] = 0;
1726 for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1727 for (i = 0; i < s; i++) {
1728 bt[i] += h * Bt[i * pinterp + j] * tt;
1729 b[i] += h * B[i * pinterp + j] * tt;
1730 }
1731 }
1732 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored");
1733 PetscCall(VecCopy(ark->Y_prev[0], X));
1734 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev));
1735 if (tab->additive) {
1736 PetscBool hasE;
1737 PetscCall(TSHasRHSFunction(ts, &hasE));
1738 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev));
1739 }
1740 PetscFunctionReturn(PETSC_SUCCESS);
1741 }
1742
TSARKIMEXTableauReset(TS ts)1743 static PetscErrorCode TSARKIMEXTableauReset(TS ts)
1744 {
1745 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1746 ARKTableau tab = ark->tableau;
1747
1748 PetscFunctionBegin;
1749 if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1750 PetscCall(PetscFree(ark->work));
1751 PetscCall(VecDestroyVecs(tab->s, &ark->Y));
1752 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI));
1753 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS));
1754 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev));
1755 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev));
1756 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev));
1757 PetscFunctionReturn(PETSC_SUCCESS);
1758 }
1759
TSReset_ARKIMEX(TS ts)1760 static PetscErrorCode TSReset_ARKIMEX(TS ts)
1761 {
1762 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1763
1764 PetscFunctionBegin;
1765 if (ark->fastslowsplit) {
1766 PetscTryMethod(ts, "TSReset_ARKIMEX_FastSlowSplit_C", (TS), (ts));
1767 } else {
1768 PetscCall(TSARKIMEXTableauReset(ts));
1769 PetscCall(VecDestroy(&ark->Ydot));
1770 PetscCall(VecDestroy(&ark->Ydot0));
1771 PetscCall(VecDestroy(&ark->Z));
1772 PetscCall(ISDestroy(&ark->alg_is));
1773 }
1774 PetscFunctionReturn(PETSC_SUCCESS);
1775 }
1776
TSAdjointReset_ARKIMEX(TS ts)1777 static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts)
1778 {
1779 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1780 ARKTableau tab = ark->tableau;
1781
1782 PetscFunctionBegin;
1783 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam));
1784 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp));
1785 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp));
1786 PetscFunctionReturn(PETSC_SUCCESS);
1787 }
1788
TSARKIMEXGetVecs(TS ts,DM dm,Vec * Z,Vec * Ydot)1789 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot)
1790 {
1791 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data;
1792
1793 PetscFunctionBegin;
1794 if (Z) {
1795 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z));
1796 else *Z = ax->Z;
1797 }
1798 if (Ydot) {
1799 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot));
1800 else *Ydot = ax->Ydot;
1801 }
1802 PetscFunctionReturn(PETSC_SUCCESS);
1803 }
1804
TSARKIMEXRestoreVecs(TS ts,DM dm,Vec * Z,Vec * Ydot)1805 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot)
1806 {
1807 PetscFunctionBegin;
1808 if (Z) {
1809 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z));
1810 }
1811 if (Ydot) {
1812 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot));
1813 }
1814 PetscFunctionReturn(PETSC_SUCCESS);
1815 }
1816
1817 /*
1818 DAEs need special handling for algebraic variables when restarting DIRK methods with explicit
1819 first stage. In particular, we need:
1820 - to zero the nonlinear function (in case the dual variables are not consistent in the first step)
1821 - to modify the matrix by calling MatZeroRows with identity on these variables.
1822 */
TSARKIMEXComputeAlgebraicIS(TS ts,PetscReal time,Vec X,IS * alg_is)1823 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is)
1824 {
1825 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1826 DM dm;
1827 Vec F, W, Xdot;
1828 const PetscScalar *w;
1829 PetscInt nz = 0, n, st;
1830 PetscInt *nzr;
1831
1832 PetscFunctionBegin;
1833 PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */
1834 PetscCall(DMGetGlobalVector(dm, &Xdot));
1835 PetscCall(DMGetGlobalVector(dm, &F));
1836 PetscCall(DMGetGlobalVector(dm, &W));
1837 PetscCall(VecSet(Xdot, 0.0));
1838 PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex));
1839 PetscCall(VecSetRandom(Xdot, NULL));
1840 PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex));
1841 PetscCall(VecAXPY(W, -1.0, F));
1842 PetscCall(VecGetOwnershipRange(W, &st, NULL));
1843 PetscCall(VecGetLocalSize(W, &n));
1844 PetscCall(VecGetArrayRead(W, &w));
1845 for (PetscInt i = 0; i < n; i++)
1846 if (w[i] == 0.0) nz++;
1847 PetscCall(PetscMalloc1(nz, &nzr));
1848 nz = 0;
1849 for (PetscInt i = 0; i < n; i++)
1850 if (w[i] == 0.0) nzr[nz++] = i + st;
1851 PetscCall(VecRestoreArrayRead(W, &w));
1852 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is));
1853 PetscCall(DMRestoreGlobalVector(dm, &Xdot));
1854 PetscCall(DMRestoreGlobalVector(dm, &F));
1855 PetscCall(DMRestoreGlobalVector(dm, &W));
1856 PetscFunctionReturn(PETSC_SUCCESS);
1857 }
1858
1859 /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure
1860 at the finest level, in the DM for coarser solves. */
TSARKIMEXGetAlgebraicIS(TS ts,DM dm,IS * alg_is)1861 static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is)
1862 {
1863 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data;
1864
1865 PetscFunctionBegin;
1866 if (dm && dm != ts->dm) PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is));
1867 else *alg_is = ax->alg_is;
1868 PetscFunctionReturn(PETSC_SUCCESS);
1869 }
1870
1871 /* This defines the nonlinear equation that is to be solved with SNES */
SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts)1872 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts)
1873 {
1874 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1875 DM dm, dmsave;
1876 Vec Z, Ydot;
1877 IS alg_is;
1878
1879 PetscFunctionBegin;
1880 PetscCall(SNESGetDM(snes, &dm));
1881 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot));
1882 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is));
1883
1884 dmsave = ts->dm;
1885 ts->dm = dm;
1886
1887 if (ark->scoeff == PETSC_MAX_REAL) {
1888 /* We are solving F(t_n,x_n,xdot) = 0 to start the method */
1889 if (!alg_is) {
1890 PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS");
1891 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is));
1892 PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is));
1893 PetscCall(PetscObjectDereference((PetscObject)alg_is));
1894 PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view"));
1895 }
1896 PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex));
1897 PetscCall(VecISSet(F, alg_is, 0.0));
1898 } else {
1899 PetscReal shift = ark->scoeff / ts->time_step;
1900 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */
1901 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex));
1902 }
1903
1904 ts->dm = dmsave;
1905 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot));
1906 PetscFunctionReturn(PETSC_SUCCESS);
1907 }
1908
SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)1909 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts)
1910 {
1911 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1912 DM dm, dmsave;
1913 Vec Ydot, Z;
1914 PetscReal shift;
1915 IS alg_is;
1916
1917 PetscFunctionBegin;
1918 PetscCall(SNESGetDM(snes, &dm));
1919 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */
1920 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot));
1921 /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */
1922 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is));
1923
1924 dmsave = ts->dm;
1925 ts->dm = dm;
1926
1927 if (ark->scoeff == PETSC_MAX_REAL) {
1928 PetscBool hasZeroRows;
1929
1930 /* We are solving F(t_n,x_n,xdot) = 0 to start the method
1931 We compute with a very large shift and then scale back the matrix */
1932 shift = 1.0 / PETSC_MACHINE_EPSILON;
1933 PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex));
1934 PetscCall(MatScale(B, PETSC_MACHINE_EPSILON));
1935 PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows));
1936 if (hasZeroRows) {
1937 PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS");
1938 /* the default of AIJ is to not keep the pattern! We should probably change it someday */
1939 PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
1940 PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL));
1941 }
1942 PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view"));
1943 if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON));
1944 } else {
1945 shift = ark->scoeff / ts->time_step;
1946 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex));
1947 }
1948 ts->dm = dmsave;
1949 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot));
1950 PetscFunctionReturn(PETSC_SUCCESS);
1951 }
1952
TSGetStages_ARKIMEX(TS ts,PetscInt * ns,Vec * Y[])1953 static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[])
1954 {
1955 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1956
1957 PetscFunctionBegin;
1958 if (ns) *ns = ark->tableau->s;
1959 if (Y) *Y = ark->Y;
1960 PetscFunctionReturn(PETSC_SUCCESS);
1961 }
1962
DMCoarsenHook_TSARKIMEX(DM fine,DM coarse,PetscCtx ctx)1963 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, PetscCtx ctx)
1964 {
1965 PetscFunctionBegin;
1966 PetscFunctionReturn(PETSC_SUCCESS);
1967 }
1968
DMRestrictHook_TSARKIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)1969 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
1970 {
1971 TS ts = (TS)ctx;
1972 Vec Z, Z_c;
1973
1974 PetscFunctionBegin;
1975 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL));
1976 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL));
1977 PetscCall(MatRestrict(restrct, Z, Z_c));
1978 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c));
1979 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL));
1980 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL));
1981 PetscFunctionReturn(PETSC_SUCCESS);
1982 }
1983
DMSubDomainHook_TSARKIMEX(DM dm,DM subdm,PetscCtx ctx)1984 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, PetscCtx ctx)
1985 {
1986 PetscFunctionBegin;
1987 PetscFunctionReturn(PETSC_SUCCESS);
1988 }
1989
DMSubDomainRestrictHook_TSARKIMEX(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,PetscCtx ctx)1990 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
1991 {
1992 TS ts = (TS)ctx;
1993 Vec Z, Z_c;
1994
1995 PetscFunctionBegin;
1996 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL));
1997 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL));
1998
1999 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD));
2000 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD));
2001
2002 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL));
2003 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL));
2004 PetscFunctionReturn(PETSC_SUCCESS);
2005 }
2006
TSARKIMEXTableauSetUp(TS ts)2007 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts)
2008 {
2009 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2010 ARKTableau tab = ark->tableau;
2011
2012 PetscFunctionBegin;
2013 PetscCall(PetscMalloc1(2 * tab->s, &ark->work));
2014 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y));
2015 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI));
2016 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS));
2017 if (ark->extrapolate) {
2018 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev));
2019 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev));
2020 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev));
2021 }
2022 PetscFunctionReturn(PETSC_SUCCESS);
2023 }
2024
TSSetUp_ARKIMEX(TS ts)2025 static PetscErrorCode TSSetUp_ARKIMEX(TS ts)
2026 {
2027 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2028 DM dm;
2029 SNES snes;
2030
2031 PetscFunctionBegin;
2032 if (ark->fastslowsplit) {
2033 PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts));
2034 } else {
2035 PetscCall(TSARKIMEXTableauSetUp(ts));
2036 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot));
2037 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0));
2038 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z));
2039 PetscCall(TSGetDM(ts, &dm));
2040 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts));
2041 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts));
2042 PetscCall(TSGetSNES(ts, &snes));
2043 PetscCall(SNESSetDM(snes, dm));
2044 }
2045 PetscFunctionReturn(PETSC_SUCCESS);
2046 }
2047
TSAdjointSetUp_ARKIMEX(TS ts)2048 static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts)
2049 {
2050 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2051 ARKTableau tab = ark->tableau;
2052
2053 PetscFunctionBegin;
2054 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam));
2055 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp));
2056 if (ts->vecs_sensip) PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp));
2057 if (PetscDefined(USE_DEBUG)) {
2058 PetscBool id = PETSC_FALSE;
2059 PetscCall(TSARKIMEXTestMassIdentity(ts, &id));
2060 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Adjoint ARKIMEX requires an identity mass matrix, however the TSIFunctionFn you provided does not utilize an identity mass matrix");
2061 }
2062 PetscFunctionReturn(PETSC_SUCCESS);
2063 }
2064
TSSetFromOptions_ARKIMEX(TS ts,PetscOptionItems PetscOptionsObject)2065 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems PetscOptionsObject)
2066 {
2067 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2068 PetscBool dirk;
2069
2070 PetscFunctionBegin;
2071 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2072 PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options");
2073 {
2074 ARKTableauLink link;
2075 PetscInt count, choice;
2076 PetscBool flg;
2077 const char **namelist;
2078 for (link = ARKTableauList, count = 0; link; link = link->next) {
2079 if (!dirk && link->tab.additive) count++;
2080 if (dirk && !link->tab.additive) count++;
2081 }
2082 PetscCall(PetscMalloc1(count, (char ***)&namelist));
2083 for (link = ARKTableauList, count = 0; link; link = link->next) {
2084 if (!dirk && link->tab.additive) namelist[count++] = link->tab.name;
2085 if (dirk && !link->tab.additive) namelist[count++] = link->tab.name;
2086 }
2087 if (dirk) {
2088 PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg));
2089 if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice]));
2090 } else {
2091 PetscBool fastslowsplit;
2092 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg));
2093 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice]));
2094 flg = (PetscBool)!ark->imex;
2095 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL));
2096 ark->imex = (PetscBool)!flg;
2097 PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg));
2098 if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit));
2099 PetscCall(TSARKIMEXGetFastSlowSplit(ts, &fastslowsplit));
2100 if (fastslowsplit) {
2101 SNES snes;
2102
2103 PetscCall(TSRHSSplitGetSNES(ts, &snes));
2104 PetscCall(SNESSetFromOptions(snes));
2105 }
2106 }
2107 PetscCall(PetscFree(namelist));
2108 PetscCall(PetscOptionsBool("-ts_arkimex_initial_guess_extrapolate", "Extrapolate the initial guess for the stage solution from stage values of the previous time step", "", ark->extrapolate, &ark->extrapolate, NULL));
2109 }
2110 PetscOptionsHeadEnd();
2111 PetscFunctionReturn(PETSC_SUCCESS);
2112 }
2113
TSView_ARKIMEX(TS ts,PetscViewer viewer)2114 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer)
2115 {
2116 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2117 PetscBool isascii, dirk;
2118
2119 PetscFunctionBegin;
2120 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2121 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2122 if (isascii) {
2123 PetscViewerFormat format;
2124 ARKTableau tab = ark->tableau;
2125 TSARKIMEXType arktype;
2126 char buf[2048];
2127 PetscBool flg;
2128
2129 PetscCall(TSARKIMEXGetType(ts, &arktype));
2130 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg));
2131 PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype));
2132 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct));
2133 PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf));
2134 PetscCall(PetscViewerGetFormat(viewer, &format));
2135 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2136 PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff "));
2137 for (PetscInt i = 0; i < tab->s; i++) {
2138 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s));
2139 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf));
2140 }
2141 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt));
2142 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf));
2143 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt));
2144 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf));
2145 }
2146 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no"));
2147 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no"));
2148 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no"));
2149 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no"));
2150 if (!dirk) {
2151 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c));
2152 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf));
2153 }
2154 }
2155 PetscFunctionReturn(PETSC_SUCCESS);
2156 }
2157
TSLoad_ARKIMEX(TS ts,PetscViewer viewer)2158 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer)
2159 {
2160 SNES snes;
2161 TSAdapt adapt;
2162
2163 PetscFunctionBegin;
2164 PetscCall(TSGetAdapt(ts, &adapt));
2165 PetscCall(TSAdaptLoad(adapt, viewer));
2166 PetscCall(TSGetSNES(ts, &snes));
2167 PetscCall(SNESLoad(snes, viewer));
2168 /* function and Jacobian context for SNES when used with TS is always ts object */
2169 PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
2170 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
2171 PetscFunctionReturn(PETSC_SUCCESS);
2172 }
2173
2174 /*@
2175 TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme
2176
2177 Logically Collective
2178
2179 Input Parameters:
2180 + ts - timestepping context
2181 - arktype - type of `TSARKIMEX` scheme
2182
2183 Options Database Key:
2184 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type
2185
2186 Level: intermediate
2187
2188 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`,
2189 `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5`
2190 @*/
TSARKIMEXSetType(TS ts,TSARKIMEXType arktype)2191 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype)
2192 {
2193 PetscFunctionBegin;
2194 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
2195 PetscAssertPointer(arktype, 2);
2196 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype));
2197 PetscFunctionReturn(PETSC_SUCCESS);
2198 }
2199
2200 /*@
2201 TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme
2202
2203 Logically Collective
2204
2205 Input Parameter:
2206 . ts - timestepping context
2207
2208 Output Parameter:
2209 . arktype - type of `TSARKIMEX` scheme
2210
2211 Level: intermediate
2212
2213 .seealso: [](ch_ts), `TSARKIMEX`
2214 @*/
TSARKIMEXGetType(TS ts,TSARKIMEXType * arktype)2215 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype)
2216 {
2217 PetscFunctionBegin;
2218 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
2219 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype));
2220 PetscFunctionReturn(PETSC_SUCCESS);
2221 }
2222
2223 /*@
2224 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly
2225
2226 Logically Collective
2227
2228 Input Parameters:
2229 + ts - timestepping context
2230 - flg - `PETSC_TRUE` for fully implicit
2231
2232 Options Database Key:
2233 . -ts_arkimex_fully_implicit <true,false> - Solve both parts of the equation implicitly
2234
2235 Level: intermediate
2236
2237 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()`
2238 @*/
TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg)2239 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg)
2240 {
2241 PetscFunctionBegin;
2242 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
2243 PetscValidLogicalCollectiveBool(ts, flg, 2);
2244 PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg));
2245 PetscFunctionReturn(PETSC_SUCCESS);
2246 }
2247
2248 /*@
2249 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly
2250
2251 Logically Collective
2252
2253 Input Parameter:
2254 . ts - timestepping context
2255
2256 Output Parameter:
2257 . flg - `PETSC_TRUE` for fully implicit
2258
2259 Level: intermediate
2260
2261 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`
2262 @*/
TSARKIMEXGetFullyImplicit(TS ts,PetscBool * flg)2263 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg)
2264 {
2265 PetscFunctionBegin;
2266 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
2267 PetscAssertPointer(flg, 2);
2268 PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg));
2269 PetscFunctionReturn(PETSC_SUCCESS);
2270 }
2271
TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType * arktype)2272 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype)
2273 {
2274 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2275
2276 PetscFunctionBegin;
2277 *arktype = ark->tableau->name;
2278 PetscFunctionReturn(PETSC_SUCCESS);
2279 }
2280
TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype)2281 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype)
2282 {
2283 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2284 PetscBool match;
2285 ARKTableauLink link;
2286
2287 PetscFunctionBegin;
2288 if (ark->tableau) {
2289 PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match));
2290 if (match) PetscFunctionReturn(PETSC_SUCCESS);
2291 }
2292 for (link = ARKTableauList; link; link = link->next) {
2293 PetscCall(PetscStrcmp(link->tab.name, arktype, &match));
2294 if (match) {
2295 if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts));
2296 ark->tableau = &link->tab;
2297 if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts));
2298 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
2299 PetscFunctionReturn(PETSC_SUCCESS);
2300 }
2301 }
2302 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype);
2303 }
2304
TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg)2305 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg)
2306 {
2307 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2308
2309 PetscFunctionBegin;
2310 ark->imex = (PetscBool)!flg;
2311 PetscFunctionReturn(PETSC_SUCCESS);
2312 }
2313
TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts,PetscBool * flg)2314 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg)
2315 {
2316 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2317
2318 PetscFunctionBegin;
2319 *flg = (PetscBool)!ark->imex;
2320 PetscFunctionReturn(PETSC_SUCCESS);
2321 }
2322
TSDestroy_ARKIMEX(TS ts)2323 static PetscErrorCode TSDestroy_ARKIMEX(TS ts)
2324 {
2325 PetscFunctionBegin;
2326 PetscCall(TSReset_ARKIMEX(ts));
2327 if (ts->dm) {
2328 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts));
2329 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts));
2330 }
2331 PetscCall(PetscFree(ts->data));
2332 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL));
2333 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL));
2334 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL));
2335 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL));
2336 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL));
2337 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL));
2338 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", NULL));
2339 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", NULL));
2340 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", NULL));
2341 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_ARKIMEX_FastSlowSplit_C", NULL));
2342 PetscFunctionReturn(PETSC_SUCCESS);
2343 }
2344
2345 /*MC
2346 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes
2347
2348 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
2349 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
2350 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`.
2351
2352 Options Database Keys:
2353 + -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - Set `TSARKIMEX` scheme type
2354 . -ts_dirk_type <type> - Set `TSDIRK` scheme type
2355 . -ts_arkimex_fully_implicit <true,false> - Solve both parts of the equation implicitly
2356 . -ts_arkimex_fastslowsplit <true,false> - Enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise,
2357 see `TSRHSSplitSetIS()`
2358 - -ts_arkimex_initial_guess_extrapolate - Extrapolate the initial guess for the stage solution from stage values of the previous time step
2359
2360 Level: beginner
2361
2362 Notes:
2363 The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or `-ts_arkimex_type`
2364
2365 If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information.
2366
2367 Methods with an explicit stage can only be used with ODE in which the stiff part $ G(t,X,\dot{X}) $ has the form $ \dot{X} + \hat{G}(t,X)$.
2368
2369 Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear.
2370
2371 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`,
2372 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`,
2373 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType`
2374 M*/
TSCreate_ARKIMEX(TS ts)2375 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts)
2376 {
2377 TS_ARKIMEX *ark;
2378 PetscBool dirk;
2379
2380 PetscFunctionBegin;
2381 PetscCall(TSARKIMEXInitializePackage());
2382 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2383
2384 ts->ops->reset = TSReset_ARKIMEX;
2385 ts->ops->adjointreset = TSAdjointReset_ARKIMEX;
2386 ts->ops->destroy = TSDestroy_ARKIMEX;
2387 ts->ops->view = TSView_ARKIMEX;
2388 ts->ops->load = TSLoad_ARKIMEX;
2389 ts->ops->setup = TSSetUp_ARKIMEX;
2390 ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX;
2391 ts->ops->step = TSStep_ARKIMEX;
2392 ts->ops->interpolate = TSInterpolate_ARKIMEX;
2393 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX;
2394 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX;
2395 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX;
2396 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX;
2397 ts->ops->getstages = TSGetStages_ARKIMEX;
2398 ts->ops->adjointstep = TSAdjointStep_ARKIMEX;
2399
2400 ts->usessnes = PETSC_TRUE;
2401
2402 PetscCall(PetscNew(&ark));
2403 ts->data = (void *)ark;
2404 ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE;
2405
2406 ark->VecsDeltaLam = NULL;
2407 ark->VecsSensiTemp = NULL;
2408 ark->VecsSensiPTemp = NULL;
2409
2410 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX));
2411 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX));
2412 if (!dirk) {
2413 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX));
2414 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX));
2415 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX));
2416 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX));
2417 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault));
2418 }
2419 PetscFunctionReturn(PETSC_SUCCESS);
2420 }
2421
TSDIRKSetType_DIRK(TS ts,TSDIRKType dirktype)2422 static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype)
2423 {
2424 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2425
2426 PetscFunctionBegin;
2427 PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype));
2428 PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype);
2429 PetscFunctionReturn(PETSC_SUCCESS);
2430 }
2431
2432 /*@
2433 TSDIRKSetType - Set the type of `TSDIRK` scheme
2434
2435 Logically Collective
2436
2437 Input Parameters:
2438 + ts - timestepping context
2439 - dirktype - type of `TSDIRK` scheme
2440
2441 Options Database Key:
2442 . -ts_dirkimex_type - set `TSDIRK` scheme type
2443
2444 Level: intermediate
2445
2446 .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType`
2447 @*/
TSDIRKSetType(TS ts,TSDIRKType dirktype)2448 PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype)
2449 {
2450 PetscFunctionBegin;
2451 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
2452 PetscAssertPointer(dirktype, 2);
2453 PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype));
2454 PetscFunctionReturn(PETSC_SUCCESS);
2455 }
2456
2457 /*@
2458 TSDIRKGetType - Get the type of `TSDIRK` scheme
2459
2460 Logically Collective
2461
2462 Input Parameter:
2463 . ts - timestepping context
2464
2465 Output Parameter:
2466 . dirktype - type of `TSDIRK` scheme
2467
2468 Level: intermediate
2469
2470 .seealso: [](ch_ts), `TSDIRKSetType()`
2471 @*/
TSDIRKGetType(TS ts,TSDIRKType * dirktype)2472 PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype)
2473 {
2474 PetscFunctionBegin;
2475 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
2476 PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype));
2477 PetscFunctionReturn(PETSC_SUCCESS);
2478 }
2479
2480 /*MC
2481 TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes.
2482
2483 Level: beginner
2484
2485 Notes:
2486 The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or `-ts_dirk_type`.
2487 The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with:
2488 + E - whether the method has an explicit first stage
2489 . S - whether the method is single diagonal
2490 . P - order of the advancing method
2491 . Q - order of the embedded method
2492 . S - number of stages
2493 . SA - whether the method is stiffly accurate
2494 . L - whether the method is L-stable
2495 - A - whether the method is A-stable
2496
2497 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`.
2498 M*/
TSCreate_DIRK(TS ts)2499 PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts)
2500 {
2501 PetscFunctionBegin;
2502 PetscCall(TSCreate_ARKIMEX(ts));
2503 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX));
2504 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK));
2505 PetscCall(TSDIRKSetType(ts, TSDIRKDefault));
2506 PetscFunctionReturn(PETSC_SUCCESS);
2507 }
2508
2509 /*@
2510 TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system
2511
2512 Logically Collective
2513
2514 Input Parameters:
2515 + ts - timestepping context
2516 - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise.
2517
2518 Options Database Key:
2519 . -ts_arkimex_fastslowsplit <true,false> - enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise
2520
2521 Level: intermediate
2522
2523 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()`, `TSRHSSplitSetIS()`
2524 @*/
TSARKIMEXSetFastSlowSplit(TS ts,PetscBool fastslow)2525 PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow)
2526 {
2527 PetscFunctionBegin;
2528 PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow));
2529 PetscFunctionReturn(PETSC_SUCCESS);
2530 }
2531
2532 /*@
2533 TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system
2534
2535 Not Collective
2536
2537 Input Parameter:
2538 . ts - timestepping context
2539
2540 Output Parameter:
2541 . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise
2542
2543 Level: intermediate
2544
2545 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()`
2546 @*/
TSARKIMEXGetFastSlowSplit(TS ts,PetscBool * fastslow)2547 PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow)
2548 {
2549 PetscFunctionBegin;
2550 PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow));
2551 PetscFunctionReturn(PETSC_SUCCESS);
2552 }
2553