18a381b04SJed Brown /* 2d27334e2SStefano Zampini Code for timestepping with additive Runge-Kutta IMEX method or Diagonally Implicit Runge-Kutta methods. 38a381b04SJed Brown 48a381b04SJed Brown Notes: 5d27334e2SStefano Zampini For ARK, the general system is written as 68a381b04SJed Brown 7f9c1d6abSBarry Smith F(t,U,Udot) = G(t,U) 88a381b04SJed Brown 98a381b04SJed Brown where F represents the stiff part of the physics and G represents the non-stiff part. 108a381b04SJed Brown 118a381b04SJed Brown */ 12af0996ceSBarry Smith #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 131e25c274SJed Brown #include <petscdm.h> 143a2a065fSHong Zhang #include <../src/ts/impls/arkimex/arkimex.h> 153a2a065fSHong Zhang #include <../src/ts/impls/arkimex/fsarkimex.h> 168a381b04SJed Brown 173a2a065fSHong Zhang static ARKTableauLink ARKTableauList; 1819fd82e9SBarry Smith static TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3; 193405e92cSStefano Zampini static TSDIRKType TSDIRKDefault = TSDIRKES213SAL; 208a381b04SJed Brown static PetscBool TSARKIMEXRegisterAllCalled; 218a381b04SJed Brown static PetscBool TSARKIMEXPackageInitialized; 2256dcabbaSDebojyoti Ghosh static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec); 238a381b04SJed Brown 241f80e275SEmil Constantinescu /*MC 251d27aa22SBarry Smith TSARKIMEXARS122 - Second order ARK IMEX scheme, {cite}`ascher_1997` 268a381b04SJed Brown 271f80e275SEmil Constantinescu This method has one explicit stage and one implicit stage. 281f80e275SEmil Constantinescu 29bcf0153eSBarry Smith Options Database Key: 3067b8a455SSatish Balay . -ts_arkimex_type ars122 - set arkimex type to ars122 319bd3e852SBarry Smith 32bcf0153eSBarry Smith Level: advanced 33bcf0153eSBarry Smith 341cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 351f80e275SEmil Constantinescu M*/ 36d27334e2SStefano Zampini 371f80e275SEmil Constantinescu /*MC 381f80e275SEmil Constantinescu TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part. 391f80e275SEmil Constantinescu 401f80e275SEmil Constantinescu This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu. 411f80e275SEmil Constantinescu 42bcf0153eSBarry Smith Options Database Key: 4367b8a455SSatish Balay . -ts_arkimex_type a2 - set arkimex type to a2 449bd3e852SBarry Smith 451f80e275SEmil Constantinescu Level: advanced 461f80e275SEmil Constantinescu 471cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 481f80e275SEmil Constantinescu M*/ 49d27334e2SStefano Zampini 501f80e275SEmil Constantinescu /*MC 511d27aa22SBarry Smith TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part, {cite}`pareschi_2005` 521f80e275SEmil Constantinescu 531f80e275SEmil Constantinescu This method has two implicit stages, and L-stable implicit scheme. 541f80e275SEmil Constantinescu 55bcf0153eSBarry Smith Options Database Key: 5667b8a455SSatish Balay . -ts_arkimex_type l2 - set arkimex type to l2 579bd3e852SBarry Smith 58bcf0153eSBarry Smith Level: advanced 59bcf0153eSBarry Smith 601cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 611f80e275SEmil Constantinescu M*/ 62d27334e2SStefano Zampini 631f80e275SEmil Constantinescu /*MC 64c79dcfbdSBarry Smith TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method. 65e817cc15SEmil Constantinescu 66e817cc15SEmil Constantinescu This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used. 67e817cc15SEmil Constantinescu 68bcf0153eSBarry Smith Options Database Key: 6967b8a455SSatish Balay . -ts_arkimex_type 1bee - set arkimex type to 1bee 709bd3e852SBarry Smith 71e817cc15SEmil Constantinescu Level: advanced 72e817cc15SEmil Constantinescu 731cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 74e817cc15SEmil Constantinescu M*/ 75d27334e2SStefano Zampini 76e817cc15SEmil Constantinescu /*MC 771f80e275SEmil Constantinescu TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part. 781f80e275SEmil Constantinescu 791f80e275SEmil Constantinescu 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. 801f80e275SEmil Constantinescu 81bcf0153eSBarry Smith Options Database Key: 8267b8a455SSatish Balay . -ts_arkimex_type 2c - set arkimex type to 2c 839bd3e852SBarry Smith 841f80e275SEmil Constantinescu Level: advanced 851f80e275SEmil Constantinescu 861cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 871f80e275SEmil Constantinescu M*/ 88d27334e2SStefano Zampini 8964f491ddSJed Brown /*MC 9064f491ddSJed Brown TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 9164f491ddSJed Brown 92da81f932SPierre Jolivet 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. 9364f491ddSJed Brown 94bcf0153eSBarry Smith Options Database Key: 9567b8a455SSatish Balay . -ts_arkimex_type 2d - set arkimex type to 2d 969bd3e852SBarry Smith 97b330ce4dSSatish Balay Level: advanced 98b330ce4dSSatish Balay 991cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 10064f491ddSJed Brown M*/ 101d27334e2SStefano Zampini 10264f491ddSJed Brown /*MC 10364f491ddSJed Brown TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 10464f491ddSJed Brown 10515229ffcSPierre Jolivet This method has one explicit stage and two implicit stages. It is an optimal method developed by Emil Constantinescu. 10664f491ddSJed Brown 107bcf0153eSBarry Smith Options Database Key: 10867b8a455SSatish Balay . -ts_arkimex_type 2e - set arkimex type to 2e 1099bd3e852SBarry Smith 110b330ce4dSSatish Balay Level: advanced 111b330ce4dSSatish Balay 1121cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 11364f491ddSJed Brown M*/ 114d27334e2SStefano Zampini 11564f491ddSJed Brown /*MC 1161d27aa22SBarry Smith TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme, {cite}`pareschi_2005` 1176cf0794eSJed Brown 1186cf0794eSJed Brown This method has three implicit stages. 1196cf0794eSJed Brown 1201d27aa22SBarry Smith This method is referred to as SSP2-(3,3,2) in <https://arxiv.org/abs/1110.4375> 1216cf0794eSJed Brown 122bcf0153eSBarry Smith Options Database Key: 12367b8a455SSatish Balay . -ts_arkimex_type prssp2 - set arkimex type to prssp2 1249bd3e852SBarry Smith 1256cf0794eSJed Brown Level: advanced 1266cf0794eSJed Brown 1271cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 1286cf0794eSJed Brown M*/ 129d27334e2SStefano Zampini 1306cf0794eSJed Brown /*MC 1311d27aa22SBarry Smith TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003` 13264f491ddSJed Brown 13364f491ddSJed Brown This method has one explicit stage and three implicit stages. 13464f491ddSJed Brown 135bcf0153eSBarry Smith Options Database Key: 13667b8a455SSatish Balay . -ts_arkimex_type 3 - set arkimex type to 3 1379bd3e852SBarry Smith 138bcf0153eSBarry Smith Level: advanced 139bcf0153eSBarry Smith 1401cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 14164f491ddSJed Brown M*/ 142d27334e2SStefano Zampini 14364f491ddSJed Brown /*MC 1441d27aa22SBarry Smith TSARKIMEXARS443 - Third order ARK IMEX scheme, {cite}`ascher_1997` 1456cf0794eSJed Brown 1466cf0794eSJed Brown This method has one explicit stage and four implicit stages. 1476cf0794eSJed Brown 148bcf0153eSBarry Smith Options Database Key: 14967b8a455SSatish Balay . -ts_arkimex_type ars443 - set arkimex type to ars443 1509bd3e852SBarry Smith 151bcf0153eSBarry Smith Level: advanced 152bcf0153eSBarry Smith 1531d27aa22SBarry Smith Notes: 1541d27aa22SBarry Smith This method is referred to as ARS(4,4,3) in <https://arxiv.org/abs/1110.4375> 1556cf0794eSJed Brown 1561cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 1576cf0794eSJed Brown M*/ 158d27334e2SStefano Zampini 1596cf0794eSJed Brown /*MC 1601d27aa22SBarry Smith TSARKIMEXBPR3 - Third order ARK IMEX scheme. Referred to as ARK3 in <https://arxiv.org/abs/1110.4375> 1616cf0794eSJed Brown 1626cf0794eSJed Brown This method has one explicit stage and four implicit stages. 1636cf0794eSJed Brown 164bcf0153eSBarry Smith Options Database Key: 16567b8a455SSatish Balay . -ts_arkimex_type bpr3 - set arkimex type to bpr3 1669bd3e852SBarry Smith 167bcf0153eSBarry Smith Level: advanced 168bcf0153eSBarry Smith 1691cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 1706cf0794eSJed Brown M*/ 171d27334e2SStefano Zampini 1726cf0794eSJed Brown /*MC 1731d27aa22SBarry Smith TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`. 17464f491ddSJed Brown 17564f491ddSJed Brown This method has one explicit stage and four implicit stages. 17664f491ddSJed Brown 177bcf0153eSBarry Smith Options Database Key: 17867b8a455SSatish Balay . -ts_arkimex_type 4 - set arkimex type to4 1799bd3e852SBarry Smith 180bcf0153eSBarry Smith Level: advanced 181bcf0153eSBarry Smith 1821cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 18364f491ddSJed Brown M*/ 184d27334e2SStefano Zampini 18564f491ddSJed Brown /*MC 1861d27aa22SBarry Smith TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`. 18764f491ddSJed Brown 18864f491ddSJed Brown This method has one explicit stage and five implicit stages. 18964f491ddSJed Brown 190bcf0153eSBarry Smith Options Database Key: 19167b8a455SSatish Balay . -ts_arkimex_type 5 - set arkimex type to 5 1929bd3e852SBarry Smith 193bcf0153eSBarry Smith Level: advanced 194bcf0153eSBarry Smith 1951cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 19664f491ddSJed Brown M*/ 19764f491ddSJed Brown 1983405e92cSStefano Zampini /*MC 1993405e92cSStefano Zampini TSDIRKS212 - Second order DIRK scheme. 2003405e92cSStefano Zampini 2013405e92cSStefano Zampini This method has two implicit stages with an embedded method of other 1. 2023405e92cSStefano Zampini See `TSDIRK` for additional details. 2033405e92cSStefano Zampini 2043405e92cSStefano Zampini Options Database Key: 2053405e92cSStefano Zampini . -ts_dirk_type s212 - select this method. 2063405e92cSStefano Zampini 2073405e92cSStefano Zampini Level: advanced 2083405e92cSStefano Zampini 2093405e92cSStefano Zampini Note: 2103405e92cSStefano Zampini This is the default DIRK scheme in SUNDIALS. 2113405e92cSStefano Zampini 2123405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2133405e92cSStefano Zampini M*/ 2143405e92cSStefano Zampini 2153405e92cSStefano Zampini /*MC 2161d27aa22SBarry Smith TSDIRKES122SAL - First order DIRK scheme <https://arxiv.org/abs/1803.01613> 2173405e92cSStefano Zampini 2183405e92cSStefano Zampini Uses backward Euler as advancing method and trapezoidal rule as embedded method. See `TSDIRK` for additional details. 2193405e92cSStefano Zampini 2203405e92cSStefano Zampini Options Database Key: 2213405e92cSStefano Zampini . -ts_dirk_type es122sal - select this method. 2223405e92cSStefano Zampini 2233405e92cSStefano Zampini Level: advanced 2243405e92cSStefano Zampini 2253405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2263405e92cSStefano Zampini M*/ 2273405e92cSStefano Zampini 2283405e92cSStefano Zampini /*MC 2291d27aa22SBarry Smith TSDIRKES213SAL - Second order DIRK scheme {cite}`kennedy2019diagonally`. Also known as TR-BDF2, see{cite}`hosea1996analysis` 2303405e92cSStefano Zampini 2313405e92cSStefano Zampini See `TSDIRK` for additional details. 2323405e92cSStefano Zampini 2333405e92cSStefano Zampini Options Database Key: 2343405e92cSStefano Zampini . -ts_dirk_type es213sal - select this method. 2353405e92cSStefano Zampini 2363405e92cSStefano Zampini Level: advanced 2373405e92cSStefano Zampini 2383405e92cSStefano Zampini Note: 2393405e92cSStefano Zampini This is the default DIRK scheme used in PETSc. 2403405e92cSStefano Zampini 2413405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2423405e92cSStefano Zampini M*/ 2433405e92cSStefano Zampini 2443405e92cSStefano Zampini /*MC 2451d27aa22SBarry Smith TSDIRKES324SAL - Third order DIRK scheme, {cite}`kennedy2019diagonally` 2463405e92cSStefano Zampini 2473405e92cSStefano Zampini See `TSDIRK` for additional details. 2483405e92cSStefano Zampini 2493405e92cSStefano Zampini Options Database Key: 2503405e92cSStefano Zampini . -ts_dirk_type es324sal - select this method. 2513405e92cSStefano Zampini 2523405e92cSStefano Zampini Level: advanced 2533405e92cSStefano Zampini 2543405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2553405e92cSStefano Zampini M*/ 2563405e92cSStefano Zampini 2573405e92cSStefano Zampini /*MC 2581d27aa22SBarry Smith TSDIRKES325SAL - Third order DIRK scheme {cite}`kennedy2019diagonally`. 2593405e92cSStefano Zampini 2603405e92cSStefano Zampini See `TSDIRK` for additional details. 2613405e92cSStefano Zampini 2623405e92cSStefano Zampini Options Database Key: 2633405e92cSStefano Zampini . -ts_dirk_type es325sal - select this method. 2643405e92cSStefano Zampini 2653405e92cSStefano Zampini Level: advanced 2663405e92cSStefano Zampini 2673405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2683405e92cSStefano Zampini M*/ 2693405e92cSStefano Zampini 2703405e92cSStefano Zampini /*MC 2711d27aa22SBarry Smith TSDIRK657A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 2723405e92cSStefano Zampini 2733405e92cSStefano Zampini See `TSDIRK` for additional details. 2743405e92cSStefano Zampini 2753405e92cSStefano Zampini Options Database Key: 2763405e92cSStefano Zampini . -ts_dirk_type 657a - select this method. 2773405e92cSStefano Zampini 2783405e92cSStefano Zampini Level: advanced 2793405e92cSStefano Zampini 2803405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2813405e92cSStefano Zampini M*/ 2823405e92cSStefano Zampini 2833405e92cSStefano Zampini /*MC 2841d27aa22SBarry Smith TSDIRKES648SA - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 2853405e92cSStefano Zampini 2863405e92cSStefano Zampini See `TSDIRK` for additional details. 2873405e92cSStefano Zampini 2883405e92cSStefano Zampini Options Database Key: 2893405e92cSStefano Zampini . -ts_dirk_type es648sa - select this method. 2903405e92cSStefano Zampini 2913405e92cSStefano Zampini Level: advanced 2923405e92cSStefano Zampini 2933405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 2943405e92cSStefano Zampini M*/ 2953405e92cSStefano Zampini 2963405e92cSStefano Zampini /*MC 2971d27aa22SBarry Smith TSDIRK658A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 2983405e92cSStefano Zampini 2993405e92cSStefano Zampini See `TSDIRK` for additional details. 3003405e92cSStefano Zampini 3013405e92cSStefano Zampini Options Database Key: 3023405e92cSStefano Zampini . -ts_dirk_type 658a - select this method. 3033405e92cSStefano Zampini 3043405e92cSStefano Zampini Level: advanced 3053405e92cSStefano Zampini 3063405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3073405e92cSStefano Zampini M*/ 3083405e92cSStefano Zampini 3093405e92cSStefano Zampini /*MC 3101d27aa22SBarry Smith TSDIRKS659A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3113405e92cSStefano Zampini 3123405e92cSStefano Zampini See `TSDIRK` for additional details. 3133405e92cSStefano Zampini 3143405e92cSStefano Zampini Options Database Key: 3153405e92cSStefano Zampini . -ts_dirk_type s659a - select this method. 3163405e92cSStefano Zampini 3173405e92cSStefano Zampini Level: advanced 3183405e92cSStefano Zampini 3193405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3203405e92cSStefano Zampini M*/ 3213405e92cSStefano Zampini 3223405e92cSStefano Zampini /*MC 3231d27aa22SBarry Smith TSDIRK7510SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3243405e92cSStefano Zampini 3253405e92cSStefano Zampini See `TSDIRK` for additional details. 3263405e92cSStefano Zampini 3273405e92cSStefano Zampini Options Database Key: 3283405e92cSStefano Zampini . -ts_dirk_type 7510sal - select this method. 3293405e92cSStefano Zampini 3303405e92cSStefano Zampini Level: advanced 3313405e92cSStefano Zampini 3323405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3333405e92cSStefano Zampini M*/ 3343405e92cSStefano Zampini 3353405e92cSStefano Zampini /*MC 3361d27aa22SBarry Smith TSDIRKES7510SA - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3373405e92cSStefano Zampini 3383405e92cSStefano Zampini See `TSDIRK` for additional details. 3393405e92cSStefano Zampini 3403405e92cSStefano Zampini Options Database Key: 3413405e92cSStefano Zampini . -ts_dirk_type es7510sa - select this method. 3423405e92cSStefano Zampini 3433405e92cSStefano Zampini Level: advanced 3443405e92cSStefano Zampini 3453405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3463405e92cSStefano Zampini M*/ 3473405e92cSStefano Zampini 3483405e92cSStefano Zampini /*MC 3491d27aa22SBarry Smith TSDIRK759A - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3503405e92cSStefano Zampini 3513405e92cSStefano Zampini See `TSDIRK` for additional details. 3523405e92cSStefano Zampini 3533405e92cSStefano Zampini Options Database Key: 3543405e92cSStefano Zampini . -ts_dirk_type 759a - select this method. 3553405e92cSStefano Zampini 3563405e92cSStefano Zampini Level: advanced 3573405e92cSStefano Zampini 3583405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3593405e92cSStefano Zampini M*/ 3603405e92cSStefano Zampini 3613405e92cSStefano Zampini /*MC 3621d27aa22SBarry Smith TSDIRKS7511SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3633405e92cSStefano Zampini 3643405e92cSStefano Zampini See `TSDIRK` for additional details. 3653405e92cSStefano Zampini 3663405e92cSStefano Zampini Options Database Key: 3673405e92cSStefano Zampini . -ts_dirk_type s7511sal - select this method. 3683405e92cSStefano Zampini 3693405e92cSStefano Zampini Level: advanced 3703405e92cSStefano Zampini 3713405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3723405e92cSStefano Zampini M*/ 3733405e92cSStefano Zampini 3743405e92cSStefano Zampini /*MC 3751d27aa22SBarry Smith TSDIRK8614A - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3763405e92cSStefano Zampini 3773405e92cSStefano Zampini See `TSDIRK` for additional details. 3783405e92cSStefano Zampini 3793405e92cSStefano Zampini Options Database Key: 3803405e92cSStefano Zampini . -ts_dirk_type 8614a - select this method. 3813405e92cSStefano Zampini 3823405e92cSStefano Zampini Level: advanced 3833405e92cSStefano Zampini 3843405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3853405e92cSStefano Zampini M*/ 3863405e92cSStefano Zampini 3873405e92cSStefano Zampini /*MC 3881d27aa22SBarry Smith TSDIRK8616SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 3893405e92cSStefano Zampini 3903405e92cSStefano Zampini See `TSDIRK` for additional details. 3913405e92cSStefano Zampini 3923405e92cSStefano Zampini Options Database Key: 3933405e92cSStefano Zampini . -ts_dirk_type 8616sal - select this method. 3943405e92cSStefano Zampini 3953405e92cSStefano Zampini Level: advanced 3963405e92cSStefano Zampini 3973405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 3983405e92cSStefano Zampini M*/ 3993405e92cSStefano Zampini 4003405e92cSStefano Zampini /*MC 4011d27aa22SBarry Smith TSDIRKES8516SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 4023405e92cSStefano Zampini 4033405e92cSStefano Zampini See `TSDIRK` for additional details. 4043405e92cSStefano Zampini 4053405e92cSStefano Zampini Options Database Key: 4063405e92cSStefano Zampini . -ts_dirk_type es8516sal - select this method. 4073405e92cSStefano Zampini 4083405e92cSStefano Zampini Level: advanced 4093405e92cSStefano Zampini 4103405e92cSStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 4113405e92cSStefano Zampini M*/ 4123405e92cSStefano Zampini 4133a2a065fSHong Zhang PetscErrorCode TSHasRHSFunction(TS ts, PetscBool *has) 414d27334e2SStefano Zampini { 4158434afd1SBarry Smith TSRHSFunctionFn *func; 416d27334e2SStefano Zampini 417d27334e2SStefano Zampini PetscFunctionBegin; 418d27334e2SStefano Zampini *has = PETSC_FALSE; 419d27334e2SStefano Zampini PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL)); 420d27334e2SStefano Zampini if (func) *has = PETSC_TRUE; 421d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 422d27334e2SStefano Zampini } 423d27334e2SStefano Zampini 4248a381b04SJed Brown /*@C 425bcf0153eSBarry Smith TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX` 4268a381b04SJed Brown 427fca742c7SJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 4288a381b04SJed Brown 4298a381b04SJed Brown Level: advanced 4308a381b04SJed Brown 4311cc06b55SBarry Smith .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()` 4328a381b04SJed Brown @*/ 433d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXRegisterAll(void) 434d71ae5a4SJacob Faibussowitsch { 4358a381b04SJed Brown PetscFunctionBegin; 4363ba16761SJacob Faibussowitsch if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 4378a381b04SJed Brown TSARKIMEXRegisterAllCalled = PETSC_TRUE; 438e817cc15SEmil Constantinescu 439d27334e2SStefano Zampini #define RC PetscRealConstant 440d27334e2SStefano Zampini #define s2 RC(1.414213562373095048802) /* PetscSqrtReal((PetscReal)2.0) */ 441d27334e2SStefano Zampini #define us2 RC(0.2928932188134524755992) /* 1.0-1.0/PetscSqrtReal((PetscReal)2.0); */ 442d27334e2SStefano Zampini 443d27334e2SStefano Zampini /* Diagonally implicit methods */ 444e817cc15SEmil Constantinescu { 445d27334e2SStefano Zampini /* DIRK212, default of SUNDIALS */ 446d27334e2SStefano Zampini const PetscReal A[2][2] = { 447d27334e2SStefano Zampini {RC(1.0), RC(0.0)}, 448d27334e2SStefano Zampini {RC(-1.0), RC(1.0)} 449d27334e2SStefano Zampini }; 450d27334e2SStefano Zampini const PetscReal b[2] = {RC(0.5), RC(0.5)}; 451d27334e2SStefano Zampini const PetscReal bembed[2] = {RC(1.0), RC(0.0)}; 452d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKS212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 453d27334e2SStefano Zampini } 454d27334e2SStefano Zampini 4559371c9d4SSatish Balay { 456d27334e2SStefano Zampini /* ESDIRK12 from https://arxiv.org/pdf/1803.01613.pdf */ 457d27334e2SStefano Zampini const PetscReal A[2][2] = { 458d27334e2SStefano Zampini {RC(0.0), RC(0.0)}, 459d27334e2SStefano Zampini {RC(0.0), RC(1.0)} 460d27334e2SStefano Zampini }; 461d27334e2SStefano Zampini const PetscReal b[2] = {RC(0.0), RC(1.0)}; 462d27334e2SStefano Zampini const PetscReal bembed[2] = {RC(0.5), RC(0.5)}; 4633405e92cSStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES122SAL, 1, 2, &A[0][0], b, NULL, bembed, 1, b)); 464d27334e2SStefano Zampini } 465d27334e2SStefano Zampini 466d27334e2SStefano Zampini { 467d27334e2SStefano Zampini /* ESDIRK213L[2]SA from KC16. 4683405e92cSStefano Zampini TR-BDF2 from Hosea and Shampine 4693405e92cSStefano Zampini ESDIRK23 in https://arxiv.org/pdf/1803.01613.pdf */ 470d27334e2SStefano Zampini const PetscReal A[3][3] = { 471d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0)}, 472d27334e2SStefano Zampini {us2, us2, RC(0.0)}, 473d27334e2SStefano Zampini {s2 / RC(4.0), s2 / RC(4.0), us2 }, 474d27334e2SStefano Zampini }; 475d27334e2SStefano Zampini const PetscReal b[3] = {s2 / RC(4.0), s2 / RC(4.0), us2}; 476d27334e2SStefano Zampini 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)}; 477d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES213SAL, 2, 3, &A[0][0], b, NULL, bembed, 1, b)); 478d27334e2SStefano Zampini } 479d27334e2SStefano Zampini 480d27334e2SStefano Zampini { 481d27334e2SStefano Zampini #define g RC(0.43586652150845899941601945) 482d27334e2SStefano Zampini #define g2 PetscSqr(g) 483d27334e2SStefano Zampini #define g3 g *g2 484d27334e2SStefano Zampini #define g4 PetscSqr(g2) 485d27334e2SStefano Zampini #define g5 g *g4 486d27334e2SStefano Zampini #define c3 RC(1.0) 487d27334e2SStefano Zampini #define a32 c3 *(c3 - RC(2.0) * g) / (RC(4.0) * g) 488d27334e2SStefano Zampini #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)) 489d27334e2SStefano Zampini #define b3 (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) / (RC(3.0) * c3 * (c3 - RC(2.0) * g)) 490d27334e2SStefano Zampini #if 0 491d27334e2SStefano Zampini /* This is for c3 = 3/5 */ 492d27334e2SStefano Zampini #define bh2 \ 493d27334e2SStefano Zampini 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)) 494d27334e2SStefano Zampini #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)) 495d27334e2SStefano Zampini #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) 496d27334e2SStefano Zampini #else 497d27334e2SStefano Zampini /* This is for c3 = 1.0 */ 498d27334e2SStefano Zampini #define bh2 a32 499d27334e2SStefano Zampini #define bh3 g 500d27334e2SStefano Zampini #define bh4 RC(0.0) 501d27334e2SStefano Zampini #endif 502d27334e2SStefano Zampini /* ESDIRK3(2I)4L[2]SA from KC16 with c3 = 1.0 */ 503d27334e2SStefano Zampini /* Given by Kvaerno https://link.springer.com/article/10.1023/b:bitn.0000046811.70614.38 */ 504d27334e2SStefano Zampini const PetscReal A[4][4] = { 505d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0), RC(0.0)}, 506d27334e2SStefano Zampini {g, g, RC(0.0), RC(0.0)}, 507d27334e2SStefano Zampini {c3 - a32 - g, a32, g, RC(0.0)}, 508d27334e2SStefano Zampini {RC(1.0) - b2 - b3 - g, b2, b3, g }, 509d27334e2SStefano Zampini }; 510d27334e2SStefano Zampini const PetscReal b[4] = {RC(1.0) - b2 - b3 - g, b2, b3, g}; 511d27334e2SStefano Zampini const PetscReal bembed[4] = {RC(1.0) - bh2 - bh3 - bh4, bh2, bh3, bh4}; 512d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES324SAL, 3, 4, &A[0][0], b, NULL, bembed, 1, b)); 513d27334e2SStefano Zampini #undef g 514d27334e2SStefano Zampini #undef g2 515d27334e2SStefano Zampini #undef g3 516d27334e2SStefano Zampini #undef c3 517d27334e2SStefano Zampini #undef a32 518d27334e2SStefano Zampini #undef b2 519d27334e2SStefano Zampini #undef b3 520d27334e2SStefano Zampini #undef bh2 521d27334e2SStefano Zampini #undef bh3 522d27334e2SStefano Zampini #undef bh4 523d27334e2SStefano Zampini } 524d27334e2SStefano Zampini 525d27334e2SStefano Zampini { 526d27334e2SStefano Zampini /* ESDIRK3(2I)5L[2]SA from KC16 */ 527d27334e2SStefano Zampini const PetscReal A[5][5] = { 528d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 529d27334e2SStefano Zampini {RC(9.0) / RC(40.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0), RC(0.0) }, 530d27334e2SStefano Zampini {RC(19.0) / RC(72.0), RC(14.0) / RC(45.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0) }, 531d27334e2SStefano Zampini {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) }, 532d27334e2SStefano Zampini {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)}, 533d27334e2SStefano Zampini }; 534d27334e2SStefano Zampini 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)}; 535d27334e2SStefano Zampini 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)}; 536d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES325SAL, 3, 5, &A[0][0], b, NULL, bembed, 1, b)); 537d27334e2SStefano Zampini } 538d27334e2SStefano Zampini 539d27334e2SStefano Zampini { 540d27334e2SStefano Zampini // DIRK(6,6)[1]A[(7,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 541d27334e2SStefano Zampini const PetscReal A[7][7] = { 542d27334e2SStefano Zampini {RC(0.303487844706747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 543d27334e2SStefano Zampini {RC(-0.279756492709814), RC(0.500032236020747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 544d27334e2SStefano Zampini {RC(0.280583215743895), RC(-0.438560061586751), RC(0.217250734515736), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 545d27334e2SStefano Zampini {RC(-0.0677678738539846), RC(0.984312781232293), RC(-0.266720192540149), RC(0.2476680834526), RC(0.0), RC(0.0), RC(0.0) }, 546d27334e2SStefano Zampini {RC(0.125671616147993), RC(-0.995401751002415), RC(0.761333109549059), RC(-0.210281837202208), RC(0.866743712636936), RC(0.0), RC(0.0) }, 547d27334e2SStefano Zampini {RC(-0.368056238801488), RC(-0.999928082701516), RC(0.534734253232519), RC(-0.174856916279082), RC(0.615007160285509), RC(0.696549912132029), RC(0.0) }, 548d27334e2SStefano Zampini {RC(-0.00570546839653984), RC(-0.113110431835656), RC(-0.000965563207671587), RC(-0.000130490084629567), RC(0.00111737736895673), RC(-0.279385587378871), RC(0.618455906845342)} 549d27334e2SStefano Zampini }; 550d27334e2SStefano Zampini 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)}; 551d27334e2SStefano Zampini 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)}; 552d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK657A, 6, 7, &A[0][0], b, NULL, bembed, 1, b)); 553d27334e2SStefano Zampini } 554d27334e2SStefano Zampini { 555d27334e2SStefano Zampini // ESDIRK(8,6)[2]SA[(8,4)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 556d27334e2SStefano Zampini const PetscReal A[8][8] = { 557d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 558d27334e2SStefano Zampini {RC(0.333222149217725), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 559d27334e2SStefano Zampini {RC(0.0639743773182214), RC(-0.0830330224410214), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 560d27334e2SStefano Zampini {RC(-0.728522201369326), RC(-0.210414479522485), RC(0.532519916559342), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 561d27334e2SStefano Zampini {RC(-0.175135269272067), RC(0.666675582067552), RC(-0.304400907370867), RC(0.656797712445756), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0) }, 562d27334e2SStefano Zampini {RC(0.222695802705462), RC(-0.0948971794681061), RC(-0.0234336346686545), RC(-0.45385925012042), RC(0.0283910313826958), RC(0.333222149217725), RC(0.0), RC(0.0) }, 563d27334e2SStefano Zampini {RC(-0.132534078051299), RC(0.702597935004879), RC(-0.433316453128078), RC(0.893717488547587), RC(0.057381454791406), RC(-0.207798411552402), RC(0.333222149217725), RC(0.0) }, 564d27334e2SStefano Zampini {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)} 565d27334e2SStefano Zampini }; 566d27334e2SStefano Zampini 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)}; 567d27334e2SStefano Zampini 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)}; 568d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES648SA, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 569d27334e2SStefano Zampini } 570d27334e2SStefano Zampini { 571d27334e2SStefano Zampini // DIRK(8,6)[1]SAL[(8,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 572d27334e2SStefano Zampini const PetscReal A[8][8] = { 573d27334e2SStefano Zampini {RC(0.477264457385826), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 574d27334e2SStefano Zampini {RC(-0.197052588415002), RC(0.476363428459584), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 575d27334e2SStefano Zampini {RC(-0.0347674430372966), RC(0.633051807335483), RC(0.193634310075028), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 576d27334e2SStefano Zampini {RC(0.0967797668578702), RC(-0.193533526466535), RC(-0.000207622945800473), RC(0.159572204849431), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 577d27334e2SStefano Zampini {RC(0.162527231819875), RC(-0.249672513547382), RC(-0.0459079972041795), RC(0.36579476400859), RC(0.255752838307699), RC(0.0), RC(0.0), RC(0.0) }, 578d27334e2SStefano Zampini {RC(-0.00707603197171262), RC(0.846299854860295), RC(0.344020016925018), RC(-0.0720926054548865), RC(-0.215492331980875), RC(0.104341097622161), RC(0.0), RC(0.0) }, 579d27334e2SStefano Zampini {RC(0.00176857935179744), RC(0.0779960013127515), RC(0.303333277564557), RC(0.213160806732836), RC(0.351769320319038), RC(-0.381545894386538), RC(0.433517909105558), RC(0.0) }, 580d27334e2SStefano Zampini {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)} 581d27334e2SStefano Zampini }; 582d27334e2SStefano Zampini 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)}; 583d27334e2SStefano Zampini 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)}; 584d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK658A, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 585d27334e2SStefano Zampini } 586d27334e2SStefano Zampini { 587d27334e2SStefano Zampini // SDIRK(9,6)[1]SAL[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 588d27334e2SStefano Zampini const PetscReal A[9][9] = { 589d27334e2SStefano Zampini {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) }, 590d27334e2SStefano Zampini {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) }, 591d27334e2SStefano Zampini {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) }, 592d27334e2SStefano Zampini {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) }, 593d27334e2SStefano Zampini {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) }, 594d27334e2SStefano Zampini {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) }, 595d27334e2SStefano Zampini {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) }, 596d27334e2SStefano Zampini {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) }, 597d27334e2SStefano Zampini {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)} 598d27334e2SStefano Zampini }; 599d27334e2SStefano Zampini 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)}; 600d27334e2SStefano Zampini 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)}; 601d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKS659A, 6, 9, &A[0][0], b, NULL, bembed, 1, b)); 602d27334e2SStefano Zampini } 603d27334e2SStefano Zampini { 604d27334e2SStefano Zampini // DIRK(10,7)[1]SAL[(10,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 605d27334e2SStefano Zampini const PetscReal A[10][10] = { 606d27334e2SStefano Zampini {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) }, 607d27334e2SStefano Zampini {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) }, 608d27334e2SStefano Zampini {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) }, 609d27334e2SStefano Zampini {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) }, 610d27334e2SStefano Zampini {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) }, 611d27334e2SStefano Zampini {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) }, 612d27334e2SStefano Zampini {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) }, 613d27334e2SStefano Zampini {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) }, 614d27334e2SStefano Zampini {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) }, 615d27334e2SStefano Zampini {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)} 616d27334e2SStefano Zampini }; 617d27334e2SStefano Zampini 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)}; 618d27334e2SStefano Zampini const PetscReal bembed[10] = 619d27334e2SStefano Zampini {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)}; 620d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK7510SAL, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 621d27334e2SStefano Zampini } 622d27334e2SStefano Zampini { 623d27334e2SStefano Zampini // ESDIRK(10,7)[2]SA[(10,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 624d27334e2SStefano Zampini const PetscReal A[10][10] = { 625d27334e2SStefano Zampini {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) }, 626d27334e2SStefano Zampini {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) }, 627d27334e2SStefano Zampini {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) }, 628d27334e2SStefano Zampini {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) }, 629d27334e2SStefano Zampini {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) }, 630d27334e2SStefano Zampini {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) }, 631d27334e2SStefano Zampini {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) }, 632d27334e2SStefano Zampini {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) }, 633d27334e2SStefano Zampini {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) }, 634d27334e2SStefano Zampini {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)} 635d27334e2SStefano Zampini }; 636d27334e2SStefano Zampini const PetscReal b[10] = {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), 637d27334e2SStefano Zampini RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)}; 638d27334e2SStefano Zampini const PetscReal bembed[10] = {RC(-0.015494246543626), RC(0.167657963820093), RC(0.269858958144236), RC(-0.0443258997755156), RC(0.150049236875266), 639d27334e2SStefano Zampini RC(0.259452082755846), RC(0.244624573502521), RC(-0.215528446920284), RC(0.0487601760292619), RC(0.134945602112201)}; 640d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES7510SA, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 641d27334e2SStefano Zampini } 642d27334e2SStefano Zampini { 643d27334e2SStefano Zampini // DIRK(9,7)[1]A[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 644d27334e2SStefano Zampini const PetscReal A[9][9] = { 645d27334e2SStefano Zampini {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) }, 646d27334e2SStefano Zampini {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) }, 647d27334e2SStefano Zampini {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) }, 648d27334e2SStefano Zampini {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) }, 649d27334e2SStefano Zampini {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) }, 650d27334e2SStefano Zampini {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) }, 651d27334e2SStefano Zampini {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) }, 652d27334e2SStefano Zampini {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) }, 653d27334e2SStefano Zampini {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)} 654d27334e2SStefano Zampini }; 655d27334e2SStefano Zampini 656d27334e2SStefano Zampini 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)}; 657d27334e2SStefano Zampini 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)}; 658d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK759A, 7, 9, &A[0][0], b, NULL, bembed, 1, b)); 659d27334e2SStefano Zampini } 660d27334e2SStefano Zampini { 661d27334e2SStefano Zampini // SDIRK(11,7)[1]SAL[(11,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 662d27334e2SStefano Zampini const PetscReal A[11][11] = { 663d27334e2SStefano Zampini {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) }, 664d27334e2SStefano Zampini {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) }, 665d27334e2SStefano Zampini {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) }, 666d27334e2SStefano Zampini {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) }, 667d27334e2SStefano Zampini {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) }, 668d27334e2SStefano Zampini {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) }, 669d27334e2SStefano Zampini {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) }, 670d27334e2SStefano Zampini {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) }, 671d27334e2SStefano Zampini {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) }, 672d27334e2SStefano Zampini {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) }, 673d27334e2SStefano Zampini {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)} 674d27334e2SStefano Zampini }; 675d27334e2SStefano Zampini const PetscReal b[11] = 676d27334e2SStefano Zampini {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)}; 677d27334e2SStefano Zampini const PetscReal bembed[11] = 678d27334e2SStefano Zampini {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)}; 679d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKS7511SAL, 7, 11, &A[0][0], b, NULL, bembed, 1, b)); 680d27334e2SStefano Zampini } 681d27334e2SStefano Zampini { 682d27334e2SStefano Zampini // DIRK(13,8)[1]A[(14,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 683d27334e2SStefano Zampini const PetscReal A[14][14] = { 684d27334e2SStefano Zampini {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) }, 685d27334e2SStefano Zampini {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) }, 686d27334e2SStefano Zampini {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) }, 687d27334e2SStefano Zampini {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) }, 688d27334e2SStefano Zampini {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) }, 689d27334e2SStefano Zampini {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) }, 690d27334e2SStefano Zampini {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) }, 691d27334e2SStefano Zampini {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) }, 692d27334e2SStefano Zampini {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) }, 693d27334e2SStefano Zampini {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) }, 694d27334e2SStefano Zampini {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) }, 695d27334e2SStefano Zampini {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) }, 696d27334e2SStefano Zampini {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) }, 697d27334e2SStefano Zampini {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)} 698d27334e2SStefano Zampini }; 699d27334e2SStefano Zampini 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)}; 700d27334e2SStefano Zampini 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), 701d27334e2SStefano Zampini RC(0.0417664613347638), RC(0.223636394275293), RC(0.231569156867596), RC(0.240526201277663), RC(-0.222933582911926), RC(-0.0111479879597561), RC(0.19915314335888)}; 702d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK8614A, 8, 14, &A[0][0], b, NULL, bembed, 1, b)); 703d27334e2SStefano Zampini } 704d27334e2SStefano Zampini { 705d27334e2SStefano Zampini // DIRK(15,8)[1]SAL[(16,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 706d27334e2SStefano Zampini const PetscReal A[16][16] = { 707d27334e2SStefano Zampini {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) }, 708d27334e2SStefano Zampini {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) }, 709d27334e2SStefano Zampini {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) }, 710d27334e2SStefano Zampini {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) }, 711d27334e2SStefano Zampini {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) }, 712d27334e2SStefano Zampini {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) }, 713d27334e2SStefano Zampini {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) }, 714d27334e2SStefano Zampini {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) }, 715d27334e2SStefano Zampini {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) }, 716d27334e2SStefano Zampini {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) }, 717d27334e2SStefano Zampini {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) }, 718d27334e2SStefano Zampini {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) }, 719d27334e2SStefano Zampini {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) }, 720d27334e2SStefano Zampini {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) }, 721d27334e2SStefano Zampini {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) }, 722d27334e2SStefano Zampini {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)} 723d27334e2SStefano Zampini }; 724d27334e2SStefano Zampini 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)}; 725d27334e2SStefano Zampini 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), 726d27334e2SStefano Zampini 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)}; 727d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK8616SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 728d27334e2SStefano Zampini } 729d27334e2SStefano Zampini { 730d27334e2SStefano Zampini // ESDIRK(16,8)[2]SAL[(16,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 731d27334e2SStefano Zampini const PetscReal A[16][16] = { 732d27334e2SStefano Zampini {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) }, 733d27334e2SStefano Zampini {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) }, 734d27334e2SStefano Zampini {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) }, 735d27334e2SStefano Zampini {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) }, 736d27334e2SStefano Zampini {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) }, 737d27334e2SStefano Zampini {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) }, 738d27334e2SStefano Zampini {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) }, 739d27334e2SStefano Zampini {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) }, 740d27334e2SStefano Zampini {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) }, 741d27334e2SStefano Zampini {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) }, 742d27334e2SStefano Zampini {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) }, 743d27334e2SStefano Zampini {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) }, 744d27334e2SStefano Zampini {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) }, 745d27334e2SStefano Zampini {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) }, 746d27334e2SStefano Zampini {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) }, 747d27334e2SStefano Zampini {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)} 748d27334e2SStefano Zampini }; 749d27334e2SStefano Zampini 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), 750d27334e2SStefano Zampini RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)}; 751d27334e2SStefano Zampini 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), 752d27334e2SStefano Zampini RC(0.000900215457460583), RC(0.0547959465692338), RC(-0.334995726863153), RC(0.0464409662093384), RC(0.301388101652194), RC(0.00524851570622031), RC(0.229538601845236), RC(0.124643044573514)}; 753d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES8516SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 754d27334e2SStefano Zampini } 755d27334e2SStefano Zampini 756d27334e2SStefano Zampini /* Additive methods */ 757d27334e2SStefano Zampini { 758d27334e2SStefano Zampini const PetscReal A[3][3] = { 759e817cc15SEmil Constantinescu {0.0, 0.0, 0.0}, 7609371c9d4SSatish Balay {0.0, 0.0, 0.0}, 7619371c9d4SSatish Balay {0.0, 0.5, 0.0} 762d27334e2SStefano Zampini }; 763d27334e2SStefano Zampini const PetscReal At[3][3] = { 764d27334e2SStefano Zampini {1.0, 0.0, 0.0}, 765d27334e2SStefano Zampini {0.0, 0.5, 0.0}, 766d27334e2SStefano Zampini {0.0, 0.5, 0.5} 767d27334e2SStefano Zampini }; 768d27334e2SStefano Zampini const PetscReal b[3] = {0.0, 0.5, 0.5}; 769d27334e2SStefano Zampini const PetscReal bembedt[3] = {1.0, 0.0, 0.0}; 7709566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 771e817cc15SEmil Constantinescu } 7728a381b04SJed Brown { 773d27334e2SStefano Zampini const PetscReal A[2][2] = { 7749371c9d4SSatish Balay {0.0, 0.0}, 7759371c9d4SSatish Balay {0.5, 0.0} 776d27334e2SStefano Zampini }; 777d27334e2SStefano Zampini const PetscReal At[2][2] = { 778d27334e2SStefano Zampini {0.0, 0.0}, 779d27334e2SStefano Zampini {0.0, 0.5} 780d27334e2SStefano Zampini }; 781d27334e2SStefano Zampini const PetscReal b[2] = {0.0, 1.0}; 782d27334e2SStefano Zampini const PetscReal bembedt[2] = {0.5, 0.5}; 7831f80e275SEmil Constantinescu /* 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 */ 7849566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 7851f80e275SEmil Constantinescu } 7861f80e275SEmil Constantinescu { 787d27334e2SStefano Zampini const PetscReal A[2][2] = { 7889371c9d4SSatish Balay {0.0, 0.0}, 7899371c9d4SSatish Balay {1.0, 0.0} 790d27334e2SStefano Zampini }; 791d27334e2SStefano Zampini const PetscReal At[2][2] = { 792d27334e2SStefano Zampini {0.0, 0.0}, 793d27334e2SStefano Zampini {0.5, 0.5} 794d27334e2SStefano Zampini }; 795d27334e2SStefano Zampini const PetscReal b[2] = {0.5, 0.5}; 796d27334e2SStefano Zampini const PetscReal bembedt[2] = {0.0, 1.0}; 7971f80e275SEmil Constantinescu /* 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 */ 7989566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 7991f80e275SEmil Constantinescu } 8001f80e275SEmil Constantinescu { 801d27334e2SStefano Zampini const PetscReal A[2][2] = { 8029371c9d4SSatish Balay {0.0, 0.0}, 8039371c9d4SSatish Balay {1.0, 0.0} 804d27334e2SStefano Zampini }; 805d27334e2SStefano Zampini const PetscReal At[2][2] = { 806d27334e2SStefano Zampini {us2, 0.0}, 807d27334e2SStefano Zampini {1.0 - 2.0 * us2, us2} 808d27334e2SStefano Zampini }; 809d27334e2SStefano Zampini const PetscReal b[2] = {0.5, 0.5}; 810d27334e2SStefano Zampini const PetscReal bembedt[2] = {0.0, 1.0}; 811d27334e2SStefano Zampini const PetscReal binterpt[2][2] = { 812d27334e2SStefano Zampini {(us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))}, 813d27334e2SStefano Zampini {1 - (us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))} 814d27334e2SStefano Zampini }; 815d27334e2SStefano Zampini const PetscReal binterp[2][2] = { 816d27334e2SStefano Zampini {1.0, -0.5}, 817d27334e2SStefano Zampini {0.0, 0.5 } 818d27334e2SStefano Zampini }; 8199566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 8201f80e275SEmil Constantinescu } 8211f80e275SEmil Constantinescu { 822d27334e2SStefano Zampini const PetscReal A[3][3] = { 8239371c9d4SSatish Balay {0, 0, 0}, 824d27334e2SStefano Zampini {2 - s2, 0, 0}, 8259371c9d4SSatish Balay {0.5, 0.5, 0} 826d27334e2SStefano Zampini }; 827d27334e2SStefano Zampini const PetscReal At[3][3] = { 828d27334e2SStefano Zampini {0, 0, 0 }, 829d27334e2SStefano Zampini {1 - 1 / s2, 1 - 1 / s2, 0 }, 830d27334e2SStefano Zampini {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 831d27334e2SStefano Zampini }; 832d27334e2SStefano Zampini const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 833d27334e2SStefano Zampini const PetscReal binterpt[3][2] = { 834d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 835d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 836d27334e2SStefano Zampini {1.0 - s2, 1.0 / s2 } 837d27334e2SStefano Zampini }; 8389566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 8391f80e275SEmil Constantinescu } 8401f80e275SEmil Constantinescu { 841d27334e2SStefano Zampini const PetscReal A[3][3] = { 8429371c9d4SSatish Balay {0, 0, 0}, 843d27334e2SStefano Zampini {2 - s2, 0, 0}, 8449371c9d4SSatish Balay {0.75, 0.25, 0} 845d27334e2SStefano Zampini }; 846d27334e2SStefano Zampini const PetscReal At[3][3] = { 847d27334e2SStefano Zampini {0, 0, 0 }, 848d27334e2SStefano Zampini {1 - 1 / s2, 1 - 1 / s2, 0 }, 849d27334e2SStefano Zampini {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 850d27334e2SStefano Zampini }; 851d27334e2SStefano Zampini const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 852d27334e2SStefano Zampini const PetscReal binterpt[3][2] = { 853d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 854d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 855d27334e2SStefano Zampini {1.0 - s2, 1.0 / s2 } 856d27334e2SStefano Zampini }; 8579566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 8588a381b04SJed Brown } 85906db7b1cSJed Brown { /* Optimal for linear implicit part */ 860d27334e2SStefano Zampini const PetscReal A[3][3] = { 8619371c9d4SSatish Balay {0, 0, 0}, 862d27334e2SStefano Zampini {2 - s2, 0, 0}, 863d27334e2SStefano Zampini {(3 - 2 * s2) / 6, (3 + 2 * s2) / 6, 0} 864d27334e2SStefano Zampini }; 865d27334e2SStefano Zampini const PetscReal At[3][3] = { 866d27334e2SStefano Zampini {0, 0, 0 }, 867d27334e2SStefano Zampini {1 - 1 / s2, 1 - 1 / s2, 0 }, 868d27334e2SStefano Zampini {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 869d27334e2SStefano Zampini }; 870d27334e2SStefano Zampini const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 871d27334e2SStefano Zampini const PetscReal binterpt[3][2] = { 872d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 873d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 874d27334e2SStefano Zampini {1.0 - s2, 1.0 / s2 } 875d27334e2SStefano Zampini }; 8769566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 877a3a57f36SJed Brown } 8786cf0794eSJed Brown { /* Optimal for linear implicit part */ 879d27334e2SStefano Zampini const PetscReal A[3][3] = { 8809371c9d4SSatish Balay {0, 0, 0}, 8816cf0794eSJed Brown {0.5, 0, 0}, 8829371c9d4SSatish Balay {0.5, 0.5, 0} 883d27334e2SStefano Zampini }; 884d27334e2SStefano Zampini const PetscReal At[3][3] = { 885d27334e2SStefano Zampini {0.25, 0, 0 }, 886d27334e2SStefano Zampini {0, 0.25, 0 }, 887d27334e2SStefano Zampini {1. / 3, 1. / 3, 1. / 3} 888d27334e2SStefano Zampini }; 8899566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 8906cf0794eSJed Brown } 891a3a57f36SJed Brown { 892d27334e2SStefano Zampini const PetscReal A[4][4] = { 8939371c9d4SSatish Balay {0, 0, 0, 0}, 8944040e9f2SJed Brown {1767732205903. / 2027836641118., 0, 0, 0}, 8954040e9f2SJed Brown {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 8969371c9d4SSatish Balay {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 897d27334e2SStefano Zampini }; 898d27334e2SStefano Zampini const PetscReal At[4][4] = { 899d27334e2SStefano Zampini {0, 0, 0, 0 }, 900d27334e2SStefano Zampini {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0 }, 901d27334e2SStefano Zampini {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0 }, 902d27334e2SStefano Zampini {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.} 903d27334e2SStefano Zampini }; 904d27334e2SStefano Zampini const PetscReal bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}; 905d27334e2SStefano Zampini const PetscReal binterpt[4][2] = { 906d27334e2SStefano Zampini {4655552711362. / 22874653954995., -215264564351. / 13552729205753. }, 907d27334e2SStefano Zampini {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119. }, 908d27334e2SStefano Zampini {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, 909d27334e2SStefano Zampini {584795268549. / 6622622206610., 2508943948391. / 7218656332882. } 910d27334e2SStefano Zampini }; 9119566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 912a3a57f36SJed Brown } 913a3a57f36SJed Brown { 914d27334e2SStefano Zampini const PetscReal A[5][5] = { 9159371c9d4SSatish Balay {0, 0, 0, 0, 0}, 9166cf0794eSJed Brown {1. / 2, 0, 0, 0, 0}, 9176cf0794eSJed Brown {11. / 18, 1. / 18, 0, 0, 0}, 9186cf0794eSJed Brown {5. / 6, -5. / 6, .5, 0, 0}, 9199371c9d4SSatish Balay {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 920d27334e2SStefano Zampini }; 921d27334e2SStefano Zampini const PetscReal At[5][5] = { 922d27334e2SStefano Zampini {0, 0, 0, 0, 0 }, 923d27334e2SStefano Zampini {0, 1. / 2, 0, 0, 0 }, 924d27334e2SStefano Zampini {0, 1. / 6, 1. / 2, 0, 0 }, 925d27334e2SStefano Zampini {0, -1. / 2, 1. / 2, 1. / 2, 0 }, 926d27334e2SStefano Zampini {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2} 927d27334e2SStefano Zampini }; 928d27334e2SStefano Zampini PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 9296cf0794eSJed Brown } 9306cf0794eSJed Brown { 931d27334e2SStefano Zampini const PetscReal A[5][5] = { 9329371c9d4SSatish Balay {0, 0, 0, 0, 0}, 9336cf0794eSJed Brown {1, 0, 0, 0, 0}, 9346cf0794eSJed Brown {4. / 9, 2. / 9, 0, 0, 0}, 9356cf0794eSJed Brown {1. / 4, 0, 3. / 4, 0, 0}, 9369371c9d4SSatish Balay {1. / 4, 0, 3. / 5, 0, 0} 937d27334e2SStefano Zampini }; 938d27334e2SStefano Zampini const PetscReal At[5][5] = { 939d27334e2SStefano Zampini {0, 0, 0, 0, 0 }, 940d27334e2SStefano Zampini {.5, .5, 0, 0, 0 }, 941d27334e2SStefano Zampini {5. / 18, -1. / 9, .5, 0, 0 }, 942d27334e2SStefano Zampini {.5, 0, 0, .5, 0 }, 943d27334e2SStefano Zampini {.25, 0, .75, -.5, .5} 944d27334e2SStefano Zampini }; 945d27334e2SStefano Zampini PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 9466cf0794eSJed Brown } 9476cf0794eSJed Brown { 948d27334e2SStefano Zampini const PetscReal A[6][6] = { 9499371c9d4SSatish Balay {0, 0, 0, 0, 0, 0}, 950a3a57f36SJed Brown {1. / 2, 0, 0, 0, 0, 0}, 9514040e9f2SJed Brown {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 9524040e9f2SJed Brown {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 9534040e9f2SJed Brown {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 9549371c9d4SSatish Balay {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 955d27334e2SStefano Zampini }; 956d27334e2SStefano Zampini const PetscReal At[6][6] = { 957d27334e2SStefano Zampini {0, 0, 0, 0, 0, 0 }, 958d27334e2SStefano Zampini {1. / 4, 1. / 4, 0, 0, 0, 0 }, 959d27334e2SStefano Zampini {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0 }, 960d27334e2SStefano Zampini {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0 }, 961d27334e2SStefano Zampini {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0 }, 962d27334e2SStefano Zampini {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4} 963d27334e2SStefano Zampini }; 964d27334e2SStefano Zampini const PetscReal bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}; 965d27334e2SStefano Zampini const PetscReal binterpt[6][3] = { 966d27334e2SStefano Zampini {6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025. }, 967d27334e2SStefano Zampini {0, 0, 0 }, 968d27334e2SStefano Zampini {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035. }, 969d27334e2SStefano Zampini {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 970d27334e2SStefano Zampini {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469. }, 971d27334e2SStefano Zampini {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.} 972d27334e2SStefano Zampini }; 9739566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 974a3a57f36SJed Brown } 975a3a57f36SJed Brown { 976d27334e2SStefano Zampini const PetscReal A[8][8] = { 9779371c9d4SSatish Balay {0, 0, 0, 0, 0, 0, 0, 0}, 978a3a57f36SJed Brown {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 9794040e9f2SJed Brown {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 9804040e9f2SJed Brown {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 9814040e9f2SJed Brown {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 9824040e9f2SJed Brown {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 9834040e9f2SJed Brown {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 9849371c9d4SSatish Balay {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 985d27334e2SStefano Zampini }; 986d27334e2SStefano Zampini const PetscReal At[8][8] = { 987d27334e2SStefano Zampini {0, 0, 0, 0, 0, 0, 0, 0 }, 988d27334e2SStefano Zampini {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0 }, 989d27334e2SStefano Zampini {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0 }, 990d27334e2SStefano Zampini {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0 }, 991d27334e2SStefano Zampini {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0 }, 992d27334e2SStefano Zampini {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0 }, 993d27334e2SStefano Zampini {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0 }, 994d27334e2SStefano Zampini {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.} 995d27334e2SStefano Zampini }; 996d27334e2SStefano Zampini const PetscReal bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}; 997d27334e2SStefano Zampini const PetscReal binterpt[8][3] = { 998d27334e2SStefano Zampini {-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439. }, 999d27334e2SStefano Zampini {0, 0, 0 }, 1000d27334e2SStefano Zampini {0, 0, 0 }, 1001d27334e2SStefano Zampini {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, 1002d27334e2SStefano Zampini {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, 1003d27334e2SStefano Zampini {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, 1004d27334e2SStefano Zampini {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, 1005d27334e2SStefano Zampini {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460. } 1006d27334e2SStefano Zampini }; 10079566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 1008a3a57f36SJed Brown } 1009d27334e2SStefano Zampini #undef RC 1010d27334e2SStefano Zampini #undef us2 1011d27334e2SStefano Zampini #undef s2 10123ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10138a381b04SJed Brown } 10148a381b04SJed Brown 10158a381b04SJed Brown /*@C 1016bcf0153eSBarry Smith TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`. 10178a381b04SJed Brown 10188a381b04SJed Brown Not Collective 10198a381b04SJed Brown 10208a381b04SJed Brown Level: advanced 10218a381b04SJed Brown 10221cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 10238a381b04SJed Brown @*/ 1024d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXRegisterDestroy(void) 1025d71ae5a4SJacob Faibussowitsch { 10268a381b04SJed Brown ARKTableauLink link; 10278a381b04SJed Brown 10288a381b04SJed Brown PetscFunctionBegin; 10298a381b04SJed Brown while ((link = ARKTableauList)) { 10308a381b04SJed Brown ARKTableau t = &link->tab; 10318a381b04SJed Brown ARKTableauList = link->next; 10329566063dSJacob Faibussowitsch PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 10339566063dSJacob Faibussowitsch PetscCall(PetscFree2(t->bembedt, t->bembed)); 10349566063dSJacob Faibussowitsch PetscCall(PetscFree2(t->binterpt, t->binterp)); 10359566063dSJacob Faibussowitsch PetscCall(PetscFree(t->name)); 10369566063dSJacob Faibussowitsch PetscCall(PetscFree(link)); 10378a381b04SJed Brown } 10388a381b04SJed Brown TSARKIMEXRegisterAllCalled = PETSC_FALSE; 10393ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10408a381b04SJed Brown } 10418a381b04SJed Brown 10428a381b04SJed Brown /*@C 1043bcf0153eSBarry Smith TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called 1044bcf0153eSBarry Smith from `TSInitializePackage()`. 10458a381b04SJed Brown 10468a381b04SJed Brown Level: developer 10478a381b04SJed Brown 10481cc06b55SBarry Smith .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()` 10498a381b04SJed Brown @*/ 1050d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXInitializePackage(void) 1051d71ae5a4SJacob Faibussowitsch { 10528a381b04SJed Brown PetscFunctionBegin; 10533ba16761SJacob Faibussowitsch if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 10548a381b04SJed Brown TSARKIMEXPackageInitialized = PETSC_TRUE; 10559566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegisterAll()); 10569566063dSJacob Faibussowitsch PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 10573ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10588a381b04SJed Brown } 10598a381b04SJed Brown 10608a381b04SJed Brown /*@C 1061bcf0153eSBarry Smith TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is 1062bcf0153eSBarry Smith called from `PetscFinalize()`. 10638a381b04SJed Brown 10648a381b04SJed Brown Level: developer 10658a381b04SJed Brown 10661cc06b55SBarry Smith .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()` 10678a381b04SJed Brown @*/ 1068d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXFinalizePackage(void) 1069d71ae5a4SJacob Faibussowitsch { 10708a381b04SJed Brown PetscFunctionBegin; 10718a381b04SJed Brown TSARKIMEXPackageInitialized = PETSC_FALSE; 10729566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegisterDestroy()); 10733ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10748a381b04SJed Brown } 10758a381b04SJed Brown 1076cd652676SJed Brown /*@C 1077bcf0153eSBarry Smith TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 1078cd652676SJed Brown 1079cc4c1da9SBarry Smith Logically Collective 1080cd652676SJed Brown 1081cd652676SJed Brown Input Parameters: 1082cd652676SJed Brown + name - identifier for method 1083cd652676SJed Brown . order - approximation order of method 1084cd652676SJed Brown . s - number of stages, this is the dimension of the matrices below 1085cd652676SJed Brown . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 10860298fd71SBarry Smith . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 10870298fd71SBarry Smith . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 1088cd652676SJed Brown . A - Non-stiff stage coefficients (dimension s*s, row-major) 10890298fd71SBarry Smith . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 10900298fd71SBarry Smith . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 10910298fd71SBarry Smith . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 10920298fd71SBarry Smith . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 1093cd652676SJed Brown . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 1094cd652676SJed Brown . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 10950298fd71SBarry Smith - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 1096cd652676SJed Brown 1097cd652676SJed Brown Level: advanced 1098cd652676SJed Brown 1099bcf0153eSBarry Smith Note: 1100bcf0153eSBarry Smith Several `TSARKIMEX` methods are provided, this function is only needed to create new methods. 1101bcf0153eSBarry Smith 11021cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS` 1103cd652676SJed Brown @*/ 1104d71ae5a4SJacob Faibussowitsch 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[]) 1105d71ae5a4SJacob Faibussowitsch { 11068a381b04SJed Brown ARKTableauLink link; 11078a381b04SJed Brown ARKTableau t; 11088a381b04SJed Brown PetscInt i, j; 11098a381b04SJed Brown 11108a381b04SJed Brown PetscFunctionBegin; 1111a748edf9SJed Brown PetscCheck(s > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected number of stages s %" PetscInt_FMT " > 0", s); 11129566063dSJacob Faibussowitsch PetscCall(TSARKIMEXInitializePackage()); 1113d27334e2SStefano Zampini for (link = ARKTableauList; link; link = link->next) { 1114d27334e2SStefano Zampini PetscBool match; 1115d27334e2SStefano Zampini 1116d27334e2SStefano Zampini PetscCall(PetscStrcmp(link->tab.name, name, &match)); 1117d27334e2SStefano Zampini PetscCheck(!match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Method with name \"%s\" already registered", name); 1118d27334e2SStefano Zampini } 11199566063dSJacob Faibussowitsch PetscCall(PetscNew(&link)); 11208a381b04SJed Brown t = &link->tab; 11219566063dSJacob Faibussowitsch PetscCall(PetscStrallocpy(name, &t->name)); 11228a381b04SJed Brown t->order = order; 11238a381b04SJed Brown t->s = s; 11249566063dSJacob Faibussowitsch PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 11259566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->At, At, s * s)); 1126d27334e2SStefano Zampini if (A) { 11279566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->A, A, s * s)); 1128d27334e2SStefano Zampini t->additive = PETSC_TRUE; 1129d27334e2SStefano Zampini } 1130d27334e2SStefano Zampini 11319566063dSJacob Faibussowitsch if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 11329371c9d4SSatish Balay else 11339371c9d4SSatish Balay for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 1134d27334e2SStefano Zampini 1135d27334e2SStefano Zampini if (t->additive) { 11369566063dSJacob Faibussowitsch if (b) PetscCall(PetscArraycpy(t->b, b, s)); 11379371c9d4SSatish Balay else 11389371c9d4SSatish Balay for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 1139d27334e2SStefano Zampini } else PetscCall(PetscArrayzero(t->b, s)); 1140d27334e2SStefano Zampini 11419566063dSJacob Faibussowitsch if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 11429371c9d4SSatish Balay else 11439371c9d4SSatish Balay for (i = 0; i < s; i++) 11449371c9d4SSatish Balay for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 1145d27334e2SStefano Zampini 1146d27334e2SStefano Zampini if (t->additive) { 11479566063dSJacob Faibussowitsch if (c) PetscCall(PetscArraycpy(t->c, c, s)); 11489371c9d4SSatish Balay else 11499371c9d4SSatish Balay for (i = 0; i < s; i++) 11509371c9d4SSatish Balay for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 1151d27334e2SStefano Zampini } else PetscCall(PetscArrayzero(t->c, s)); 1152d27334e2SStefano Zampini 1153e817cc15SEmil Constantinescu t->stiffly_accurate = PETSC_TRUE; 11549371c9d4SSatish Balay for (i = 0; i < s; i++) 11559371c9d4SSatish Balay if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 1156d27334e2SStefano Zampini 1157e817cc15SEmil Constantinescu t->explicit_first_stage = PETSC_TRUE; 11589371c9d4SSatish Balay for (i = 0; i < s; i++) 11599371c9d4SSatish Balay if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 1160d27334e2SStefano Zampini 1161e817cc15SEmil Constantinescu /* def of FSAL can be made more precise */ 11624e9d4bf5SJed Brown t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 1163d27334e2SStefano Zampini 1164108c343cSJed Brown if (bembedt) { 11659566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 11669566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 11679566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 1168108c343cSJed Brown } 1169108c343cSJed Brown 11704f385281SJed Brown t->pinterp = pinterp; 11719566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 11729566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 11739566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 1174d27334e2SStefano Zampini 11758a381b04SJed Brown link->next = ARKTableauList; 11768a381b04SJed Brown ARKTableauList = link; 11773ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 11788a381b04SJed Brown } 11798a381b04SJed Brown 1180d27334e2SStefano Zampini /*@C 1181d27334e2SStefano Zampini TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation 1182d27334e2SStefano Zampini 1183d27334e2SStefano Zampini Logically Collective. 1184d27334e2SStefano Zampini 1185d27334e2SStefano Zampini Input Parameters: 1186d27334e2SStefano Zampini + name - identifier for method 1187d27334e2SStefano Zampini . order - approximation order of method 1188d27334e2SStefano Zampini . s - number of stages, this is the dimension of the matrices below 1189d27334e2SStefano Zampini . At - Butcher table of stage coefficients (dimension `s`*`s`, row-major order) 1190d27334e2SStefano Zampini . bt - Butcher table for completing the step (dimension `s`; pass `NULL` to use the last row of `At`) 1191d27334e2SStefano Zampini . ct - Abscissa of each stage (dimension s, NULL to use row sums of At) 1192d27334e2SStefano Zampini . bembedt - Stiff part of completion table for embedded method (dimension s; `NULL` if not available) 1193d27334e2SStefano Zampini . pinterp - Order of the interpolation scheme, equal to the number of columns of `binterpt` and `binterp` 1194d27334e2SStefano Zampini - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 1195d27334e2SStefano Zampini 1196d27334e2SStefano Zampini Level: advanced 1197d27334e2SStefano Zampini 1198d27334e2SStefano Zampini Note: 1199d27334e2SStefano Zampini Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods. 1200d27334e2SStefano Zampini 1201d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS` 1202d27334e2SStefano Zampini @*/ 1203d27334e2SStefano Zampini 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[]) 1204d27334e2SStefano Zampini { 1205d27334e2SStefano Zampini PetscFunctionBegin; 1206d27334e2SStefano Zampini PetscCall(TSARKIMEXRegister(name, order, s, At, bt, ct, NULL, NULL, NULL, bembedt, NULL, pinterp, binterpt, NULL)); 1207d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 1208d27334e2SStefano Zampini } 1209d27334e2SStefano Zampini 1210108c343cSJed Brown /* 1211108c343cSJed Brown The step completion formula is 1212108c343cSJed Brown 1213108c343cSJed Brown x1 = x0 - h bt^T YdotI + h b^T YdotRHS 1214108c343cSJed Brown 1215108c343cSJed Brown This function can be called before or after ts->vec_sol has been updated. 1216108c343cSJed Brown Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 1217108c343cSJed Brown We can write 1218108c343cSJed Brown 1219108c343cSJed Brown x1e = x0 - h bet^T YdotI + h be^T YdotRHS 1220108c343cSJed Brown = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 1221108c343cSJed Brown = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 1222108c343cSJed Brown 1223108c343cSJed Brown so we can evaluate the method with different order even after the step has been optimistically completed. 1224108c343cSJed Brown */ 1225d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) 1226d71ae5a4SJacob Faibussowitsch { 1227108c343cSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1228108c343cSJed Brown ARKTableau tab = ark->tableau; 1229108c343cSJed Brown PetscScalar *w = ark->work; 1230108c343cSJed Brown PetscReal h; 1231108c343cSJed Brown PetscInt s = tab->s, j; 1232d27334e2SStefano Zampini PetscBool hasE; 1233108c343cSJed Brown 1234108c343cSJed Brown PetscFunctionBegin; 1235108c343cSJed Brown switch (ark->status) { 1236108c343cSJed Brown case TS_STEP_INCOMPLETE: 1237d71ae5a4SJacob Faibussowitsch case TS_STEP_PENDING: 1238d71ae5a4SJacob Faibussowitsch h = ts->time_step; 1239d71ae5a4SJacob Faibussowitsch break; 1240d71ae5a4SJacob Faibussowitsch case TS_STEP_COMPLETE: 1241d71ae5a4SJacob Faibussowitsch h = ts->ptime - ts->ptime_prev; 1242d71ae5a4SJacob Faibussowitsch break; 1243d71ae5a4SJacob Faibussowitsch default: 1244d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1245108c343cSJed Brown } 1246108c343cSJed Brown if (order == tab->order) { 1247e817cc15SEmil Constantinescu if (ark->status == TS_STEP_INCOMPLETE) { 1248740132f1SEmil Constantinescu if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 12499566063dSJacob Faibussowitsch PetscCall(VecCopy(ark->Y[s - 1], X)); 1250e817cc15SEmil Constantinescu } else { /* Use the standard completion formula (bt,b) */ 12519566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, X)); 1252e817cc15SEmil Constantinescu for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 12539566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1254d27334e2SStefano Zampini if (tab->additive && ark->imex) { /* Method is IMEX, complete the explicit formula */ 1255d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1256d27334e2SStefano Zampini if (hasE) { 1257108c343cSJed Brown for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 12589566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1259e817cc15SEmil Constantinescu } 1260e817cc15SEmil Constantinescu } 1261d27334e2SStefano Zampini } 12629566063dSJacob Faibussowitsch } else PetscCall(VecCopy(ts->vec_sol, X)); 1263108c343cSJed Brown if (done) *done = PETSC_TRUE; 12643ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1265108c343cSJed Brown } else if (order == tab->order - 1) { 1266108c343cSJed Brown if (!tab->bembedt) goto unavailable; 1267108c343cSJed Brown if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 12689566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, X)); 1269e817cc15SEmil Constantinescu for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 12709566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1271d27334e2SStefano Zampini if (tab->additive) { 1272d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1273d27334e2SStefano Zampini if (hasE) { 1274108c343cSJed Brown for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 12759566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1276d27334e2SStefano Zampini } 1277d27334e2SStefano Zampini } 1278108c343cSJed Brown } else { /* Rollback and re-complete using (bet-be,be-b) */ 12799566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, X)); 1280e817cc15SEmil Constantinescu for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 12819566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 1282d27334e2SStefano Zampini if (tab->additive) { 1283d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1284d27334e2SStefano Zampini if (hasE) { 1285108c343cSJed Brown for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 12869566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1287108c343cSJed Brown } 1288d27334e2SStefano Zampini } 1289d27334e2SStefano Zampini } 1290108c343cSJed Brown if (done) *done = PETSC_TRUE; 12913ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1292108c343cSJed Brown } 1293108c343cSJed Brown unavailable: 1294108c343cSJed Brown if (done) *done = PETSC_FALSE; 12959371c9d4SSatish Balay else 12969371c9d4SSatish Balay 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, 12979371c9d4SSatish Balay tab->order, order); 12983ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1299108c343cSJed Brown } 1300108c343cSJed Brown 1301d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 1302d71ae5a4SJacob Faibussowitsch { 1303c79dcfbdSBarry Smith Vec Udot, Y1, Y2; 1304c79dcfbdSBarry Smith TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1305c79dcfbdSBarry Smith PetscReal norm; 1306c79dcfbdSBarry Smith 1307c79dcfbdSBarry Smith PetscFunctionBegin; 13089566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 13099566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 13109566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 13119566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 13129566063dSJacob Faibussowitsch PetscCall(VecSetRandom(Udot, NULL)); 13139566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 13149566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y2, -1.0, Y1)); 13159566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y2, -1.0, Udot)); 13169566063dSJacob Faibussowitsch PetscCall(VecNorm(Y2, NORM_2, &norm)); 1317c79dcfbdSBarry Smith if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 1318c79dcfbdSBarry Smith *id = PETSC_TRUE; 1319c79dcfbdSBarry Smith } else { 1320c79dcfbdSBarry Smith *id = PETSC_FALSE; 1321835f2295SStefano Zampini 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)); 1322c79dcfbdSBarry Smith } 13239566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Udot)); 13249566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Y1)); 13259566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Y2)); 13263ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1327c79dcfbdSBarry Smith } 1328c79dcfbdSBarry Smith 13290467964bSStefano Zampini static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS, PetscReal, Vec, IS *); 13300467964bSStefano Zampini 1331d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSStep_ARKIMEX(TS ts) 1332d71ae5a4SJacob Faibussowitsch { 13338a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 13348a381b04SJed Brown ARKTableau tab = ark->tableau; 13358a381b04SJed Brown const PetscInt s = tab->s; 133624655328SShri const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 1337406d0ec2SJed Brown PetscScalar *w = ark->work; 13381297b224SEmil Constantinescu Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 133996400bd6SLisandro Dalcin PetscBool extrapolate = ark->extrapolate; 1340108c343cSJed Brown TSAdapt adapt; 13418a381b04SJed Brown SNES snes; 1342fecfb714SLisandro Dalcin PetscInt i, j, its, lits; 1343be5899b3SLisandro Dalcin PetscInt rejections = 0; 1344d27334e2SStefano Zampini PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE; 134596400bd6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 13468a381b04SJed Brown 13478a381b04SJed Brown PetscFunctionBegin; 134896400bd6SLisandro Dalcin if (ark->extrapolate && !ark->Y_prev) { 13499566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 13509566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1351d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 135296400bd6SLisandro Dalcin } 135396400bd6SLisandro Dalcin 1354d27334e2SStefano Zampini if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE)); 1355d27334e2SStefano Zampini if (!hasE) dirk = PETSC_TRUE; 1356d27334e2SStefano Zampini 135796400bd6SLisandro Dalcin if (!ts->steprollback) { 1358d27334e2SStefano Zampini if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 13599566063dSJacob Faibussowitsch PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 136096400bd6SLisandro Dalcin } 1361fecfb714SLisandro Dalcin if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 136296400bd6SLisandro Dalcin for (i = 0; i < s; i++) { 13639566063dSJacob Faibussowitsch PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 13649566063dSJacob Faibussowitsch PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 1365d27334e2SStefano Zampini if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 136696400bd6SLisandro Dalcin } 136796400bd6SLisandro Dalcin } 136896400bd6SLisandro Dalcin } 136996400bd6SLisandro Dalcin 13703b98415fSStefano Zampini /* 13713b98415fSStefano Zampini For fully implicit formulations we must solve the equations 13723b98415fSStefano Zampini 13733b98415fSStefano Zampini F(t_n,x_n,xdot) = 0 13743b98415fSStefano Zampini 13753b98415fSStefano Zampini for the explicit first stage. 13763b98415fSStefano Zampini Here we call SNESSolve using PETSC_MAX_REAL as shift to flag it. 13773b98415fSStefano Zampini Special handling is inside SNESTSFormFunction_ARKIMEX and SNESTSFormJacobian_ARKIMEX 1378c6bf8827SStefano Zampini We compute Ydot0 if we restart the step or if we resized the problem after remeshing 13793b98415fSStefano Zampini */ 1380c6bf8827SStefano Zampini if (dirk && tab->explicit_first_stage && (ts->steprestart || ts->stepresize)) { 138198940a98SHong Zhang ark->scoeff = PETSC_MAX_REAL; 1382d27334e2SStefano Zampini PetscCall(VecCopy(ts->vec_sol, Z)); 13830467964bSStefano Zampini if (!ark->alg_is) { 13840467964bSStefano Zampini PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ts->ptime, Z, &ark->alg_is)); 13850467964bSStefano Zampini PetscCall(ISViewFromOptions(ark->alg_is, (PetscObject)ts, "-ts_arkimex_algebraic_is_view")); 13860467964bSStefano Zampini } 1387d27334e2SStefano Zampini PetscCall(TSGetSNES(ts, &snes)); 13880467964bSStefano Zampini PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, 1)); 1389d27334e2SStefano Zampini PetscCall(SNESSolve(snes, NULL, Ydot0)); 1390c6bf8827SStefano Zampini if (ark->alg_is) PetscCall(VecISSet(Ydot0, ark->alg_is, 0.0)); 13910467964bSStefano Zampini PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, -1)); 1392d27334e2SStefano Zampini } 1393d27334e2SStefano Zampini 1394d27334e2SStefano Zampini /* For IMEX we compute a step */ 1395d27334e2SStefano Zampini if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 139696400bd6SLisandro Dalcin TS ts_start; 1397d27334e2SStefano Zampini if (PetscDefined(USE_DEBUG) && hasE) { 1398c79dcfbdSBarry Smith PetscBool id = PETSC_FALSE; 13999566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 14008434afd1SBarry Smith 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"); 1401c79dcfbdSBarry Smith } 14029566063dSJacob Faibussowitsch PetscCall(TSClone(ts, &ts_start)); 14039566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 14049566063dSJacob Faibussowitsch PetscCall(TSSetTime(ts_start, ts->ptime)); 14059566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 14069566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 14079566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 14089566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 14099566063dSJacob Faibussowitsch PetscCall(TSSetType(ts_start, TSARKIMEX)); 14109566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 14119566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 141234561852SEmil Constantinescu 14139566063dSJacob Faibussowitsch PetscCall(TSRestartStep(ts_start)); 14149566063dSJacob Faibussowitsch PetscCall(TSSolve(ts_start, ts->vec_sol)); 14159566063dSJacob Faibussowitsch PetscCall(TSGetTime(ts_start, &ts->ptime)); 14169566063dSJacob Faibussowitsch PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 1417bbd56ea5SKarl Rupp 141885fc7851SLisandro Dalcin { /* Save the initial slope for the next step */ 141985fc7851SLisandro Dalcin TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 14209566063dSJacob Faibussowitsch PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 142185fc7851SLisandro Dalcin } 142296400bd6SLisandro Dalcin ts->steps++; 142334561852SEmil Constantinescu 1424d15a3a53SEmil Constantinescu /* Set the correct TS in SNES */ 1425d15a3a53SEmil Constantinescu /* We'll try to bypass this by changing the method on the fly */ 142696400bd6SLisandro Dalcin { 14279566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 14289566063dSJacob Faibussowitsch PetscCall(TSSetSNES(ts, snes)); 142996400bd6SLisandro Dalcin } 14309566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts_start)); 1431e817cc15SEmil Constantinescu } 1432e817cc15SEmil Constantinescu 1433108c343cSJed Brown ark->status = TS_STEP_INCOMPLETE; 143496400bd6SLisandro Dalcin while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 143596400bd6SLisandro Dalcin PetscReal t = ts->ptime; 1436108c343cSJed Brown PetscReal h = ts->time_step; 14378a381b04SJed Brown for (i = 0; i < s; i++) { 14389be3e283SDebojyoti Ghosh ark->stage_time = t + h * ct[i]; 14399566063dSJacob Faibussowitsch PetscCall(TSPreStage(ts, ark->stage_time)); 14408a381b04SJed Brown if (At[i * s + i] == 0) { /* This stage is explicit */ 14413c633725SBarry Smith 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"); 14429566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, Y[i])); 1443e817cc15SEmil Constantinescu for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 14449566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 1445d27334e2SStefano Zampini if (tab->additive && hasE) { 14468a381b04SJed Brown for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 14479566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 1448d27334e2SStefano Zampini } 144912b1dd1aSStefano Zampini PetscCall(TSGetSNES(ts, &snes)); 145012b1dd1aSStefano Zampini PetscCall(SNESResetCounters(snes)); 14518a381b04SJed Brown } else { 1452b296d7d5SJed Brown ark->scoeff = 1. / At[i * s + i]; 14538a381b04SJed Brown /* Ydot = shift*(Y-Z) */ 14549566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, Z)); 1455e817cc15SEmil Constantinescu for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 14569566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Z, i, w, YdotI)); 1457d27334e2SStefano Zampini if (tab->additive && hasE) { 1458c58d1302SEmil Constantinescu for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 14599566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 1460d27334e2SStefano Zampini } 1461fecfb714SLisandro Dalcin if (extrapolate && !ts->steprestart) { 146256dcabbaSDebojyoti Ghosh /* Initial guess extrapolated from previous time step stage values */ 14639566063dSJacob Faibussowitsch PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 146456dcabbaSDebojyoti Ghosh } else { 14658a381b04SJed Brown /* Initial guess taken from last stage */ 14669566063dSJacob Faibussowitsch PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 146756dcabbaSDebojyoti Ghosh } 14689566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 14699566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, Y[i])); 14709566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its)); 14719566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 14729371c9d4SSatish Balay ts->snes_its += its; 14739371c9d4SSatish Balay ts->ksp_its += lits; 14749566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 14759566063dSJacob Faibussowitsch PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 147696400bd6SLisandro Dalcin if (!stageok) { 14771be93e3eSJed Brown /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 14781be93e3eSJed Brown * use extrapolation to initialize the solves on the next attempt. */ 147996400bd6SLisandro Dalcin extrapolate = PETSC_FALSE; 14801be93e3eSJed Brown goto reject_step; 14811be93e3eSJed Brown } 14828a381b04SJed Brown } 1483d27334e2SStefano Zampini if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { 1484e817cc15SEmil Constantinescu if (i == 0 && tab->explicit_first_stage) { 1485d27334e2SStefano Zampini 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", 1486d27334e2SStefano Zampini ((PetscObject)ts)->type_name, ark->tableau->name); 14879566063dSJacob Faibussowitsch PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 1488e817cc15SEmil Constantinescu } else { 14899566063dSJacob Faibussowitsch PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1490e817cc15SEmil Constantinescu } 1491e817cc15SEmil Constantinescu } else { 14925eca1a21SEmil Constantinescu if (i == 0 && tab->explicit_first_stage) { 14939566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Ydot)); 14949566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 14959566063dSJacob Faibussowitsch PetscCall(VecScale(YdotI[i], -1.0)); 14965eca1a21SEmil Constantinescu } else { 14979566063dSJacob Faibussowitsch PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 14985eca1a21SEmil Constantinescu } 1499d27334e2SStefano Zampini if (hasE) { 15004cc180ffSJed Brown if (ark->imex) { 15019566063dSJacob Faibussowitsch PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 15024cc180ffSJed Brown } else { 15039566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(YdotRHS[i])); 15044cc180ffSJed Brown } 15058a381b04SJed Brown } 1506d27334e2SStefano Zampini } 15079566063dSJacob Faibussowitsch PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 1508e817cc15SEmil Constantinescu } 150996400bd6SLisandro Dalcin 1510be5899b3SLisandro Dalcin ark->status = TS_STEP_INCOMPLETE; 15119566063dSJacob Faibussowitsch PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 1512108c343cSJed Brown ark->status = TS_STEP_PENDING; 15139566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 15149566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidatesClear(adapt)); 15159566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 15169566063dSJacob Faibussowitsch PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 151796400bd6SLisandro Dalcin ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 151896400bd6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1519c61711c8SStefano Zampini PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1520be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1521be5899b3SLisandro Dalcin goto reject_step; 152296400bd6SLisandro Dalcin } 152396400bd6SLisandro Dalcin 15248a381b04SJed Brown ts->ptime += ts->time_step; 1525cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 1526108c343cSJed Brown break; 152796400bd6SLisandro Dalcin 152896400bd6SLisandro Dalcin reject_step: 15299371c9d4SSatish Balay ts->reject++; 15309371c9d4SSatish Balay accept = PETSC_FALSE; 1531be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 153296400bd6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 153363a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1534108c343cSJed Brown } 1535f85781f1SEmil Constantinescu } 15363ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 15378a381b04SJed Brown } 1538d27334e2SStefano Zampini 153980ab5e31SHong Zhang /* 154080ab5e31SHong Zhang This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 154180ab5e31SHong Zhang Udot = H(t,U) + G(t,U) 154280ab5e31SHong Zhang This corresponds to F(t,U,Udot) = Udot-H(t,U). 154380ab5e31SHong Zhang 154480ab5e31SHong Zhang The complete adjoint equations are 154580ab5e31SHong Zhang (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 15465d7a6ebeSHong Zhang dGdU (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 15475d7a6ebeSHong Zhang + dHdU (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 154880ab5e31SHong Zhang lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 154980ab5e31SHong Zhang mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 15505d7a6ebeSHong Zhang + dGdP (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 15515d7a6ebeSHong Zhang + dHdP (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 155280ab5e31SHong Zhang mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 155380ab5e31SHong Zhang where shift = 1/(at[i][i]*h) 155480ab5e31SHong Zhang 155580ab5e31SHong Zhang If at[i][i] is 0, the first equation falls back to 155680ab5e31SHong Zhang lambda_s[i] = h ( 155780ab5e31SHong Zhang (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 155880ab5e31SHong Zhang + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 155980ab5e31SHong Zhang 156080ab5e31SHong Zhang */ 156180ab5e31SHong Zhang static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 156280ab5e31SHong Zhang { 156380ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 156480ab5e31SHong Zhang ARKTableau tab = ark->tableau; 156580ab5e31SHong Zhang const PetscInt s = tab->s; 156680ab5e31SHong Zhang const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 156780ab5e31SHong Zhang PetscScalar *w = ark->work; 156880ab5e31SHong Zhang Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 156980ab5e31SHong Zhang Mat Jex, Jim, Jimpre; 157080ab5e31SHong Zhang PetscInt i, j, nadj; 157180ab5e31SHong Zhang PetscReal t = ts->ptime, stage_time_ex; 157280ab5e31SHong Zhang PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 157380ab5e31SHong Zhang KSP ksp; 157480ab5e31SHong Zhang 157580ab5e31SHong Zhang PetscFunctionBegin; 157680ab5e31SHong Zhang ark->status = TS_STEP_INCOMPLETE; 157780ab5e31SHong Zhang PetscCall(SNESGetKSP(ts->snes, &ksp)); 157880ab5e31SHong Zhang PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 157980ab5e31SHong Zhang PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 158080ab5e31SHong Zhang 158180ab5e31SHong Zhang for (i = s - 1; i >= 0; i--) { 158280ab5e31SHong Zhang ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 158380ab5e31SHong Zhang stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 158480ab5e31SHong Zhang if (At[i * s + i] == 0) { // This stage is explicit 158580ab5e31SHong Zhang ark->scoeff = 0.; 158680ab5e31SHong Zhang } else { 158780ab5e31SHong Zhang ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 158880ab5e31SHong Zhang } 158980ab5e31SHong Zhang PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 159080ab5e31SHong Zhang PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 159180ab5e31SHong Zhang if (ts->vecs_sensip) { 159280ab5e31SHong Zhang 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 159380ab5e31SHong Zhang PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 159480ab5e31SHong Zhang } 159580ab5e31SHong Zhang /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 15965d7a6ebeSHong Zhang for (nadj = 0; nadj < ts->numcost; nadj++) { 15975d7a6ebeSHong Zhang /* build implicit part */ 15985d7a6ebeSHong Zhang PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 159980ab5e31SHong Zhang if (s - i - 1 > 0) { 160080ab5e31SHong Zhang /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 160180ab5e31SHong Zhang for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 160280ab5e31SHong Zhang PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 16035d7a6ebeSHong Zhang } 16045d7a6ebeSHong Zhang /* Temp = Temp - bt[i] lambda_{n+1} */ 16055d7a6ebeSHong Zhang PetscCall(VecAXPY(VecsSensiTemp[nadj], -bt[i], ts->vecs_sensi[nadj])); 16065d7a6ebeSHong Zhang if (bt[i] || s - i - 1 > 0) { 160780ab5e31SHong Zhang /* (shift I - dHdU) Temp */ 160880ab5e31SHong Zhang PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 160980ab5e31SHong Zhang /* cancel out shift Temp where shift=-scoeff/h */ 161080ab5e31SHong Zhang PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 161180ab5e31SHong Zhang if (ts->vecs_sensip) { 161280ab5e31SHong Zhang /* - dHdP Temp */ 161380ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 16145d7a6ebeSHong Zhang /* mu_n += -h dHdP Temp */ 16155d7a6ebeSHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 161680ab5e31SHong Zhang } 16175d7a6ebeSHong Zhang } else { 16185d7a6ebeSHong Zhang PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 16195d7a6ebeSHong Zhang } 16205d7a6ebeSHong Zhang /* build explicit part */ 16215d7a6ebeSHong Zhang PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 16225d7a6ebeSHong Zhang if (s - i - 1 > 0) { 162380ab5e31SHong Zhang /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 162480ab5e31SHong Zhang for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 162580ab5e31SHong Zhang PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 16265d7a6ebeSHong Zhang } 16275d7a6ebeSHong Zhang /* Temp = Temp + b[i] lambda_{n+1} */ 16285d7a6ebeSHong Zhang PetscCall(VecAXPY(VecsSensiTemp[nadj], b[i], ts->vecs_sensi[nadj])); 16295d7a6ebeSHong Zhang if (b[i] || s - i - 1 > 0) { 163080ab5e31SHong Zhang /* dGdU Temp */ 163180ab5e31SHong Zhang PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 163280ab5e31SHong Zhang if (ts->vecs_sensip) { 163380ab5e31SHong Zhang /* dGdP Temp */ 163480ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 163580ab5e31SHong Zhang /* mu_n += h dGdP Temp */ 163680ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 163780ab5e31SHong Zhang } 163880ab5e31SHong Zhang } 163980ab5e31SHong Zhang /* Build LHS for first-order adjoint */ 164080ab5e31SHong Zhang if (At[i * s + i] == 0) { // This stage is explicit 164180ab5e31SHong Zhang PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 164280ab5e31SHong Zhang } else { 164380ab5e31SHong Zhang KSP ksp; 164480ab5e31SHong Zhang KSPConvergedReason kspreason; 164580ab5e31SHong Zhang PetscCall(SNESGetKSP(ts->snes, &ksp)); 164680ab5e31SHong Zhang PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 164780ab5e31SHong Zhang PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 164880ab5e31SHong Zhang PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 164980ab5e31SHong Zhang PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 165080ab5e31SHong Zhang if (kspreason < 0) { 165180ab5e31SHong Zhang ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 165280ab5e31SHong Zhang PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj)); 165380ab5e31SHong Zhang } 165480ab5e31SHong Zhang if (ts->vecs_sensip) { 165580ab5e31SHong Zhang /* -dHdP lambda_s[i] */ 165680ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 165780ab5e31SHong Zhang /* mu_n += h at[i][i] dHdP lambda_s[i] */ 165880ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 165980ab5e31SHong Zhang } 166080ab5e31SHong Zhang } 166180ab5e31SHong Zhang } 166280ab5e31SHong Zhang } 166380ab5e31SHong Zhang for (j = 0; j < s; j++) w[j] = 1.0; 166480ab5e31SHong Zhang for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 166580ab5e31SHong Zhang PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 166680ab5e31SHong Zhang ark->status = TS_STEP_COMPLETE; 166780ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 166880ab5e31SHong Zhang } 16698a381b04SJed Brown 1670d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1671d71ae5a4SJacob Faibussowitsch { 1672cd652676SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1673d27334e2SStefano Zampini ARKTableau tab = ark->tableau; 1674d27334e2SStefano Zampini PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1675108c343cSJed Brown PetscReal h; 1676108c343cSJed Brown PetscReal tt, t; 1677d27334e2SStefano Zampini PetscScalar *bt = ark->work, *b = ark->work + s; 1678d27334e2SStefano Zampini const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1679cd652676SJed Brown 1680cd652676SJed Brown PetscFunctionBegin; 1681d27334e2SStefano Zampini PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name); 1682108c343cSJed Brown switch (ark->status) { 1683108c343cSJed Brown case TS_STEP_INCOMPLETE: 1684108c343cSJed Brown case TS_STEP_PENDING: 1685108c343cSJed Brown h = ts->time_step; 1686108c343cSJed Brown t = (itime - ts->ptime) / h; 1687108c343cSJed Brown break; 1688108c343cSJed Brown case TS_STEP_COMPLETE: 1689be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1690108c343cSJed Brown t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1691108c343cSJed Brown break; 1692d71ae5a4SJacob Faibussowitsch default: 1693d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1694108c343cSJed Brown } 1695cd652676SJed Brown for (i = 0; i < s; i++) bt[i] = b[i] = 0; 16964f385281SJed Brown for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1697cd652676SJed Brown for (i = 0; i < s; i++) { 1698c1758d98SDebojyoti Ghosh bt[i] += h * Bt[i * pinterp + j] * tt; 1699108c343cSJed Brown b[i] += h * B[i * pinterp + j] * tt; 1700cd652676SJed Brown } 1701cd652676SJed Brown } 17029566063dSJacob Faibussowitsch PetscCall(VecCopy(ark->Y[0], X)); 17039566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1704d27334e2SStefano Zampini if (tab->additive) { 1705d27334e2SStefano Zampini PetscBool hasE; 1706d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1707d27334e2SStefano Zampini if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1708d27334e2SStefano Zampini } 17093ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1710cd652676SJed Brown } 1711cd652676SJed Brown 1712d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1713d71ae5a4SJacob Faibussowitsch { 171456dcabbaSDebojyoti Ghosh TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1715d27334e2SStefano Zampini ARKTableau tab = ark->tableau; 1716d27334e2SStefano Zampini PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1717be5899b3SLisandro Dalcin PetscReal h, h_prev, t, tt; 1718d27334e2SStefano Zampini PetscScalar *bt = ark->work, *b = ark->work + s; 1719d27334e2SStefano Zampini const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 172056dcabbaSDebojyoti Ghosh 172156dcabbaSDebojyoti Ghosh PetscFunctionBegin; 17223c633725SBarry Smith PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 172381d12688SDebojyoti Ghosh h = ts->time_step; 1724be5899b3SLisandro Dalcin h_prev = ts->ptime - ts->ptime_prev; 1725be5899b3SLisandro Dalcin t = 1 + h / h_prev * c; 1726d27334e2SStefano Zampini for (i = 0; i < s; i++) bt[i] = b[i] = 0; 172756dcabbaSDebojyoti Ghosh for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 172856dcabbaSDebojyoti Ghosh for (i = 0; i < s; i++) { 172981d12688SDebojyoti Ghosh bt[i] += h * Bt[i * pinterp + j] * tt; 173056dcabbaSDebojyoti Ghosh b[i] += h * B[i * pinterp + j] * tt; 173156dcabbaSDebojyoti Ghosh } 173256dcabbaSDebojyoti Ghosh } 17333c633725SBarry Smith PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 17349566063dSJacob Faibussowitsch PetscCall(VecCopy(ark->Y_prev[0], X)); 17359566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1736d27334e2SStefano Zampini if (tab->additive) { 1737d27334e2SStefano Zampini PetscBool hasE; 1738d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1739d27334e2SStefano Zampini if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1740d27334e2SStefano Zampini } 17413ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 174256dcabbaSDebojyoti Ghosh } 174356dcabbaSDebojyoti Ghosh 1744d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1745d71ae5a4SJacob Faibussowitsch { 174696400bd6SLisandro Dalcin TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 174796400bd6SLisandro Dalcin ARKTableau tab = ark->tableau; 174896400bd6SLisandro Dalcin 174996400bd6SLisandro Dalcin PetscFunctionBegin; 17503ba16761SJacob Faibussowitsch if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 17519566063dSJacob Faibussowitsch PetscCall(PetscFree(ark->work)); 17529566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 17539566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 17549566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 17559566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 17569566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 17579566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 17583ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 175996400bd6SLisandro Dalcin } 176096400bd6SLisandro Dalcin 1761d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSReset_ARKIMEX(TS ts) 1762d71ae5a4SJacob Faibussowitsch { 17638a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 17648a381b04SJed Brown 17658a381b04SJed Brown PetscFunctionBegin; 17663a2a065fSHong Zhang if (ark->fastslowsplit) { 17673a2a065fSHong Zhang PetscTryMethod(ts, "TSReset_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 17683a2a065fSHong Zhang } else { 17699566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTableauReset(ts)); 17709566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ark->Ydot)); 17719566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ark->Ydot0)); 17729566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ark->Z)); 17730467964bSStefano Zampini PetscCall(ISDestroy(&ark->alg_is)); 17743a2a065fSHong Zhang } 17753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 17768a381b04SJed Brown } 17778a381b04SJed Brown 177880ab5e31SHong Zhang static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 177980ab5e31SHong Zhang { 178080ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 178180ab5e31SHong Zhang ARKTableau tab = ark->tableau; 178280ab5e31SHong Zhang 178380ab5e31SHong Zhang PetscFunctionBegin; 178480ab5e31SHong Zhang PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 178580ab5e31SHong Zhang PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 178680ab5e31SHong Zhang PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 178780ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 178880ab5e31SHong Zhang } 178980ab5e31SHong Zhang 1790d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1791d71ae5a4SJacob Faibussowitsch { 1792d5e6173cSPeter Brune TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1793d5e6173cSPeter Brune 1794d5e6173cSPeter Brune PetscFunctionBegin; 1795d5e6173cSPeter Brune if (Z) { 17964ce06fd1SHong Zhang if (dm && dm != ts->dm) { 17979566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1798d5e6173cSPeter Brune } else *Z = ax->Z; 1799d5e6173cSPeter Brune } 1800d5e6173cSPeter Brune if (Ydot) { 18014ce06fd1SHong Zhang if (dm && dm != ts->dm) { 18029566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1803d5e6173cSPeter Brune } else *Ydot = ax->Ydot; 1804d5e6173cSPeter Brune } 18053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1806d5e6173cSPeter Brune } 1807d5e6173cSPeter Brune 1808d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1809d71ae5a4SJacob Faibussowitsch { 1810d5e6173cSPeter Brune PetscFunctionBegin; 1811d5e6173cSPeter Brune if (Z) { 18124ce06fd1SHong Zhang if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1813d5e6173cSPeter Brune } 1814d5e6173cSPeter Brune if (Ydot) { 18154ce06fd1SHong Zhang if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1816d5e6173cSPeter Brune } 18173ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1818d5e6173cSPeter Brune } 1819d5e6173cSPeter Brune 18200467964bSStefano Zampini /* 18210467964bSStefano Zampini DAEs need special handling for algebraic variables when restarting DIRK methods with explicit 18220467964bSStefano Zampini first stage. In particular, we need: 18230467964bSStefano Zampini - to zero the nonlinear function (in case the dual variables are not consistent in the first step) 18240467964bSStefano Zampini - to modify the preconditioning matrix by calling MatZeroRows with identity on these variables. 18250467964bSStefano Zampini */ 18260467964bSStefano Zampini static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is) 18270467964bSStefano Zampini { 18280467964bSStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 18290467964bSStefano Zampini DM dm; 18300467964bSStefano Zampini Vec F, W, Xdot; 18310467964bSStefano Zampini const PetscScalar *w; 18320467964bSStefano Zampini PetscInt nz = 0, n, st; 18330467964bSStefano Zampini PetscInt *nzr; 18340467964bSStefano Zampini 18350467964bSStefano Zampini PetscFunctionBegin; 18360467964bSStefano Zampini PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */ 18370467964bSStefano Zampini PetscCall(DMGetGlobalVector(dm, &Xdot)); 18380467964bSStefano Zampini PetscCall(DMGetGlobalVector(dm, &F)); 18390467964bSStefano Zampini PetscCall(DMGetGlobalVector(dm, &W)); 18400467964bSStefano Zampini PetscCall(VecSet(Xdot, 0.0)); 18410467964bSStefano Zampini PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex)); 18420467964bSStefano Zampini PetscCall(VecSetRandom(Xdot, NULL)); 18430467964bSStefano Zampini PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex)); 18440467964bSStefano Zampini PetscCall(VecAXPY(W, -1.0, F)); 18450467964bSStefano Zampini PetscCall(VecGetOwnershipRange(W, &st, NULL)); 18460467964bSStefano Zampini PetscCall(VecGetLocalSize(W, &n)); 18470467964bSStefano Zampini PetscCall(VecGetArrayRead(W, &w)); 18480467964bSStefano Zampini for (PetscInt i = 0; i < n; i++) 18490467964bSStefano Zampini if (w[i] == 0.0) nz++; 18500467964bSStefano Zampini PetscCall(PetscMalloc1(nz, &nzr)); 18510467964bSStefano Zampini nz = 0; 18520467964bSStefano Zampini for (PetscInt i = 0; i < n; i++) 18530467964bSStefano Zampini if (w[i] == 0.0) nzr[nz++] = i + st; 18540467964bSStefano Zampini PetscCall(VecRestoreArrayRead(W, &w)); 18550467964bSStefano Zampini PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is)); 18560467964bSStefano Zampini PetscCall(DMRestoreGlobalVector(dm, &Xdot)); 18570467964bSStefano Zampini PetscCall(DMRestoreGlobalVector(dm, &F)); 18580467964bSStefano Zampini PetscCall(DMRestoreGlobalVector(dm, &W)); 18590467964bSStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 18600467964bSStefano Zampini } 18610467964bSStefano Zampini 18620467964bSStefano Zampini /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure 18630467964bSStefano Zampini at the finest level, in the DM for coarser solves. */ 18640467964bSStefano Zampini static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is) 18650467964bSStefano Zampini { 18660467964bSStefano Zampini TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 18670467964bSStefano Zampini 18680467964bSStefano Zampini PetscFunctionBegin; 18690467964bSStefano Zampini if (dm && dm != ts->dm) { 18700467964bSStefano Zampini PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is)); 18710467964bSStefano Zampini } else *alg_is = ax->alg_is; 18720467964bSStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 18730467964bSStefano Zampini } 18743b98415fSStefano Zampini 18753b98415fSStefano Zampini /* This defines the nonlinear equation that is to be solved with SNES */ 1876d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1877d71ae5a4SJacob Faibussowitsch { 18788a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1879d5e6173cSPeter Brune DM dm, dmsave; 18804ce06fd1SHong Zhang Vec Z, Ydot; 18810467964bSStefano Zampini IS alg_is; 18828a381b04SJed Brown 18838a381b04SJed Brown PetscFunctionBegin; 18849566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 18859566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 18860467964bSStefano Zampini if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 18870467964bSStefano Zampini 1888d5e6173cSPeter Brune dmsave = ts->dm; 1889d5e6173cSPeter Brune ts->dm = dm; 1890740132f1SEmil Constantinescu 189198940a98SHong Zhang if (ark->scoeff == PETSC_MAX_REAL) { 18923b98415fSStefano Zampini /* We are solving F(t_n,x_n,xdot) = 0 to start the method */ 18930467964bSStefano Zampini if (!alg_is) { 18940467964bSStefano Zampini PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 18950467964bSStefano Zampini PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is)); 18960467964bSStefano Zampini PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is)); 18970467964bSStefano Zampini PetscCall(PetscObjectDereference((PetscObject)alg_is)); 18980467964bSStefano Zampini PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view")); 18990467964bSStefano Zampini } 19004ce06fd1SHong Zhang PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 19010467964bSStefano Zampini PetscCall(VecISSet(F, alg_is, 0.0)); 1902d27334e2SStefano Zampini } else { 190398940a98SHong Zhang PetscReal shift = ark->scoeff / ts->time_step; 1904d27334e2SStefano Zampini PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 19054ce06fd1SHong Zhang PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1906d27334e2SStefano Zampini } 1907e817cc15SEmil Constantinescu 1908d5e6173cSPeter Brune ts->dm = dmsave; 19099566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 19103ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19118a381b04SJed Brown } 19128a381b04SJed Brown 1913d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1914d71ae5a4SJacob Faibussowitsch { 19158a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1916d5e6173cSPeter Brune DM dm, dmsave; 19174ce06fd1SHong Zhang Vec Ydot, Z; 191898940a98SHong Zhang PetscReal shift; 19190467964bSStefano Zampini IS alg_is; 19208a381b04SJed Brown 19218a381b04SJed Brown PetscFunctionBegin; 19229566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 19238a381b04SJed Brown /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 19240467964bSStefano Zampini PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 19250467964bSStefano Zampini /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */ 19260467964bSStefano Zampini if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 19270467964bSStefano Zampini 1928d5e6173cSPeter Brune dmsave = ts->dm; 1929d5e6173cSPeter Brune ts->dm = dm; 1930740132f1SEmil Constantinescu 193198940a98SHong Zhang if (ark->scoeff == PETSC_MAX_REAL) { 19323b98415fSStefano Zampini PetscBool hasZeroRows; 19333b98415fSStefano Zampini 19343b98415fSStefano Zampini /* We are solving F(t_n,x_n,xdot) = 0 to start the method 19350467964bSStefano Zampini We compute with a very large shift and then scale back the matrix */ 1936d27334e2SStefano Zampini shift = 1.0 / PETSC_MACHINE_EPSILON; 19374ce06fd1SHong Zhang PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1938d27334e2SStefano Zampini PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 19393b98415fSStefano Zampini PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows)); 19403b98415fSStefano Zampini if (hasZeroRows) { 19410467964bSStefano Zampini PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 19423b98415fSStefano Zampini /* the default of AIJ is to not keep the pattern! We should probably change it someday */ 19433b98415fSStefano Zampini PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 19443b98415fSStefano Zampini PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL)); 19453b98415fSStefano Zampini } 19463b98415fSStefano Zampini PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view")); 1947d27334e2SStefano Zampini if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1948d27334e2SStefano Zampini } else { 194998940a98SHong Zhang shift = ark->scoeff / ts->time_step; 19504ce06fd1SHong Zhang PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1951d27334e2SStefano Zampini } 1952d5e6173cSPeter Brune ts->dm = dmsave; 1953d27334e2SStefano Zampini PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 19543ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1955d5e6173cSPeter Brune } 1956d5e6173cSPeter Brune 195780ab5e31SHong Zhang static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 195880ab5e31SHong Zhang { 195980ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 196080ab5e31SHong Zhang 196180ab5e31SHong Zhang PetscFunctionBegin; 196280ab5e31SHong Zhang if (ns) *ns = ark->tableau->s; 196380ab5e31SHong Zhang if (Y) *Y = ark->Y; 196480ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 196580ab5e31SHong Zhang } 196680ab5e31SHong Zhang 1967d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1968d71ae5a4SJacob Faibussowitsch { 1969d5e6173cSPeter Brune PetscFunctionBegin; 19703ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1971d5e6173cSPeter Brune } 1972d5e6173cSPeter Brune 1973d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1974d71ae5a4SJacob Faibussowitsch { 1975d5e6173cSPeter Brune TS ts = (TS)ctx; 1976d5e6173cSPeter Brune Vec Z, Z_c; 1977d5e6173cSPeter Brune 1978d5e6173cSPeter Brune PetscFunctionBegin; 19799566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 19809566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 19819566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Z, Z_c)); 19829566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 19839566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 19849566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 19853ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19868a381b04SJed Brown } 19878a381b04SJed Brown 1988d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1989d71ae5a4SJacob Faibussowitsch { 1990cdb298fcSPeter Brune PetscFunctionBegin; 19913ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1992cdb298fcSPeter Brune } 1993cdb298fcSPeter Brune 1994d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1995d71ae5a4SJacob Faibussowitsch { 1996cdb298fcSPeter Brune TS ts = (TS)ctx; 1997cdb298fcSPeter Brune Vec Z, Z_c; 1998cdb298fcSPeter Brune 1999cdb298fcSPeter Brune PetscFunctionBegin; 20009566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 20019566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 2002cdb298fcSPeter Brune 20039566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 20049566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2005cdb298fcSPeter Brune 20069566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 20079566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 20083ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2009cdb298fcSPeter Brune } 2010cdb298fcSPeter Brune 2011d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 2012d71ae5a4SJacob Faibussowitsch { 201396400bd6SLisandro Dalcin TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 201496400bd6SLisandro Dalcin ARKTableau tab = ark->tableau; 201596400bd6SLisandro Dalcin 201696400bd6SLisandro Dalcin PetscFunctionBegin; 2017d27334e2SStefano Zampini PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 20189566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 20199566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 2020d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 202196400bd6SLisandro Dalcin if (ark->extrapolate) { 20229566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 20239566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 2024d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 202596400bd6SLisandro Dalcin } 20263ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 202796400bd6SLisandro Dalcin } 202896400bd6SLisandro Dalcin 2029d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 2030d71ae5a4SJacob Faibussowitsch { 20318a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2032d5e6173cSPeter Brune DM dm; 203396400bd6SLisandro Dalcin SNES snes; 2034f9c1d6abSBarry Smith 20358a381b04SJed Brown PetscFunctionBegin; 20363a2a065fSHong Zhang if (ark->fastslowsplit) { 20373a2a065fSHong Zhang PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 20383a2a065fSHong Zhang } else { 20399566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTableauSetUp(ts)); 20409566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 20419566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 20429566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 20439566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm)); 20449566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 20459566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 20469566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 20473a2a065fSHong Zhang PetscCall(SNESSetDM(snes, dm)); 20483a2a065fSHong Zhang } 20493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20508a381b04SJed Brown } 20518a381b04SJed Brown 205280ab5e31SHong Zhang static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 205380ab5e31SHong Zhang { 205480ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 205580ab5e31SHong Zhang ARKTableau tab = ark->tableau; 205680ab5e31SHong Zhang 205780ab5e31SHong Zhang PetscFunctionBegin; 205880ab5e31SHong Zhang PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 205980ab5e31SHong Zhang PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 206080ab5e31SHong Zhang if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); } 206180ab5e31SHong Zhang if (PetscDefined(USE_DEBUG)) { 206280ab5e31SHong Zhang PetscBool id = PETSC_FALSE; 206380ab5e31SHong Zhang PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 20648434afd1SBarry Smith 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"); 206580ab5e31SHong Zhang } 206680ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 206780ab5e31SHong Zhang } 206880ab5e31SHong Zhang 2069ce78bad3SBarry Smith static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems PetscOptionsObject) 2070d71ae5a4SJacob Faibussowitsch { 20714cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2072d27334e2SStefano Zampini PetscBool dirk; 20738a381b04SJed Brown 20748a381b04SJed Brown PetscFunctionBegin; 2075d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2076d27334e2SStefano Zampini PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 20778a381b04SJed Brown { 20788a381b04SJed Brown ARKTableauLink link; 20798a381b04SJed Brown PetscInt count, choice; 20808a381b04SJed Brown PetscBool flg; 20818a381b04SJed Brown const char **namelist; 2082d27334e2SStefano Zampini for (link = ARKTableauList, count = 0; link; link = link->next) { 2083d27334e2SStefano Zampini if (!dirk && link->tab.additive) count++; 2084d27334e2SStefano Zampini if (dirk && !link->tab.additive) count++; 2085d27334e2SStefano Zampini } 20869566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(count, (char ***)&namelist)); 2087d27334e2SStefano Zampini for (link = ARKTableauList, count = 0; link; link = link->next) { 2088d27334e2SStefano Zampini if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 2089d27334e2SStefano Zampini if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 2090d27334e2SStefano Zampini } 2091d27334e2SStefano Zampini if (dirk) { 2092d27334e2SStefano Zampini PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2093d27334e2SStefano Zampini if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 2094d27334e2SStefano Zampini } else { 20953a2a065fSHong Zhang PetscBool fastslowsplit; 20969566063dSJacob Faibussowitsch PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 20979566063dSJacob Faibussowitsch if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 20984cc180ffSJed Brown flg = (PetscBool)!ark->imex; 20999566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 21004cc180ffSJed Brown ark->imex = (PetscBool)!flg; 21013a2a065fSHong Zhang PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg)); 21023a2a065fSHong Zhang if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit)); 2103cb0310aaSHong Zhang PetscCall(TSARKIMEXGetFastSlowSplit(ts, &fastslowsplit)); 2104cb0310aaSHong Zhang if (fastslowsplit) { 2105cb0310aaSHong Zhang SNES snes; 2106cb0310aaSHong Zhang 2107cb0310aaSHong Zhang PetscCall(TSRHSSplitGetSNES(ts, &snes)); 2108cb0310aaSHong Zhang PetscCall(SNESSetFromOptions(snes)); 2109cb0310aaSHong Zhang } 2110d27334e2SStefano Zampini } 2111d27334e2SStefano Zampini PetscCall(PetscFree(namelist)); 21129566063dSJacob Faibussowitsch 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)); 21138a381b04SJed Brown } 2114d0609cedSBarry Smith PetscOptionsHeadEnd(); 21153ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 21168a381b04SJed Brown } 21178a381b04SJed Brown 2118d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 2119d71ae5a4SJacob Faibussowitsch { 21208a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2121d27334e2SStefano Zampini PetscBool iascii, dirk; 21228a381b04SJed Brown 21238a381b04SJed Brown PetscFunctionBegin; 2124d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 21259566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 21268a381b04SJed Brown if (iascii) { 2127d27334e2SStefano Zampini PetscViewerFormat format; 21289c334d8fSLisandro Dalcin ARKTableau tab = ark->tableau; 212919fd82e9SBarry Smith TSARKIMEXType arktype; 2130d27334e2SStefano Zampini char buf[2048]; 21313a28c0e5SStefano Zampini PetscBool flg; 21323a28c0e5SStefano Zampini 21339566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetType(ts, &arktype)); 21349566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 2135d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 21369566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 2137d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 2138d27334e2SStefano Zampini PetscCall(PetscViewerGetFormat(viewer, &format)); 2139d27334e2SStefano Zampini if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2140d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 2141d27334e2SStefano Zampini for (PetscInt i = 0; i < tab->s; i++) { 2142d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 2143d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 2144d27334e2SStefano Zampini } 2145d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 2146d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 2147d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 2148d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 2149d27334e2SStefano Zampini } 21509566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 21519566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 21529566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 21539566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 2154d27334e2SStefano Zampini if (!dirk) { 2155d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 21569566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 21578a381b04SJed Brown } 2158d27334e2SStefano Zampini } 21593ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 21608a381b04SJed Brown } 21618a381b04SJed Brown 2162d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 2163d71ae5a4SJacob Faibussowitsch { 2164f2c2a1b9SBarry Smith SNES snes; 21659c334d8fSLisandro Dalcin TSAdapt adapt; 2166f2c2a1b9SBarry Smith 2167f2c2a1b9SBarry Smith PetscFunctionBegin; 21689566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 21699566063dSJacob Faibussowitsch PetscCall(TSAdaptLoad(adapt, viewer)); 21709566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 21719566063dSJacob Faibussowitsch PetscCall(SNESLoad(snes, viewer)); 2172ad6bc421SBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 21739566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 21749566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 21753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2176f2c2a1b9SBarry Smith } 2177f2c2a1b9SBarry Smith 2178cc4c1da9SBarry Smith /*@ 2179bcf0153eSBarry Smith TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 21808a381b04SJed Brown 218120f4b53cSBarry Smith Logically Collective 21828a381b04SJed Brown 2183d8d19677SJose E. Roman Input Parameters: 21848a381b04SJed Brown + ts - timestepping context 2185bcf0153eSBarry Smith - arktype - type of `TSARKIMEX` scheme 21868a381b04SJed Brown 2187bcf0153eSBarry Smith Options Database Key: 2188bcf0153eSBarry Smith . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 21899bd3e852SBarry Smith 21908a381b04SJed Brown Level: intermediate 21918a381b04SJed Brown 21921cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 2193db781477SPatrick Sanan `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 21948a381b04SJed Brown @*/ 2195d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 2196d71ae5a4SJacob Faibussowitsch { 21978a381b04SJed Brown PetscFunctionBegin; 21988a381b04SJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 21994f572ea9SToby Isaac PetscAssertPointer(arktype, 2); 2200cac4c232SBarry Smith PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 22013ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22028a381b04SJed Brown } 22038a381b04SJed Brown 2204cc4c1da9SBarry Smith /*@ 2205bcf0153eSBarry Smith TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 22068a381b04SJed Brown 220720f4b53cSBarry Smith Logically Collective 22088a381b04SJed Brown 22098a381b04SJed Brown Input Parameter: 22108a381b04SJed Brown . ts - timestepping context 22118a381b04SJed Brown 22128a381b04SJed Brown Output Parameter: 2213bcf0153eSBarry Smith . arktype - type of `TSARKIMEX` scheme 22148a381b04SJed Brown 22158a381b04SJed Brown Level: intermediate 22168a381b04SJed Brown 2217*bfe80ac4SPierre Jolivet .seealso: [](ch_ts), `TSARKIMEX` 22188a381b04SJed Brown @*/ 2219d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2220d71ae5a4SJacob Faibussowitsch { 22218a381b04SJed Brown PetscFunctionBegin; 22228a381b04SJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2223cac4c232SBarry Smith PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 22243ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22258a381b04SJed Brown } 22268a381b04SJed Brown 222716353aafSBarry Smith /*@ 2228bcf0153eSBarry Smith TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 22294cc180ffSJed Brown 223020f4b53cSBarry Smith Logically Collective 22314cc180ffSJed Brown 2232d8d19677SJose E. Roman Input Parameters: 22334cc180ffSJed Brown + ts - timestepping context 2234bcf0153eSBarry Smith - flg - `PETSC_TRUE` for fully implicit 22354cc180ffSJed Brown 22364cc180ffSJed Brown Level: intermediate 22374cc180ffSJed Brown 22381cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 22394cc180ffSJed Brown @*/ 2240d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 2241d71ae5a4SJacob Faibussowitsch { 22424cc180ffSJed Brown PetscFunctionBegin; 22434cc180ffSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 22443a28c0e5SStefano Zampini PetscValidLogicalCollectiveBool(ts, flg, 2); 2245cac4c232SBarry Smith PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 22463ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22474cc180ffSJed Brown } 22484cc180ffSJed Brown 22493a28c0e5SStefano Zampini /*@ 22503a28c0e5SStefano Zampini TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 22513a28c0e5SStefano Zampini 225220f4b53cSBarry Smith Logically Collective 22533a28c0e5SStefano Zampini 22543a28c0e5SStefano Zampini Input Parameter: 22553a28c0e5SStefano Zampini . ts - timestepping context 22563a28c0e5SStefano Zampini 22577a7aea1fSJed Brown Output Parameter: 2258bcf0153eSBarry Smith . flg - `PETSC_TRUE` for fully implicit 22593a28c0e5SStefano Zampini 22603a28c0e5SStefano Zampini Level: intermediate 22613a28c0e5SStefano Zampini 22621cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 22633a28c0e5SStefano Zampini @*/ 2264d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 2265d71ae5a4SJacob Faibussowitsch { 22663a28c0e5SStefano Zampini PetscFunctionBegin; 22673a28c0e5SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 22684f572ea9SToby Isaac PetscAssertPointer(flg, 2); 2269cac4c232SBarry Smith PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 22703ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22713a28c0e5SStefano Zampini } 22723a28c0e5SStefano Zampini 2273d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 2274d71ae5a4SJacob Faibussowitsch { 22758a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 22768a381b04SJed Brown 22778a381b04SJed Brown PetscFunctionBegin; 22788a381b04SJed Brown *arktype = ark->tableau->name; 22793ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22808a381b04SJed Brown } 2281d27334e2SStefano Zampini 2282d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 2283d71ae5a4SJacob Faibussowitsch { 22848a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 22858a381b04SJed Brown PetscBool match; 22868a381b04SJed Brown ARKTableauLink link; 22878a381b04SJed Brown 22888a381b04SJed Brown PetscFunctionBegin; 22898a381b04SJed Brown if (ark->tableau) { 22909566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 22913ba16761SJacob Faibussowitsch if (match) PetscFunctionReturn(PETSC_SUCCESS); 22928a381b04SJed Brown } 22938a381b04SJed Brown for (link = ARKTableauList; link; link = link->next) { 22949566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 22958a381b04SJed Brown if (match) { 22969566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 22978a381b04SJed Brown ark->tableau = &link->tab; 22989566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 22992ffb9264SLisandro Dalcin ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 23003ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 23018a381b04SJed Brown } 23028a381b04SJed Brown } 230398921bdaSJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 23048a381b04SJed Brown } 2305e0877f53SBarry Smith 2306d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 2307d71ae5a4SJacob Faibussowitsch { 23084cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 23094cc180ffSJed Brown 23104cc180ffSJed Brown PetscFunctionBegin; 23114cc180ffSJed Brown ark->imex = (PetscBool)!flg; 23123ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 23134cc180ffSJed Brown } 23148a381b04SJed Brown 2315d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 2316d71ae5a4SJacob Faibussowitsch { 23173a28c0e5SStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 23183a28c0e5SStefano Zampini 23193a28c0e5SStefano Zampini PetscFunctionBegin; 23203a28c0e5SStefano Zampini *flg = (PetscBool)!ark->imex; 23213ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 23223a28c0e5SStefano Zampini } 23233a28c0e5SStefano Zampini 2324d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 2325d71ae5a4SJacob Faibussowitsch { 2326b3a6b972SJed Brown PetscFunctionBegin; 23279566063dSJacob Faibussowitsch PetscCall(TSReset_ARKIMEX(ts)); 2328b3a6b972SJed Brown if (ts->dm) { 23299566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 23309566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2331b3a6b972SJed Brown } 23329566063dSJacob Faibussowitsch PetscCall(PetscFree(ts->data)); 2333d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL)); 2334d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL)); 23359566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 23369566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 23379566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 23382e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 23393a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", NULL)); 23403a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", NULL)); 23413a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", NULL)); 23423a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_ARKIMEX_FastSlowSplit_C", NULL)); 23433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2344b3a6b972SJed Brown } 2345b3a6b972SJed Brown 23468a381b04SJed Brown /* ------------------------------------------------------------ */ 23478a381b04SJed Brown /*MC 2348c79dcfbdSBarry Smith TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 23498a381b04SJed Brown 2350fca742c7SJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2351fca742c7SJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2352bcf0153eSBarry Smith of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2353d0685a90SJed Brown 23548a381b04SJed Brown Level: beginner 23558a381b04SJed Brown 2356bcf0153eSBarry Smith Notes: 2357bcf0153eSBarry Smith The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 23588a381b04SJed Brown 2359bcf0153eSBarry Smith If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2360bcf0153eSBarry Smith 2361bcf0153eSBarry Smith 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). 2362bcf0153eSBarry Smith 2363bcf0153eSBarry Smith Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2364bcf0153eSBarry Smith 23651cc06b55SBarry Smith .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2366bcf0153eSBarry Smith `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2367bcf0153eSBarry Smith `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 23688a381b04SJed Brown M*/ 2369d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2370d71ae5a4SJacob Faibussowitsch { 237180ab5e31SHong Zhang TS_ARKIMEX *ark; 2372d27334e2SStefano Zampini PetscBool dirk; 23738a381b04SJed Brown 23748a381b04SJed Brown PetscFunctionBegin; 23759566063dSJacob Faibussowitsch PetscCall(TSARKIMEXInitializePackage()); 2376d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 23778a381b04SJed Brown 23788a381b04SJed Brown ts->ops->reset = TSReset_ARKIMEX; 237980ab5e31SHong Zhang ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 23808a381b04SJed Brown ts->ops->destroy = TSDestroy_ARKIMEX; 23818a381b04SJed Brown ts->ops->view = TSView_ARKIMEX; 2382f2c2a1b9SBarry Smith ts->ops->load = TSLoad_ARKIMEX; 23838a381b04SJed Brown ts->ops->setup = TSSetUp_ARKIMEX; 238480ab5e31SHong Zhang ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 23858a381b04SJed Brown ts->ops->step = TSStep_ARKIMEX; 2386cd652676SJed Brown ts->ops->interpolate = TSInterpolate_ARKIMEX; 2387108c343cSJed Brown ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 23888a381b04SJed Brown ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 23898a381b04SJed Brown ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 23908a381b04SJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 239180ab5e31SHong Zhang ts->ops->getstages = TSGetStages_ARKIMEX; 239280ab5e31SHong Zhang ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 23938a381b04SJed Brown 2394efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 2395efd4aadfSBarry Smith 239680ab5e31SHong Zhang PetscCall(PetscNew(&ark)); 239780ab5e31SHong Zhang ts->data = (void *)ark; 2398d27334e2SStefano Zampini ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 239980ab5e31SHong Zhang 240080ab5e31SHong Zhang ark->VecsDeltaLam = NULL; 240180ab5e31SHong Zhang ark->VecsSensiTemp = NULL; 240280ab5e31SHong Zhang ark->VecsSensiPTemp = NULL; 24038a381b04SJed Brown 24049566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2405d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2406d27334e2SStefano Zampini if (!dirk) { 24079566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 24089566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 24093a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX)); 24103a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX)); 24119566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2412d27334e2SStefano Zampini } 2413d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2414d27334e2SStefano Zampini } 2415d27334e2SStefano Zampini 2416d27334e2SStefano Zampini /* ------------------------------------------------------------ */ 2417d27334e2SStefano Zampini 2418d27334e2SStefano Zampini static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2419d27334e2SStefano Zampini { 2420d27334e2SStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2421d27334e2SStefano Zampini 2422d27334e2SStefano Zampini PetscFunctionBegin; 2423d27334e2SStefano Zampini PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2424d27334e2SStefano Zampini PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2425d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2426d27334e2SStefano Zampini } 2427d27334e2SStefano Zampini 2428cc4c1da9SBarry Smith /*@ 2429d27334e2SStefano Zampini TSDIRKSetType - Set the type of `TSDIRK` scheme 2430d27334e2SStefano Zampini 2431d27334e2SStefano Zampini Logically Collective 2432d27334e2SStefano Zampini 2433d27334e2SStefano Zampini Input Parameters: 2434d27334e2SStefano Zampini + ts - timestepping context 2435d27334e2SStefano Zampini - dirktype - type of `TSDIRK` scheme 2436d27334e2SStefano Zampini 2437d27334e2SStefano Zampini Options Database Key: 2438d27334e2SStefano Zampini . -ts_dirkimex_type - set `TSDIRK` scheme type 2439d27334e2SStefano Zampini 2440d27334e2SStefano Zampini Level: intermediate 2441d27334e2SStefano Zampini 2442d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2443d27334e2SStefano Zampini @*/ 2444d27334e2SStefano Zampini PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2445d27334e2SStefano Zampini { 2446d27334e2SStefano Zampini PetscFunctionBegin; 2447d27334e2SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2448d27334e2SStefano Zampini PetscAssertPointer(dirktype, 2); 2449d27334e2SStefano Zampini PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2450d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2451d27334e2SStefano Zampini } 2452d27334e2SStefano Zampini 2453cc4c1da9SBarry Smith /*@ 2454d27334e2SStefano Zampini TSDIRKGetType - Get the type of `TSDIRK` scheme 2455d27334e2SStefano Zampini 2456d27334e2SStefano Zampini Logically Collective 2457d27334e2SStefano Zampini 2458d27334e2SStefano Zampini Input Parameter: 2459d27334e2SStefano Zampini . ts - timestepping context 2460d27334e2SStefano Zampini 2461d27334e2SStefano Zampini Output Parameter: 2462d27334e2SStefano Zampini . dirktype - type of `TSDIRK` scheme 2463d27334e2SStefano Zampini 2464d27334e2SStefano Zampini Level: intermediate 2465d27334e2SStefano Zampini 2466d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRKSetType()` 2467d27334e2SStefano Zampini @*/ 2468d27334e2SStefano Zampini PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2469d27334e2SStefano Zampini { 2470d27334e2SStefano Zampini PetscFunctionBegin; 2471d27334e2SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2472d27334e2SStefano Zampini PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2473d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2474d27334e2SStefano Zampini } 2475d27334e2SStefano Zampini 2476d27334e2SStefano Zampini /*MC 24773405e92cSStefano Zampini TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes. 2478d27334e2SStefano Zampini 2479d27334e2SStefano Zampini Level: beginner 2480d27334e2SStefano Zampini 2481d27334e2SStefano Zampini Notes: 24823405e92cSStefano Zampini The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type. 24833405e92cSStefano Zampini The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with: 24843405e92cSStefano Zampini + E - whether the method has an explicit first stage 24853405e92cSStefano Zampini . S - whether the method is single diagonal 24863405e92cSStefano Zampini . P - order of the advancing method 24873405e92cSStefano Zampini . Q - order of the embedded method 24883405e92cSStefano Zampini . S - number of stages 24893405e92cSStefano Zampini . SA - whether the method is stiffly accurate 24903405e92cSStefano Zampini . L - whether the method is L-stable 24913405e92cSStefano Zampini - A - whether the method is A-stable 2492d27334e2SStefano Zampini 2493d27334e2SStefano Zampini .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2494d27334e2SStefano Zampini M*/ 2495d27334e2SStefano Zampini PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2496d27334e2SStefano Zampini { 2497d27334e2SStefano Zampini PetscFunctionBegin; 2498d27334e2SStefano Zampini PetscCall(TSCreate_ARKIMEX(ts)); 2499d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2500d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2501d27334e2SStefano Zampini PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 25023ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 25038a381b04SJed Brown } 25043a2a065fSHong Zhang 25053a2a065fSHong Zhang /*@ 25063a2a065fSHong Zhang TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system 25073a2a065fSHong Zhang 25083a2a065fSHong Zhang Logically Collective 25093a2a065fSHong Zhang 25103a2a065fSHong Zhang Input Parameters: 25113a2a065fSHong Zhang + ts - timestepping context 25123a2a065fSHong Zhang - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise. 25133a2a065fSHong Zhang 25143a2a065fSHong Zhang Options Database Key: 25153a2a065fSHong Zhang . -ts_arkimex_fastslowsplit - <true,false> 25163a2a065fSHong Zhang 25173a2a065fSHong Zhang Level: intermediate 25183a2a065fSHong Zhang 25193a2a065fSHong Zhang .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()` 25203a2a065fSHong Zhang @*/ 25213a2a065fSHong Zhang PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow) 25223a2a065fSHong Zhang { 25233a2a065fSHong Zhang PetscFunctionBegin; 25243a2a065fSHong Zhang PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow)); 25253a2a065fSHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 25263a2a065fSHong Zhang } 25273a2a065fSHong Zhang 25283a2a065fSHong Zhang /*@ 25293a2a065fSHong Zhang TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system 25303a2a065fSHong Zhang 25313a2a065fSHong Zhang Not Collective 25323a2a065fSHong Zhang 25333a2a065fSHong Zhang Input Parameter: 25343a2a065fSHong Zhang . ts - timestepping context 25353a2a065fSHong Zhang 25363a2a065fSHong Zhang Output Parameter: 25373a2a065fSHong Zhang . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise 25383a2a065fSHong Zhang 25393a2a065fSHong Zhang Level: intermediate 25403a2a065fSHong Zhang 25413a2a065fSHong Zhang .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()` 25423a2a065fSHong Zhang @*/ 25433a2a065fSHong Zhang PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow) 25443a2a065fSHong Zhang { 25453a2a065fSHong Zhang PetscFunctionBegin; 25463a2a065fSHong Zhang PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow)); 25473a2a065fSHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 25483a2a065fSHong Zhang } 2549