Actual source code: arkimex.c
1: /*
2: Code for timestepping with additive Runge-Kutta IMEX method or Diagonally Implicit Runge-Kutta methods.
4: Notes:
5: For ARK, the general system is written as
7: F(t,U,Udot) = G(t,U)
9: where F represents the stiff part of the physics and G represents the non-stiff part.
11: */
12: #include <petsc/private/tsimpl.h>
13: #include <petscdm.h>
14: #include <../src/ts/impls/arkimex/arkimex.h>
15: #include <../src/ts/impls/arkimex/fsarkimex.h>
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);
24: /*MC
25: TSARKIMEXARS122 - Second order ARK IMEX scheme, {cite}`ascher_1997`
27: This method has one explicit stage and one implicit stage.
29: Options Database Key:
30: . -ts_arkimex_type ars122 - set arkimex type to ars122
32: Level: advanced
34: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
35: M*/
37: /*MC
38: TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part.
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.
42: Options Database Key:
43: . -ts_arkimex_type a2 - set arkimex type to a2
45: Level: advanced
47: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
48: M*/
50: /*MC
51: TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part, {cite}`pareschi_2005`
53: This method has two implicit stages, and L-stable implicit scheme.
55: Options Database Key:
56: . -ts_arkimex_type l2 - set arkimex type to l2
58: Level: advanced
60: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
61: M*/
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.
66: This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used.
68: Options Database Key:
69: . -ts_arkimex_type 1bee - set arkimex type to 1bee
71: Level: advanced
73: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
74: M*/
76: /*MC
77: TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part.
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.
81: Options Database Key:
82: . -ts_arkimex_type 2c - set arkimex type to 2c
84: Level: advanced
86: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
87: M*/
89: /*MC
90: TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part.
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.
94: Options Database Key:
95: . -ts_arkimex_type 2d - set arkimex type to 2d
97: Level: advanced
99: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
100: M*/
102: /*MC
103: TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part.
105: This method has one explicit stage and two implicit stages. It is an optimal method developed by Emil Constantinescu.
107: Options Database Key:
108: . -ts_arkimex_type 2e - set arkimex type to 2e
110: Level: advanced
112: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
113: M*/
115: /*MC
116: TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme, {cite}`pareschi_2005`
118: This method has three implicit stages.
120: This method is referred to as SSP2-(3,3,2) in <https://arxiv.org/abs/1110.4375>
122: Options Database Key:
123: . -ts_arkimex_type prssp2 - set arkimex type to prssp2
125: Level: advanced
127: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
128: M*/
130: /*MC
131: TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`
133: This method has one explicit stage and three implicit stages.
135: Options Database Key:
136: . -ts_arkimex_type 3 - set arkimex type to 3
138: Level: advanced
140: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
141: M*/
143: /*MC
144: TSARKIMEXARS443 - Third order ARK IMEX scheme, {cite}`ascher_1997`
146: This method has one explicit stage and four implicit stages.
148: Options Database Key:
149: . -ts_arkimex_type ars443 - set arkimex type to ars443
151: Level: advanced
153: Notes:
154: This method is referred to as ARS(4,4,3) in <https://arxiv.org/abs/1110.4375>
156: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
157: M*/
159: /*MC
160: TSARKIMEXBPR3 - Third order ARK IMEX scheme. Referred to as ARK3 in <https://arxiv.org/abs/1110.4375>
162: This method has one explicit stage and four implicit stages.
164: Options Database Key:
165: . -ts_arkimex_type bpr3 - set arkimex type to bpr3
167: Level: advanced
169: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
170: M*/
172: /*MC
173: TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`.
175: This method has one explicit stage and four implicit stages.
177: Options Database Key:
178: . -ts_arkimex_type 4 - set arkimex type to4
180: Level: advanced
182: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
183: M*/
185: /*MC
186: TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`.
188: This method has one explicit stage and five implicit stages.
190: Options Database Key:
191: . -ts_arkimex_type 5 - set arkimex type to 5
193: Level: advanced
195: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
196: M*/
198: /*MC
199: TSDIRKS212 - Second order DIRK scheme.
201: This method has two implicit stages with an embedded method of other 1.
202: See `TSDIRK` for additional details.
204: Options Database Key:
205: . -ts_dirk_type s212 - select this method.
207: Level: advanced
209: Note:
210: This is the default DIRK scheme in SUNDIALS.
212: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
213: M*/
215: /*MC
216: TSDIRKES122SAL - First order DIRK scheme <https://arxiv.org/abs/1803.01613>
218: Uses backward Euler as advancing method and trapezoidal rule as embedded method. See `TSDIRK` for additional details.
220: Options Database Key:
221: . -ts_dirk_type es122sal - select this method.
223: Level: advanced
225: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
226: M*/
228: /*MC
229: TSDIRKES213SAL - Second order DIRK scheme {cite}`kennedy2019diagonally`. Also known as TR-BDF2, see{cite}`hosea1996analysis`
231: See `TSDIRK` for additional details.
233: Options Database Key:
234: . -ts_dirk_type es213sal - select this method.
236: Level: advanced
238: Note:
239: This is the default DIRK scheme used in PETSc.
241: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
242: M*/
244: /*MC
245: TSDIRKES324SAL - Third order DIRK scheme, {cite}`kennedy2019diagonally`
247: See `TSDIRK` for additional details.
249: Options Database Key:
250: . -ts_dirk_type es324sal - select this method.
252: Level: advanced
254: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
255: M*/
257: /*MC
258: TSDIRKES325SAL - Third order DIRK scheme {cite}`kennedy2019diagonally`.
260: See `TSDIRK` for additional details.
262: Options Database Key:
263: . -ts_dirk_type es325sal - select this method.
265: Level: advanced
267: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
268: M*/
270: /*MC
271: TSDIRK657A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
273: See `TSDIRK` for additional details.
275: Options Database Key:
276: . -ts_dirk_type 657a - select this method.
278: Level: advanced
280: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
281: M*/
283: /*MC
284: TSDIRKES648SA - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
286: See `TSDIRK` for additional details.
288: Options Database Key:
289: . -ts_dirk_type es648sa - select this method.
291: Level: advanced
293: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
294: M*/
296: /*MC
297: TSDIRK658A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
299: See `TSDIRK` for additional details.
301: Options Database Key:
302: . -ts_dirk_type 658a - select this method.
304: Level: advanced
306: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
307: M*/
309: /*MC
310: TSDIRKS659A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
312: See `TSDIRK` for additional details.
314: Options Database Key:
315: . -ts_dirk_type s659a - select this method.
317: Level: advanced
319: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
320: M*/
322: /*MC
323: TSDIRK7510SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
325: See `TSDIRK` for additional details.
327: Options Database Key:
328: . -ts_dirk_type 7510sal - select this method.
330: Level: advanced
332: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
333: M*/
335: /*MC
336: TSDIRKES7510SA - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
338: See `TSDIRK` for additional details.
340: Options Database Key:
341: . -ts_dirk_type es7510sa - select this method.
343: Level: advanced
345: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
346: M*/
348: /*MC
349: TSDIRK759A - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
351: See `TSDIRK` for additional details.
353: Options Database Key:
354: . -ts_dirk_type 759a - select this method.
356: Level: advanced
358: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
359: M*/
361: /*MC
362: TSDIRKS7511SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
364: See `TSDIRK` for additional details.
366: Options Database Key:
367: . -ts_dirk_type s7511sal - select this method.
369: Level: advanced
371: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
372: M*/
374: /*MC
375: TSDIRK8614A - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
377: See `TSDIRK` for additional details.
379: Options Database Key:
380: . -ts_dirk_type 8614a - select this method.
382: Level: advanced
384: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
385: M*/
387: /*MC
388: TSDIRK8616SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
390: See `TSDIRK` for additional details.
392: Options Database Key:
393: . -ts_dirk_type 8616sal - select this method.
395: Level: advanced
397: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
398: M*/
400: /*MC
401: TSDIRKES8516SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs>
403: See `TSDIRK` for additional details.
405: Options Database Key:
406: . -ts_dirk_type es8516sal - select this method.
408: Level: advanced
410: .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()`
411: M*/
413: PetscErrorCode TSHasRHSFunction(TS ts, PetscBool *has)
414: {
415: TSRHSFunctionFn *func;
417: PetscFunctionBegin;
418: *has = PETSC_FALSE;
419: PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL));
420: if (func) *has = PETSC_TRUE;
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@C
425: TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX`
427: Not Collective, but should be called by all processes which will need the schemes to be registered
429: Level: advanced
431: .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()`
432: @*/
433: PetscErrorCode TSARKIMEXRegisterAll(void)
434: {
435: PetscFunctionBegin;
436: if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
437: TSARKIMEXRegisterAllCalled = PETSC_TRUE;
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); */
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: }
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: }
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: }
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: }
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: }
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: };
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: }
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: }
1015: /*@C
1016: TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`.
1018: Not Collective
1020: Level: advanced
1022: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()`
1023: @*/
1024: PetscErrorCode TSARKIMEXRegisterDestroy(void)
1025: {
1026: ARKTableauLink link;
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: }
1042: /*@C
1043: TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called
1044: from `TSInitializePackage()`.
1046: Level: developer
1048: .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()`
1049: @*/
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: }
1060: /*@C
1061: TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is
1062: called from `PetscFinalize()`.
1064: Level: developer
1066: .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()`
1067: @*/
1068: PetscErrorCode TSARKIMEXFinalizePackage(void)
1069: {
1070: PetscFunctionBegin;
1071: TSARKIMEXPackageInitialized = PETSC_FALSE;
1072: PetscCall(TSARKIMEXRegisterDestroy());
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }
1076: /*@C
1077: TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
1079: Logically Collective
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)
1097: Level: advanced
1099: Note:
1100: Several `TSARKIMEX` methods are provided, this function is only needed to create new methods.
1102: .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS`
1103: @*/
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;
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;
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: }
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];
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));
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];
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));
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;
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;
1161: /* def of FSAL can be made more precise */
1162: t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate);
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: }
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));
1175: link->next = ARKTableauList;
1176: ARKTableauList = link;
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: /*@C
1181: TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation
1183: Logically Collective.
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)
1196: Level: advanced
1198: Note:
1199: Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods.
1201: .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS`
1202: @*/
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: }
1210: /*
1211: The step completion formula is
1213: x1 = x0 - h bt^T YdotI + h b^T YdotRHS
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
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
1223: so we can evaluate the method with different order even after the step has been optimistically completed.
1224: */
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;
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: if (done) *done = PETSC_FALSE;
1295: else
1296: SETERRQ(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.", tab->name,
1297: tab->order, order);
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1301: static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id)
1302: {
1303: Vec Udot, Y1, Y2;
1304: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1305: PetscReal norm;
1307: PetscFunctionBegin;
1308: PetscCall(VecDuplicate(ts->vec_sol, &Udot));
1309: PetscCall(VecDuplicate(ts->vec_sol, &Y1));
1310: PetscCall(VecDuplicate(ts->vec_sol, &Y2));
1311: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex));
1312: PetscCall(VecSetRandom(Udot, NULL));
1313: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex));
1314: PetscCall(VecAXPY(Y2, -1.0, Y1));
1315: PetscCall(VecAXPY(Y2, -1.0, Udot));
1316: PetscCall(VecNorm(Y2, NORM_2, &norm));
1317: if (norm < 100.0 * PETSC_MACHINE_EPSILON) {
1318: *id = PETSC_TRUE;
1319: } else {
1320: *id = PETSC_FALSE;
1321: 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));
1322: }
1323: PetscCall(VecDestroy(&Udot));
1324: PetscCall(VecDestroy(&Y1));
1325: PetscCall(VecDestroy(&Y2));
1326: PetscFunctionReturn(PETSC_SUCCESS);
1327: }
1329: static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS, PetscReal, Vec, IS *);
1331: static PetscErrorCode TSStep_ARKIMEX(TS ts)
1332: {
1333: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1334: ARKTableau tab = ark->tableau;
1335: const PetscInt s = tab->s;
1336: const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c;
1337: PetscScalar *w = ark->work;
1338: Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z;
1339: PetscBool extrapolate = ark->extrapolate;
1340: TSAdapt adapt;
1341: SNES snes;
1342: PetscInt i, j, its, lits;
1343: PetscInt rejections = 0;
1344: PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE;
1345: PetscReal next_time_step = ts->time_step;
1347: PetscFunctionBegin;
1348: if (ark->extrapolate && !ark->Y_prev) {
1349: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev));
1350: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev));
1351: if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev));
1352: }
1354: if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE));
1355: if (!hasE) dirk = PETSC_TRUE;
1357: if (!ts->steprollback) {
1358: if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */
1359: PetscCall(VecCopy(YdotI[s - 1], Ydot0));
1360: }
1361: if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */
1362: for (i = 0; i < s; i++) {
1363: PetscCall(VecCopy(Y[i], ark->Y_prev[i]));
1364: PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i]));
1365: if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i]));
1366: }
1367: }
1368: }
1370: /*
1371: For fully implicit formulations we must solve the equations
1373: F(t_n,x_n,xdot) = 0
1375: for the explicit first stage.
1376: Here we call SNESSolve using PETSC_MAX_REAL as shift to flag it.
1377: Special handling is inside SNESTSFormFunction_ARKIMEX and SNESTSFormJacobian_ARKIMEX
1378: We compute Ydot0 if we restart the step or if we resized the problem after remeshing
1379: */
1380: if (dirk && tab->explicit_first_stage && (ts->steprestart || ts->stepresize)) {
1381: ark->scoeff = PETSC_MAX_REAL;
1382: PetscCall(VecCopy(ts->vec_sol, Z));
1383: if (!ark->alg_is) {
1384: PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ts->ptime, Z, &ark->alg_is));
1385: PetscCall(ISViewFromOptions(ark->alg_is, (PetscObject)ts, "-ts_arkimex_algebraic_is_view"));
1386: }
1387: PetscCall(TSGetSNES(ts, &snes));
1388: PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, 1));
1389: PetscCall(SNESSolve(snes, NULL, Ydot0));
1390: if (ark->alg_is) PetscCall(VecISSet(Ydot0, ark->alg_is, 0.0));
1391: PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, -1));
1392: }
1394: /* For IMEX we compute a step */
1395: if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) {
1396: TS ts_start;
1397: if (PetscDefined(USE_DEBUG) && hasE) {
1398: PetscBool id = PETSC_FALSE;
1399: PetscCall(TSARKIMEXTestMassIdentity(ts, &id));
1400: 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");
1401: }
1402: PetscCall(TSClone(ts, &ts_start));
1403: PetscCall(TSSetSolution(ts_start, ts->vec_sol));
1404: PetscCall(TSSetTime(ts_start, ts->ptime));
1405: PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1));
1406: PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step));
1407: PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER));
1408: PetscCall(TSSetTimeStep(ts_start, ts->time_step));
1409: PetscCall(TSSetType(ts_start, TSARKIMEX));
1410: PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE));
1411: PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE));
1413: PetscCall(TSRestartStep(ts_start));
1414: PetscCall(TSSolve(ts_start, ts->vec_sol));
1415: PetscCall(TSGetTime(ts_start, &ts->ptime));
1416: PetscCall(TSGetTimeStep(ts_start, &ts->time_step));
1418: { /* Save the initial slope for the next step */
1419: TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data;
1420: PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0));
1421: }
1422: ts->steps++;
1424: /* Set the correct TS in SNES */
1425: /* We'll try to bypass this by changing the method on the fly */
1426: {
1427: PetscCall(TSGetSNES(ts, &snes));
1428: PetscCall(TSSetSNES(ts, snes));
1429: }
1430: PetscCall(TSDestroy(&ts_start));
1431: }
1433: ark->status = TS_STEP_INCOMPLETE;
1434: while (!ts->reason && ark->status != TS_STEP_COMPLETE) {
1435: PetscReal t = ts->ptime;
1436: PetscReal h = ts->time_step;
1437: for (i = 0; i < s; i++) {
1438: ark->stage_time = t + h * ct[i];
1439: PetscCall(TSPreStage(ts, ark->stage_time));
1440: if (At[i * s + i] == 0) { /* This stage is explicit */
1441: 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");
1442: PetscCall(VecCopy(ts->vec_sol, Y[i]));
1443: for (j = 0; j < i; j++) w[j] = h * At[i * s + j];
1444: PetscCall(VecMAXPY(Y[i], i, w, YdotI));
1445: if (tab->additive && hasE) {
1446: for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
1447: PetscCall(VecMAXPY(Y[i], i, w, YdotRHS));
1448: }
1449: PetscCall(TSGetSNES(ts, &snes));
1450: PetscCall(SNESResetCounters(snes));
1451: } else {
1452: ark->scoeff = 1. / At[i * s + i];
1453: /* Ydot = shift*(Y-Z) */
1454: PetscCall(VecCopy(ts->vec_sol, Z));
1455: for (j = 0; j < i; j++) w[j] = h * At[i * s + j];
1456: PetscCall(VecMAXPY(Z, i, w, YdotI));
1457: if (tab->additive && hasE) {
1458: for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
1459: PetscCall(VecMAXPY(Z, i, w, YdotRHS));
1460: }
1461: if (extrapolate && !ts->steprestart) {
1462: /* Initial guess extrapolated from previous time step stage values */
1463: PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i]));
1464: } else {
1465: /* Initial guess taken from last stage */
1466: PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i]));
1467: }
1468: PetscCall(TSGetSNES(ts, &snes));
1469: PetscCall(SNESSolve(snes, NULL, Y[i]));
1470: PetscCall(SNESGetIterationNumber(snes, &its));
1471: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1472: ts->snes_its += its;
1473: ts->ksp_its += lits;
1474: PetscCall(TSGetAdapt(ts, &adapt));
1475: PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok));
1476: if (!stageok) {
1477: /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to
1478: * use extrapolation to initialize the solves on the next attempt. */
1479: extrapolate = PETSC_FALSE;
1480: goto reject_step;
1481: }
1482: }
1483: if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) {
1484: if (i == 0 && tab->explicit_first_stage) {
1485: 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",
1486: ((PetscObject)ts)->type_name, ark->tableau->name);
1487: PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */
1488: } else {
1489: PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */
1490: }
1491: } else {
1492: if (i == 0 && tab->explicit_first_stage) {
1493: PetscCall(VecZeroEntries(Ydot));
1494: PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */
1495: PetscCall(VecScale(YdotI[i], -1.0));
1496: } else {
1497: PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */
1498: }
1499: if (hasE) {
1500: if (ark->imex) {
1501: PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i]));
1502: } else {
1503: PetscCall(VecZeroEntries(YdotRHS[i]));
1504: }
1505: }
1506: }
1507: PetscCall(TSPostStage(ts, ark->stage_time, i, Y));
1508: }
1510: ark->status = TS_STEP_INCOMPLETE;
1511: PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL));
1512: ark->status = TS_STEP_PENDING;
1513: PetscCall(TSGetAdapt(ts, &adapt));
1514: PetscCall(TSAdaptCandidatesClear(adapt));
1515: PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
1516: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
1517: ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1518: if (!accept) { /* Roll back the current step */
1519: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
1520: ts->time_step = next_time_step;
1521: goto reject_step;
1522: }
1524: ts->ptime += ts->time_step;
1525: ts->time_step = next_time_step;
1526: break;
1528: reject_step:
1529: ts->reject++;
1530: accept = PETSC_FALSE;
1531: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1532: ts->reason = TS_DIVERGED_STEP_REJECTED;
1533: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
1534: }
1535: }
1536: PetscFunctionReturn(PETSC_SUCCESS);
1537: }
1539: /*
1540: This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form
1541: Udot = H(t,U) + G(t,U)
1542: This corresponds to F(t,U,Udot) = Udot-H(t,U).
1544: The complete adjoint equations are
1545: (shift*I - dHdu) lambda_s[i] = 1/at[i][i] (
1546: dGdU (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j])
1547: + dHdU (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0
1548: lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j]
1549: mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i]
1550: + dGdP (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j])
1551: + dHdP (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0
1552: mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j]
1553: where shift = 1/(at[i][i]*h)
1555: If at[i][i] is 0, the first equation falls back to
1556: lambda_s[i] = h (
1557: (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j]
1558: + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j])
1560: */
1561: static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts)
1562: {
1563: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1564: ARKTableau tab = ark->tableau;
1565: const PetscInt s = tab->s;
1566: const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt;
1567: PetscScalar *w = ark->work;
1568: Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp;
1569: Mat Jex, Jim, Jimpre;
1570: PetscInt i, j, nadj;
1571: PetscReal t = ts->ptime, stage_time_ex;
1572: PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */
1573: KSP ksp;
1575: PetscFunctionBegin;
1576: ark->status = TS_STEP_INCOMPLETE;
1577: PetscCall(SNESGetKSP(ts->snes, &ksp));
1578: PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL));
1579: PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL));
1581: for (i = s - 1; i >= 0; i--) {
1582: ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]);
1583: stage_time_ex = t - adjoint_time_step * (1.0 - c[i]);
1584: if (At[i * s + i] == 0) { // This stage is explicit
1585: ark->scoeff = 0.;
1586: } else {
1587: ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative
1588: }
1589: PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre));
1590: PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex));
1591: if (ts->vecs_sensip) {
1592: 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
1593: PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP
1594: }
1595: /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */
1596: for (nadj = 0; nadj < ts->numcost; nadj++) {
1597: /* build implicit part */
1598: PetscCall(VecSet(VecsSensiTemp[nadj], 0));
1599: if (s - i - 1 > 0) {
1600: /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */
1601: for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i];
1602: PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
1603: }
1604: /* Temp = Temp - bt[i] lambda_{n+1} */
1605: PetscCall(VecAXPY(VecsSensiTemp[nadj], -bt[i], ts->vecs_sensi[nadj]));
1606: if (bt[i] || s - i - 1 > 0) {
1607: /* (shift I - dHdU) Temp */
1608: PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i]));
1609: /* cancel out shift Temp where shift=-scoeff/h */
1610: PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj]));
1611: if (ts->vecs_sensip) {
1612: /* - dHdP Temp */
1613: PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj]));
1614: /* mu_n += -h dHdP Temp */
1615: PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj]));
1616: }
1617: } else {
1618: PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized
1619: }
1620: /* build explicit part */
1621: PetscCall(VecSet(VecsSensiTemp[nadj], 0));
1622: if (s - i - 1 > 0) {
1623: /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */
1624: for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i];
1625: PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
1626: }
1627: /* Temp = Temp + b[i] lambda_{n+1} */
1628: PetscCall(VecAXPY(VecsSensiTemp[nadj], b[i], ts->vecs_sensi[nadj]));
1629: if (b[i] || s - i - 1 > 0) {
1630: /* dGdU Temp */
1631: PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i]));
1632: if (ts->vecs_sensip) {
1633: /* dGdP Temp */
1634: PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj]));
1635: /* mu_n += h dGdP Temp */
1636: PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj]));
1637: }
1638: }
1639: /* Build LHS for first-order adjoint */
1640: if (At[i * s + i] == 0) { // This stage is explicit
1641: PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step));
1642: } else {
1643: KSP ksp;
1644: KSPConvergedReason kspreason;
1645: PetscCall(SNESGetKSP(ts->snes, &ksp));
1646: PetscCall(KSPSetOperators(ksp, Jim, Jimpre));
1647: PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i]));
1648: PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i]));
1649: PetscCall(KSPGetConvergedReason(ksp, &kspreason));
1650: if (kspreason < 0) {
1651: ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
1652: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj));
1653: }
1654: if (ts->vecs_sensip) {
1655: /* -dHdP lambda_s[i] */
1656: PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj]));
1657: /* mu_n += h at[i][i] dHdP lambda_s[i] */
1658: PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj]));
1659: }
1660: }
1661: }
1662: }
1663: for (j = 0; j < s; j++) w[j] = 1.0;
1664: for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's
1665: PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s]));
1666: ark->status = TS_STEP_COMPLETE;
1667: PetscFunctionReturn(PETSC_SUCCESS);
1668: }
1670: static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X)
1671: {
1672: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1673: ARKTableau tab = ark->tableau;
1674: PetscInt s = tab->s, pinterp = tab->pinterp, i, j;
1675: PetscReal h;
1676: PetscReal tt, t;
1677: PetscScalar *bt = ark->work, *b = ark->work + s;
1678: const PetscReal *Bt = tab->binterpt, *B = tab->binterp;
1680: PetscFunctionBegin;
1681: PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name);
1682: switch (ark->status) {
1683: case TS_STEP_INCOMPLETE:
1684: case TS_STEP_PENDING:
1685: h = ts->time_step;
1686: t = (itime - ts->ptime) / h;
1687: break;
1688: case TS_STEP_COMPLETE:
1689: h = ts->ptime - ts->ptime_prev;
1690: t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1691: break;
1692: default:
1693: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1694: }
1695: for (i = 0; i < s; i++) bt[i] = b[i] = 0;
1696: for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1697: for (i = 0; i < s; i++) {
1698: bt[i] += h * Bt[i * pinterp + j] * tt;
1699: b[i] += h * B[i * pinterp + j] * tt;
1700: }
1701: }
1702: PetscCall(VecCopy(ark->Y[0], X));
1703: PetscCall(VecMAXPY(X, s, bt, ark->YdotI));
1704: if (tab->additive) {
1705: PetscBool hasE;
1706: PetscCall(TSHasRHSFunction(ts, &hasE));
1707: if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS));
1708: }
1709: PetscFunctionReturn(PETSC_SUCCESS);
1710: }
1712: static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X)
1713: {
1714: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1715: ARKTableau tab = ark->tableau;
1716: PetscInt s = tab->s, pinterp = tab->pinterp, i, j;
1717: PetscReal h, h_prev, t, tt;
1718: PetscScalar *bt = ark->work, *b = ark->work + s;
1719: const PetscReal *Bt = tab->binterpt, *B = tab->binterp;
1721: PetscFunctionBegin;
1722: PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name);
1723: h = ts->time_step;
1724: h_prev = ts->ptime - ts->ptime_prev;
1725: t = 1 + h / h_prev * c;
1726: for (i = 0; i < s; i++) bt[i] = b[i] = 0;
1727: for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1728: for (i = 0; i < s; i++) {
1729: bt[i] += h * Bt[i * pinterp + j] * tt;
1730: b[i] += h * B[i * pinterp + j] * tt;
1731: }
1732: }
1733: PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored");
1734: PetscCall(VecCopy(ark->Y_prev[0], X));
1735: PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev));
1736: if (tab->additive) {
1737: PetscBool hasE;
1738: PetscCall(TSHasRHSFunction(ts, &hasE));
1739: if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev));
1740: }
1741: PetscFunctionReturn(PETSC_SUCCESS);
1742: }
1744: static PetscErrorCode TSARKIMEXTableauReset(TS ts)
1745: {
1746: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1747: ARKTableau tab = ark->tableau;
1749: PetscFunctionBegin;
1750: if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1751: PetscCall(PetscFree(ark->work));
1752: PetscCall(VecDestroyVecs(tab->s, &ark->Y));
1753: PetscCall(VecDestroyVecs(tab->s, &ark->YdotI));
1754: PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS));
1755: PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev));
1756: PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev));
1757: PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev));
1758: PetscFunctionReturn(PETSC_SUCCESS);
1759: }
1761: static PetscErrorCode TSReset_ARKIMEX(TS ts)
1762: {
1763: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1765: PetscFunctionBegin;
1766: if (ark->fastslowsplit) {
1767: PetscTryMethod(ts, "TSReset_ARKIMEX_FastSlowSplit_C", (TS), (ts));
1768: } else {
1769: PetscCall(TSARKIMEXTableauReset(ts));
1770: PetscCall(VecDestroy(&ark->Ydot));
1771: PetscCall(VecDestroy(&ark->Ydot0));
1772: PetscCall(VecDestroy(&ark->Z));
1773: PetscCall(ISDestroy(&ark->alg_is));
1774: }
1775: PetscFunctionReturn(PETSC_SUCCESS);
1776: }
1778: static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts)
1779: {
1780: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1781: ARKTableau tab = ark->tableau;
1783: PetscFunctionBegin;
1784: PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam));
1785: PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp));
1786: PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp));
1787: PetscFunctionReturn(PETSC_SUCCESS);
1788: }
1790: static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot)
1791: {
1792: TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data;
1794: PetscFunctionBegin;
1795: if (Z) {
1796: if (dm && dm != ts->dm) {
1797: PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z));
1798: } else *Z = ax->Z;
1799: }
1800: if (Ydot) {
1801: if (dm && dm != ts->dm) {
1802: PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot));
1803: } else *Ydot = ax->Ydot;
1804: }
1805: PetscFunctionReturn(PETSC_SUCCESS);
1806: }
1808: static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot)
1809: {
1810: PetscFunctionBegin;
1811: if (Z) {
1812: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z));
1813: }
1814: if (Ydot) {
1815: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot));
1816: }
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: /*
1821: DAEs need special handling for algebraic variables when restarting DIRK methods with explicit
1822: first stage. In particular, we need:
1823: - to zero the nonlinear function (in case the dual variables are not consistent in the first step)
1824: - to modify the preconditioning matrix by calling MatZeroRows with identity on these variables.
1825: */
1826: static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is)
1827: {
1828: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1829: DM dm;
1830: Vec F, W, Xdot;
1831: const PetscScalar *w;
1832: PetscInt nz = 0, n, st;
1833: PetscInt *nzr;
1835: PetscFunctionBegin;
1836: PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */
1837: PetscCall(DMGetGlobalVector(dm, &Xdot));
1838: PetscCall(DMGetGlobalVector(dm, &F));
1839: PetscCall(DMGetGlobalVector(dm, &W));
1840: PetscCall(VecSet(Xdot, 0.0));
1841: PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex));
1842: PetscCall(VecSetRandom(Xdot, NULL));
1843: PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex));
1844: PetscCall(VecAXPY(W, -1.0, F));
1845: PetscCall(VecGetOwnershipRange(W, &st, NULL));
1846: PetscCall(VecGetLocalSize(W, &n));
1847: PetscCall(VecGetArrayRead(W, &w));
1848: for (PetscInt i = 0; i < n; i++)
1849: if (w[i] == 0.0) nz++;
1850: PetscCall(PetscMalloc1(nz, &nzr));
1851: nz = 0;
1852: for (PetscInt i = 0; i < n; i++)
1853: if (w[i] == 0.0) nzr[nz++] = i + st;
1854: PetscCall(VecRestoreArrayRead(W, &w));
1855: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is));
1856: PetscCall(DMRestoreGlobalVector(dm, &Xdot));
1857: PetscCall(DMRestoreGlobalVector(dm, &F));
1858: PetscCall(DMRestoreGlobalVector(dm, &W));
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1862: /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure
1863: at the finest level, in the DM for coarser solves. */
1864: static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is)
1865: {
1866: TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data;
1868: PetscFunctionBegin;
1869: if (dm && dm != ts->dm) {
1870: PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is));
1871: } else *alg_is = ax->alg_is;
1872: PetscFunctionReturn(PETSC_SUCCESS);
1873: }
1875: /* This defines the nonlinear equation that is to be solved with SNES */
1876: static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts)
1877: {
1878: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1879: DM dm, dmsave;
1880: Vec Z, Ydot;
1881: IS alg_is;
1883: PetscFunctionBegin;
1884: PetscCall(SNESGetDM(snes, &dm));
1885: PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot));
1886: if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is));
1888: dmsave = ts->dm;
1889: ts->dm = dm;
1891: if (ark->scoeff == PETSC_MAX_REAL) {
1892: /* We are solving F(t_n,x_n,xdot) = 0 to start the method */
1893: if (!alg_is) {
1894: PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS");
1895: PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is));
1896: PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is));
1897: PetscCall(PetscObjectDereference((PetscObject)alg_is));
1898: PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view"));
1899: }
1900: PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex));
1901: PetscCall(VecISSet(F, alg_is, 0.0));
1902: } else {
1903: PetscReal shift = ark->scoeff / ts->time_step;
1904: PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */
1905: PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex));
1906: }
1908: ts->dm = dmsave;
1909: PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot));
1910: PetscFunctionReturn(PETSC_SUCCESS);
1911: }
1913: static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts)
1914: {
1915: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1916: DM dm, dmsave;
1917: Vec Ydot, Z;
1918: PetscReal shift;
1919: IS alg_is;
1921: PetscFunctionBegin;
1922: PetscCall(SNESGetDM(snes, &dm));
1923: /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */
1924: PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot));
1925: /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */
1926: if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is));
1928: dmsave = ts->dm;
1929: ts->dm = dm;
1931: if (ark->scoeff == PETSC_MAX_REAL) {
1932: PetscBool hasZeroRows;
1934: /* We are solving F(t_n,x_n,xdot) = 0 to start the method
1935: We compute with a very large shift and then scale back the matrix */
1936: shift = 1.0 / PETSC_MACHINE_EPSILON;
1937: PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex));
1938: PetscCall(MatScale(B, PETSC_MACHINE_EPSILON));
1939: PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows));
1940: if (hasZeroRows) {
1941: PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS");
1942: /* the default of AIJ is to not keep the pattern! We should probably change it someday */
1943: PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
1944: PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL));
1945: }
1946: PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view"));
1947: if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON));
1948: } else {
1949: shift = ark->scoeff / ts->time_step;
1950: PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex));
1951: }
1952: ts->dm = dmsave;
1953: PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot));
1954: PetscFunctionReturn(PETSC_SUCCESS);
1955: }
1957: static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[])
1958: {
1959: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1961: PetscFunctionBegin;
1962: if (ns) *ns = ark->tableau->s;
1963: if (Y) *Y = ark->Y;
1964: PetscFunctionReturn(PETSC_SUCCESS);
1965: }
1967: static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx)
1968: {
1969: PetscFunctionBegin;
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
1974: {
1975: TS ts = (TS)ctx;
1976: Vec Z, Z_c;
1978: PetscFunctionBegin;
1979: PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL));
1980: PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL));
1981: PetscCall(MatRestrict(restrct, Z, Z_c));
1982: PetscCall(VecPointwiseMult(Z_c, rscale, Z_c));
1983: PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL));
1984: PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL));
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx)
1989: {
1990: PetscFunctionBegin;
1991: PetscFunctionReturn(PETSC_SUCCESS);
1992: }
1994: static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
1995: {
1996: TS ts = (TS)ctx;
1997: Vec Z, Z_c;
1999: PetscFunctionBegin;
2000: PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL));
2001: PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL));
2003: PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD));
2004: PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD));
2006: PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL));
2007: PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL));
2008: PetscFunctionReturn(PETSC_SUCCESS);
2009: }
2011: static PetscErrorCode TSARKIMEXTableauSetUp(TS ts)
2012: {
2013: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2014: ARKTableau tab = ark->tableau;
2016: PetscFunctionBegin;
2017: PetscCall(PetscMalloc1(2 * tab->s, &ark->work));
2018: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y));
2019: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI));
2020: if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS));
2021: if (ark->extrapolate) {
2022: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev));
2023: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev));
2024: if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev));
2025: }
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: static PetscErrorCode TSSetUp_ARKIMEX(TS ts)
2030: {
2031: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2032: DM dm;
2033: SNES snes;
2035: PetscFunctionBegin;
2036: if (ark->fastslowsplit) {
2037: PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts));
2038: } else {
2039: PetscCall(TSARKIMEXTableauSetUp(ts));
2040: PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot));
2041: PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0));
2042: PetscCall(VecDuplicate(ts->vec_sol, &ark->Z));
2043: PetscCall(TSGetDM(ts, &dm));
2044: PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts));
2045: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts));
2046: PetscCall(TSGetSNES(ts, &snes));
2047: PetscCall(SNESSetDM(snes, dm));
2048: }
2049: PetscFunctionReturn(PETSC_SUCCESS);
2050: }
2052: static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts)
2053: {
2054: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2055: ARKTableau tab = ark->tableau;
2057: PetscFunctionBegin;
2058: PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam));
2059: PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp));
2060: if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); }
2061: if (PetscDefined(USE_DEBUG)) {
2062: PetscBool id = PETSC_FALSE;
2063: PetscCall(TSARKIMEXTestMassIdentity(ts, &id));
2064: 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");
2065: }
2066: PetscFunctionReturn(PETSC_SUCCESS);
2067: }
2069: static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject)
2070: {
2071: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2072: PetscBool dirk;
2074: PetscFunctionBegin;
2075: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2076: PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options");
2077: {
2078: ARKTableauLink link;
2079: PetscInt count, choice;
2080: PetscBool flg;
2081: const char **namelist;
2082: for (link = ARKTableauList, count = 0; link; link = link->next) {
2083: if (!dirk && link->tab.additive) count++;
2084: if (dirk && !link->tab.additive) count++;
2085: }
2086: PetscCall(PetscMalloc1(count, (char ***)&namelist));
2087: for (link = ARKTableauList, count = 0; link; link = link->next) {
2088: if (!dirk && link->tab.additive) namelist[count++] = link->tab.name;
2089: if (dirk && !link->tab.additive) namelist[count++] = link->tab.name;
2090: }
2091: if (dirk) {
2092: PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg));
2093: if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice]));
2094: } else {
2095: PetscBool fastslowsplit;
2096: PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg));
2097: if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice]));
2098: flg = (PetscBool)!ark->imex;
2099: PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL));
2100: ark->imex = (PetscBool)!flg;
2101: PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg));
2102: if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit));
2103: PetscCall(TSARKIMEXGetFastSlowSplit(ts, &fastslowsplit));
2104: if (fastslowsplit) {
2105: SNES snes;
2107: PetscCall(TSRHSSplitGetSNES(ts, &snes));
2108: PetscCall(SNESSetFromOptions(snes));
2109: }
2110: }
2111: PetscCall(PetscFree(namelist));
2112: 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));
2113: }
2114: PetscOptionsHeadEnd();
2115: PetscFunctionReturn(PETSC_SUCCESS);
2116: }
2118: static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer)
2119: {
2120: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2121: PetscBool iascii, dirk;
2123: PetscFunctionBegin;
2124: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2125: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2126: if (iascii) {
2127: PetscViewerFormat format;
2128: ARKTableau tab = ark->tableau;
2129: TSARKIMEXType arktype;
2130: char buf[2048];
2131: PetscBool flg;
2133: PetscCall(TSARKIMEXGetType(ts, &arktype));
2134: PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg));
2135: PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype));
2136: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct));
2137: PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf));
2138: PetscCall(PetscViewerGetFormat(viewer, &format));
2139: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2140: PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff "));
2141: for (PetscInt i = 0; i < tab->s; i++) {
2142: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s));
2143: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf));
2144: }
2145: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt));
2146: PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf));
2147: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt));
2148: PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf));
2149: }
2150: PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no"));
2151: PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no"));
2152: PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no"));
2153: PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no"));
2154: if (!dirk) {
2155: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c));
2156: PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf));
2157: }
2158: }
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer)
2163: {
2164: SNES snes;
2165: TSAdapt adapt;
2167: PetscFunctionBegin;
2168: PetscCall(TSGetAdapt(ts, &adapt));
2169: PetscCall(TSAdaptLoad(adapt, viewer));
2170: PetscCall(TSGetSNES(ts, &snes));
2171: PetscCall(SNESLoad(snes, viewer));
2172: /* function and Jacobian context for SNES when used with TS is always ts object */
2173: PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
2174: PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
2175: PetscFunctionReturn(PETSC_SUCCESS);
2176: }
2178: /*@
2179: TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme
2181: Logically Collective
2183: Input Parameters:
2184: + ts - timestepping context
2185: - arktype - type of `TSARKIMEX` scheme
2187: Options Database Key:
2188: . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type
2190: Level: intermediate
2192: .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`,
2193: `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5`
2194: @*/
2195: PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype)
2196: {
2197: PetscFunctionBegin;
2199: PetscAssertPointer(arktype, 2);
2200: PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype));
2201: PetscFunctionReturn(PETSC_SUCCESS);
2202: }
2204: /*@
2205: TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme
2207: Logically Collective
2209: Input Parameter:
2210: . ts - timestepping context
2212: Output Parameter:
2213: . arktype - type of `TSARKIMEX` scheme
2215: Level: intermediate
2217: .seealso: [](ch_ts), `TSARKIMEXc`
2218: @*/
2219: PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype)
2220: {
2221: PetscFunctionBegin;
2223: PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype));
2224: PetscFunctionReturn(PETSC_SUCCESS);
2225: }
2227: /*@
2228: TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly
2230: Logically Collective
2232: Input Parameters:
2233: + ts - timestepping context
2234: - flg - `PETSC_TRUE` for fully implicit
2236: Level: intermediate
2238: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()`
2239: @*/
2240: PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg)
2241: {
2242: PetscFunctionBegin;
2245: PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg));
2246: PetscFunctionReturn(PETSC_SUCCESS);
2247: }
2249: /*@
2250: TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly
2252: Logically Collective
2254: Input Parameter:
2255: . ts - timestepping context
2257: Output Parameter:
2258: . flg - `PETSC_TRUE` for fully implicit
2260: Level: intermediate
2262: .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`
2263: @*/
2264: PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg)
2265: {
2266: PetscFunctionBegin;
2268: PetscAssertPointer(flg, 2);
2269: PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg));
2270: PetscFunctionReturn(PETSC_SUCCESS);
2271: }
2273: static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype)
2274: {
2275: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2277: PetscFunctionBegin;
2278: *arktype = ark->tableau->name;
2279: PetscFunctionReturn(PETSC_SUCCESS);
2280: }
2282: static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype)
2283: {
2284: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2285: PetscBool match;
2286: ARKTableauLink link;
2288: PetscFunctionBegin;
2289: if (ark->tableau) {
2290: PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match));
2291: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2292: }
2293: for (link = ARKTableauList; link; link = link->next) {
2294: PetscCall(PetscStrcmp(link->tab.name, arktype, &match));
2295: if (match) {
2296: if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts));
2297: ark->tableau = &link->tab;
2298: if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts));
2299: ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2302: }
2303: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype);
2304: }
2306: static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg)
2307: {
2308: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2310: PetscFunctionBegin;
2311: ark->imex = (PetscBool)!flg;
2312: PetscFunctionReturn(PETSC_SUCCESS);
2313: }
2315: static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg)
2316: {
2317: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2319: PetscFunctionBegin;
2320: *flg = (PetscBool)!ark->imex;
2321: PetscFunctionReturn(PETSC_SUCCESS);
2322: }
2324: static PetscErrorCode TSDestroy_ARKIMEX(TS ts)
2325: {
2326: PetscFunctionBegin;
2327: PetscCall(TSReset_ARKIMEX(ts));
2328: if (ts->dm) {
2329: PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts));
2330: PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts));
2331: }
2332: PetscCall(PetscFree(ts->data));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL));
2335: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL));
2336: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL));
2337: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL));
2338: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL));
2339: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", NULL));
2340: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", NULL));
2341: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", NULL));
2342: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_ARKIMEX_FastSlowSplit_C", NULL));
2343: PetscFunctionReturn(PETSC_SUCCESS);
2344: }
2346: /* ------------------------------------------------------------ */
2347: /*MC
2348: TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes
2350: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
2351: nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
2352: of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`.
2354: Level: beginner
2356: Notes:
2357: The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type
2359: If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information.
2361: Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).
2363: Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear.
2365: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`,
2366: `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`,
2367: `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType`
2368: M*/
2369: PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts)
2370: {
2371: TS_ARKIMEX *ark;
2372: PetscBool dirk;
2374: PetscFunctionBegin;
2375: PetscCall(TSARKIMEXInitializePackage());
2376: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk));
2378: ts->ops->reset = TSReset_ARKIMEX;
2379: ts->ops->adjointreset = TSAdjointReset_ARKIMEX;
2380: ts->ops->destroy = TSDestroy_ARKIMEX;
2381: ts->ops->view = TSView_ARKIMEX;
2382: ts->ops->load = TSLoad_ARKIMEX;
2383: ts->ops->setup = TSSetUp_ARKIMEX;
2384: ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX;
2385: ts->ops->step = TSStep_ARKIMEX;
2386: ts->ops->interpolate = TSInterpolate_ARKIMEX;
2387: ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX;
2388: ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX;
2389: ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX;
2390: ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX;
2391: ts->ops->getstages = TSGetStages_ARKIMEX;
2392: ts->ops->adjointstep = TSAdjointStep_ARKIMEX;
2394: ts->usessnes = PETSC_TRUE;
2396: PetscCall(PetscNew(&ark));
2397: ts->data = (void *)ark;
2398: ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE;
2400: ark->VecsDeltaLam = NULL;
2401: ark->VecsSensiTemp = NULL;
2402: ark->VecsSensiPTemp = NULL;
2404: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX));
2405: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX));
2406: if (!dirk) {
2407: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX));
2408: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX));
2409: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX));
2410: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX));
2411: PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault));
2412: }
2413: PetscFunctionReturn(PETSC_SUCCESS);
2414: }
2416: /* ------------------------------------------------------------ */
2418: static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype)
2419: {
2420: TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
2422: PetscFunctionBegin;
2423: PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype));
2424: PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype);
2425: PetscFunctionReturn(PETSC_SUCCESS);
2426: }
2428: /*@
2429: TSDIRKSetType - Set the type of `TSDIRK` scheme
2431: Logically Collective
2433: Input Parameters:
2434: + ts - timestepping context
2435: - dirktype - type of `TSDIRK` scheme
2437: Options Database Key:
2438: . -ts_dirkimex_type - set `TSDIRK` scheme type
2440: Level: intermediate
2442: .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType`
2443: @*/
2444: PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype)
2445: {
2446: PetscFunctionBegin;
2448: PetscAssertPointer(dirktype, 2);
2449: PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype));
2450: PetscFunctionReturn(PETSC_SUCCESS);
2451: }
2453: /*@
2454: TSDIRKGetType - Get the type of `TSDIRK` scheme
2456: Logically Collective
2458: Input Parameter:
2459: . ts - timestepping context
2461: Output Parameter:
2462: . dirktype - type of `TSDIRK` scheme
2464: Level: intermediate
2466: .seealso: [](ch_ts), `TSDIRKSetType()`
2467: @*/
2468: PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype)
2469: {
2470: PetscFunctionBegin;
2472: PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype));
2473: PetscFunctionReturn(PETSC_SUCCESS);
2474: }
2476: /*MC
2477: TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes.
2479: Level: beginner
2481: Notes:
2482: The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type.
2483: The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with:
2484: + E - whether the method has an explicit first stage
2485: . S - whether the method is single diagonal
2486: . P - order of the advancing method
2487: . Q - order of the embedded method
2488: . S - number of stages
2489: . SA - whether the method is stiffly accurate
2490: . L - whether the method is L-stable
2491: - A - whether the method is A-stable
2493: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`.
2494: M*/
2495: PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts)
2496: {
2497: PetscFunctionBegin;
2498: PetscCall(TSCreate_ARKIMEX(ts));
2499: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX));
2500: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK));
2501: PetscCall(TSDIRKSetType(ts, TSDIRKDefault));
2502: PetscFunctionReturn(PETSC_SUCCESS);
2503: }
2505: /*@
2506: TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system
2508: Logically Collective
2510: Input Parameters:
2511: + ts - timestepping context
2512: - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise.
2514: Options Database Key:
2515: . -ts_arkimex_fastslowsplit - <true,false>
2517: Level: intermediate
2519: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()`
2520: @*/
2521: PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow)
2522: {
2523: PetscFunctionBegin;
2524: PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow));
2525: PetscFunctionReturn(PETSC_SUCCESS);
2526: }
2528: /*@
2529: TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system
2531: Not Collective
2533: Input Parameter:
2534: . ts - timestepping context
2536: Output Parameter:
2537: . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise
2539: Level: intermediate
2541: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()`
2542: @*/
2543: PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow)
2544: {
2545: PetscFunctionBegin;
2546: PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow));
2547: PetscFunctionReturn(PETSC_SUCCESS);
2548: }