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> 14*3a2a065fSHong Zhang #include <../src/ts/impls/arkimex/arkimex.h> 15*3a2a065fSHong Zhang #include <../src/ts/impls/arkimex/fsarkimex.h> 168a381b04SJed Brown 17*3a2a065fSHong 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 413*3a2a065fSHong 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; 13219566063dSJacob Faibussowitsch PetscCall(PetscInfo((PetscObject)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; 1766*3a2a065fSHong Zhang if (ark->fastslowsplit) { 1767*3a2a065fSHong Zhang PetscTryMethod(ts, "TSReset_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 1768*3a2a065fSHong 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)); 1774*3a2a065fSHong 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) { 1796*3a2a065fSHong Zhang if (dm && dm != ts->dm && !ax->fastslowsplit) { 17979566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1798d5e6173cSPeter Brune } else *Z = ax->Z; 1799d5e6173cSPeter Brune } 1800d5e6173cSPeter Brune if (Ydot) { 1801*3a2a065fSHong Zhang if (dm && dm != ts->dm && !ax->fastslowsplit) { 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 { 1810*3a2a065fSHong Zhang TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1811*3a2a065fSHong Zhang 1812d5e6173cSPeter Brune PetscFunctionBegin; 1813d5e6173cSPeter Brune if (Z) { 1814*3a2a065fSHong Zhang if (dm && dm != ts->dm && !ax->fastslowsplit) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1815d5e6173cSPeter Brune } 1816d5e6173cSPeter Brune if (Ydot) { 1817*3a2a065fSHong Zhang if (dm && dm != ts->dm && !ax->fastslowsplit) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1818d5e6173cSPeter Brune } 18193ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1820d5e6173cSPeter Brune } 1821d5e6173cSPeter Brune 18220467964bSStefano Zampini /* 18230467964bSStefano Zampini DAEs need special handling for algebraic variables when restarting DIRK methods with explicit 18240467964bSStefano Zampini first stage. In particular, we need: 18250467964bSStefano Zampini - to zero the nonlinear function (in case the dual variables are not consistent in the first step) 18260467964bSStefano Zampini - to modify the preconditioning matrix by calling MatZeroRows with identity on these variables. 18270467964bSStefano Zampini */ 18280467964bSStefano Zampini static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is) 18290467964bSStefano Zampini { 18300467964bSStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 18310467964bSStefano Zampini DM dm; 18320467964bSStefano Zampini Vec F, W, Xdot; 18330467964bSStefano Zampini const PetscScalar *w; 18340467964bSStefano Zampini PetscInt nz = 0, n, st; 18350467964bSStefano Zampini PetscInt *nzr; 18360467964bSStefano Zampini 18370467964bSStefano Zampini PetscFunctionBegin; 18380467964bSStefano Zampini PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */ 18390467964bSStefano Zampini PetscCall(DMGetGlobalVector(dm, &Xdot)); 18400467964bSStefano Zampini PetscCall(DMGetGlobalVector(dm, &F)); 18410467964bSStefano Zampini PetscCall(DMGetGlobalVector(dm, &W)); 18420467964bSStefano Zampini PetscCall(VecSet(Xdot, 0.0)); 18430467964bSStefano Zampini PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex)); 18440467964bSStefano Zampini PetscCall(VecSetRandom(Xdot, NULL)); 18450467964bSStefano Zampini PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex)); 18460467964bSStefano Zampini PetscCall(VecAXPY(W, -1.0, F)); 18470467964bSStefano Zampini PetscCall(VecGetOwnershipRange(W, &st, NULL)); 18480467964bSStefano Zampini PetscCall(VecGetLocalSize(W, &n)); 18490467964bSStefano Zampini PetscCall(VecGetArrayRead(W, &w)); 18500467964bSStefano Zampini for (PetscInt i = 0; i < n; i++) 18510467964bSStefano Zampini if (w[i] == 0.0) nz++; 18520467964bSStefano Zampini PetscCall(PetscMalloc1(nz, &nzr)); 18530467964bSStefano Zampini nz = 0; 18540467964bSStefano Zampini for (PetscInt i = 0; i < n; i++) 18550467964bSStefano Zampini if (w[i] == 0.0) nzr[nz++] = i + st; 18560467964bSStefano Zampini PetscCall(VecRestoreArrayRead(W, &w)); 18570467964bSStefano Zampini PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is)); 18580467964bSStefano Zampini PetscCall(DMRestoreGlobalVector(dm, &Xdot)); 18590467964bSStefano Zampini PetscCall(DMRestoreGlobalVector(dm, &F)); 18600467964bSStefano Zampini PetscCall(DMRestoreGlobalVector(dm, &W)); 18610467964bSStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 18620467964bSStefano Zampini } 18630467964bSStefano Zampini 18640467964bSStefano Zampini /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure 18650467964bSStefano Zampini at the finest level, in the DM for coarser solves. */ 18660467964bSStefano Zampini static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is) 18670467964bSStefano Zampini { 18680467964bSStefano Zampini TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 18690467964bSStefano Zampini 18700467964bSStefano Zampini PetscFunctionBegin; 18710467964bSStefano Zampini if (dm && dm != ts->dm) { 18720467964bSStefano Zampini PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is)); 18730467964bSStefano Zampini } else *alg_is = ax->alg_is; 18740467964bSStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 18750467964bSStefano Zampini } 18763b98415fSStefano Zampini 18773b98415fSStefano Zampini /* This defines the nonlinear equation that is to be solved with SNES */ 1878d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1879d71ae5a4SJacob Faibussowitsch { 18808a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1881d5e6173cSPeter Brune DM dm, dmsave; 1882*3a2a065fSHong Zhang Vec Z, Ydot, Y = ark->Y_snes; 18830467964bSStefano Zampini IS alg_is; 1884*3a2a065fSHong Zhang TS subts = ark->subts_fast ? ark->subts_fast : ts; 18858a381b04SJed Brown 18868a381b04SJed Brown PetscFunctionBegin; 18879566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 18889566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 18890467964bSStefano Zampini if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 18900467964bSStefano Zampini 1891d5e6173cSPeter Brune dmsave = ts->dm; 1892d5e6173cSPeter Brune ts->dm = dm; 1893740132f1SEmil Constantinescu 189498940a98SHong Zhang if (ark->scoeff == PETSC_MAX_REAL) { 18953b98415fSStefano Zampini /* We are solving F(t_n,x_n,xdot) = 0 to start the method */ 18960467964bSStefano Zampini if (!alg_is) { 18970467964bSStefano Zampini PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 18980467964bSStefano Zampini PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is)); 18990467964bSStefano Zampini PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is)); 19000467964bSStefano Zampini PetscCall(PetscObjectDereference((PetscObject)alg_is)); 19010467964bSStefano Zampini PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view")); 19020467964bSStefano Zampini } 1903*3a2a065fSHong Zhang if (ark->fastslowsplit && ark->is_slow) PetscCall(VecISCopy(Y, ark->is_fast, SCATTER_FORWARD, Z)); 1904*3a2a065fSHong Zhang else Y = Z; 1905*3a2a065fSHong Zhang PetscCall(TSComputeIFunction(subts, ark->stage_time, Y, X, F, ark->imex)); 19060467964bSStefano Zampini PetscCall(VecISSet(F, alg_is, 0.0)); 1907d27334e2SStefano Zampini } else { 190898940a98SHong Zhang PetscReal shift = ark->scoeff / ts->time_step; 1909d27334e2SStefano Zampini PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1910*3a2a065fSHong Zhang if (ark->fastslowsplit && ark->is_slow) PetscCall(VecISCopy(Y, ark->is_fast, SCATTER_FORWARD, X)); 1911*3a2a065fSHong Zhang else Y = X; 1912*3a2a065fSHong Zhang PetscCall(TSComputeIFunction(subts, ark->stage_time, Y, Ydot, F, ark->imex)); 1913d27334e2SStefano Zampini } 1914e817cc15SEmil Constantinescu 1915d5e6173cSPeter Brune ts->dm = dmsave; 19169566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 19173ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19188a381b04SJed Brown } 19198a381b04SJed Brown 1920d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1921d71ae5a4SJacob Faibussowitsch { 19228a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1923d5e6173cSPeter Brune DM dm, dmsave; 1924*3a2a065fSHong Zhang Vec Ydot, Z, Y = ark->Y_snes; 192598940a98SHong Zhang PetscReal shift; 19260467964bSStefano Zampini IS alg_is; 1927*3a2a065fSHong Zhang TS subts = ark->subts_fast ? ark->subts_fast : ts; 19288a381b04SJed Brown 19298a381b04SJed Brown PetscFunctionBegin; 19309566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 19318a381b04SJed Brown /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 19320467964bSStefano Zampini PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 19330467964bSStefano Zampini /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */ 19340467964bSStefano Zampini if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 19350467964bSStefano Zampini 1936d5e6173cSPeter Brune dmsave = ts->dm; 1937d5e6173cSPeter Brune ts->dm = dm; 1938740132f1SEmil Constantinescu 193998940a98SHong Zhang if (ark->scoeff == PETSC_MAX_REAL) { 19403b98415fSStefano Zampini PetscBool hasZeroRows; 19413b98415fSStefano Zampini 19423b98415fSStefano Zampini /* We are solving F(t_n,x_n,xdot) = 0 to start the method 19430467964bSStefano Zampini We compute with a very large shift and then scale back the matrix */ 1944d27334e2SStefano Zampini shift = 1.0 / PETSC_MACHINE_EPSILON; 1945*3a2a065fSHong Zhang if (ark->fastslowsplit && ark->is_slow) PetscCall(VecISCopy(Y, ark->is_fast, SCATTER_FORWARD, Z)); 1946*3a2a065fSHong Zhang else Y = Z; 1947*3a2a065fSHong Zhang PetscCall(TSComputeIJacobian(subts, ark->stage_time, Y, X, shift, A, B, ark->imex)); 1948d27334e2SStefano Zampini PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 19493b98415fSStefano Zampini PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows)); 19503b98415fSStefano Zampini if (hasZeroRows) { 19510467964bSStefano Zampini PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 19523b98415fSStefano Zampini /* the default of AIJ is to not keep the pattern! We should probably change it someday */ 19533b98415fSStefano Zampini PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 19543b98415fSStefano Zampini PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL)); 19553b98415fSStefano Zampini } 19563b98415fSStefano Zampini PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view")); 1957d27334e2SStefano Zampini if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1958d27334e2SStefano Zampini } else { 195998940a98SHong Zhang shift = ark->scoeff / ts->time_step; 1960*3a2a065fSHong Zhang if (ark->fastslowsplit && ark->is_slow) PetscCall(VecISCopy(Y, ark->is_fast, SCATTER_FORWARD, X)); 1961*3a2a065fSHong Zhang else Y = X; 1962*3a2a065fSHong Zhang PetscCall(TSComputeIJacobian(subts, ark->stage_time, Y, Ydot, shift, A, B, ark->imex)); 1963d27334e2SStefano Zampini } 1964d5e6173cSPeter Brune ts->dm = dmsave; 1965d27334e2SStefano Zampini PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 19663ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1967d5e6173cSPeter Brune } 1968d5e6173cSPeter Brune 196980ab5e31SHong Zhang static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 197080ab5e31SHong Zhang { 197180ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 197280ab5e31SHong Zhang 197380ab5e31SHong Zhang PetscFunctionBegin; 197480ab5e31SHong Zhang if (ns) *ns = ark->tableau->s; 197580ab5e31SHong Zhang if (Y) *Y = ark->Y; 197680ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 197780ab5e31SHong Zhang } 197880ab5e31SHong Zhang 1979d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1980d71ae5a4SJacob Faibussowitsch { 1981d5e6173cSPeter Brune PetscFunctionBegin; 19823ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1983d5e6173cSPeter Brune } 1984d5e6173cSPeter Brune 1985d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1986d71ae5a4SJacob Faibussowitsch { 1987d5e6173cSPeter Brune TS ts = (TS)ctx; 1988d5e6173cSPeter Brune Vec Z, Z_c; 1989d5e6173cSPeter Brune 1990d5e6173cSPeter Brune PetscFunctionBegin; 19919566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 19929566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 19939566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Z, Z_c)); 19949566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 19959566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 19969566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 19973ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19988a381b04SJed Brown } 19998a381b04SJed Brown 2000d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 2001d71ae5a4SJacob Faibussowitsch { 2002cdb298fcSPeter Brune PetscFunctionBegin; 20033ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2004cdb298fcSPeter Brune } 2005cdb298fcSPeter Brune 2006d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 2007d71ae5a4SJacob Faibussowitsch { 2008cdb298fcSPeter Brune TS ts = (TS)ctx; 2009cdb298fcSPeter Brune Vec Z, Z_c; 2010cdb298fcSPeter Brune 2011cdb298fcSPeter Brune PetscFunctionBegin; 20129566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 20139566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 2014cdb298fcSPeter Brune 20159566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 20169566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2017cdb298fcSPeter Brune 20189566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 20199566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 20203ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2021cdb298fcSPeter Brune } 2022cdb298fcSPeter Brune 2023d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 2024d71ae5a4SJacob Faibussowitsch { 202596400bd6SLisandro Dalcin TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 202696400bd6SLisandro Dalcin ARKTableau tab = ark->tableau; 202796400bd6SLisandro Dalcin 202896400bd6SLisandro Dalcin PetscFunctionBegin; 2029d27334e2SStefano Zampini PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 20309566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 20319566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 2032d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 203396400bd6SLisandro Dalcin if (ark->extrapolate) { 20349566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 20359566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 2036d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 203796400bd6SLisandro Dalcin } 20383ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 203996400bd6SLisandro Dalcin } 204096400bd6SLisandro Dalcin 2041d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 2042d71ae5a4SJacob Faibussowitsch { 20438a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2044d5e6173cSPeter Brune DM dm; 204596400bd6SLisandro Dalcin SNES snes; 2046f9c1d6abSBarry Smith 20478a381b04SJed Brown PetscFunctionBegin; 2048*3a2a065fSHong Zhang if (ark->fastslowsplit) { 2049*3a2a065fSHong Zhang PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 2050*3a2a065fSHong Zhang } else { 20519566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTableauSetUp(ts)); 20529566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 20539566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 20549566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 20559566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm)); 20569566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 20579566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 20589566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 2059*3a2a065fSHong Zhang PetscCall(SNESSetDM(snes, dm)); 2060*3a2a065fSHong Zhang } 20613ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20628a381b04SJed Brown } 20638a381b04SJed Brown 206480ab5e31SHong Zhang static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 206580ab5e31SHong Zhang { 206680ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 206780ab5e31SHong Zhang ARKTableau tab = ark->tableau; 206880ab5e31SHong Zhang 206980ab5e31SHong Zhang PetscFunctionBegin; 207080ab5e31SHong Zhang PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 207180ab5e31SHong Zhang PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 207280ab5e31SHong Zhang if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); } 207380ab5e31SHong Zhang if (PetscDefined(USE_DEBUG)) { 207480ab5e31SHong Zhang PetscBool id = PETSC_FALSE; 207580ab5e31SHong Zhang PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 20768434afd1SBarry 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"); 207780ab5e31SHong Zhang } 207880ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 207980ab5e31SHong Zhang } 208080ab5e31SHong Zhang 2081d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject) 2082d71ae5a4SJacob Faibussowitsch { 20834cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2084d27334e2SStefano Zampini PetscBool dirk; 20858a381b04SJed Brown 20868a381b04SJed Brown PetscFunctionBegin; 2087d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2088d27334e2SStefano Zampini PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 20898a381b04SJed Brown { 20908a381b04SJed Brown ARKTableauLink link; 20918a381b04SJed Brown PetscInt count, choice; 20928a381b04SJed Brown PetscBool flg; 20938a381b04SJed Brown const char **namelist; 2094d27334e2SStefano Zampini for (link = ARKTableauList, count = 0; link; link = link->next) { 2095d27334e2SStefano Zampini if (!dirk && link->tab.additive) count++; 2096d27334e2SStefano Zampini if (dirk && !link->tab.additive) count++; 2097d27334e2SStefano Zampini } 20989566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(count, (char ***)&namelist)); 2099d27334e2SStefano Zampini for (link = ARKTableauList, count = 0; link; link = link->next) { 2100d27334e2SStefano Zampini if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 2101d27334e2SStefano Zampini if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 2102d27334e2SStefano Zampini } 2103d27334e2SStefano Zampini if (dirk) { 2104d27334e2SStefano Zampini PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2105d27334e2SStefano Zampini if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 2106d27334e2SStefano Zampini } else { 2107*3a2a065fSHong Zhang PetscBool fastslowsplit; 21089566063dSJacob Faibussowitsch PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 21099566063dSJacob Faibussowitsch if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 21104cc180ffSJed Brown flg = (PetscBool)!ark->imex; 21119566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 21124cc180ffSJed Brown ark->imex = (PetscBool)!flg; 2113*3a2a065fSHong Zhang PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg)); 2114*3a2a065fSHong Zhang if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit)); 2115d27334e2SStefano Zampini } 2116d27334e2SStefano Zampini PetscCall(PetscFree(namelist)); 21179566063dSJacob 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)); 21188a381b04SJed Brown } 2119d0609cedSBarry Smith PetscOptionsHeadEnd(); 21203ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 21218a381b04SJed Brown } 21228a381b04SJed Brown 2123d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 2124d71ae5a4SJacob Faibussowitsch { 21258a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2126d27334e2SStefano Zampini PetscBool iascii, dirk; 21278a381b04SJed Brown 21288a381b04SJed Brown PetscFunctionBegin; 2129d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 21309566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 21318a381b04SJed Brown if (iascii) { 2132d27334e2SStefano Zampini PetscViewerFormat format; 21339c334d8fSLisandro Dalcin ARKTableau tab = ark->tableau; 213419fd82e9SBarry Smith TSARKIMEXType arktype; 2135d27334e2SStefano Zampini char buf[2048]; 21363a28c0e5SStefano Zampini PetscBool flg; 21373a28c0e5SStefano Zampini 21389566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetType(ts, &arktype)); 21399566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 2140d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 21419566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 2142d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 2143d27334e2SStefano Zampini PetscCall(PetscViewerGetFormat(viewer, &format)); 2144d27334e2SStefano Zampini if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2145d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 2146d27334e2SStefano Zampini for (PetscInt i = 0; i < tab->s; i++) { 2147d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 2148d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 2149d27334e2SStefano Zampini } 2150d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 2151d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 2152d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 2153d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 2154d27334e2SStefano Zampini } 21559566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 21569566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 21579566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 21589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 2159d27334e2SStefano Zampini if (!dirk) { 2160d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 21619566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 21628a381b04SJed Brown } 2163d27334e2SStefano Zampini } 21643ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 21658a381b04SJed Brown } 21668a381b04SJed Brown 2167d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 2168d71ae5a4SJacob Faibussowitsch { 2169f2c2a1b9SBarry Smith SNES snes; 21709c334d8fSLisandro Dalcin TSAdapt adapt; 2171f2c2a1b9SBarry Smith 2172f2c2a1b9SBarry Smith PetscFunctionBegin; 21739566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 21749566063dSJacob Faibussowitsch PetscCall(TSAdaptLoad(adapt, viewer)); 21759566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 21769566063dSJacob Faibussowitsch PetscCall(SNESLoad(snes, viewer)); 2177ad6bc421SBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 21789566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 21799566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 21803ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2181f2c2a1b9SBarry Smith } 2182f2c2a1b9SBarry Smith 2183cc4c1da9SBarry Smith /*@ 2184bcf0153eSBarry Smith TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 21858a381b04SJed Brown 218620f4b53cSBarry Smith Logically Collective 21878a381b04SJed Brown 2188d8d19677SJose E. Roman Input Parameters: 21898a381b04SJed Brown + ts - timestepping context 2190bcf0153eSBarry Smith - arktype - type of `TSARKIMEX` scheme 21918a381b04SJed Brown 2192bcf0153eSBarry Smith Options Database Key: 2193bcf0153eSBarry Smith . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 21949bd3e852SBarry Smith 21958a381b04SJed Brown Level: intermediate 21968a381b04SJed Brown 21971cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 2198db781477SPatrick Sanan `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 21998a381b04SJed Brown @*/ 2200d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 2201d71ae5a4SJacob Faibussowitsch { 22028a381b04SJed Brown PetscFunctionBegin; 22038a381b04SJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 22044f572ea9SToby Isaac PetscAssertPointer(arktype, 2); 2205cac4c232SBarry Smith PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 22063ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22078a381b04SJed Brown } 22088a381b04SJed Brown 2209cc4c1da9SBarry Smith /*@ 2210bcf0153eSBarry Smith TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 22118a381b04SJed Brown 221220f4b53cSBarry Smith Logically Collective 22138a381b04SJed Brown 22148a381b04SJed Brown Input Parameter: 22158a381b04SJed Brown . ts - timestepping context 22168a381b04SJed Brown 22178a381b04SJed Brown Output Parameter: 2218bcf0153eSBarry Smith . arktype - type of `TSARKIMEX` scheme 22198a381b04SJed Brown 22208a381b04SJed Brown Level: intermediate 22218a381b04SJed Brown 222242747ad1SJacob Faibussowitsch .seealso: [](ch_ts), `TSARKIMEXc` 22238a381b04SJed Brown @*/ 2224d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2225d71ae5a4SJacob Faibussowitsch { 22268a381b04SJed Brown PetscFunctionBegin; 22278a381b04SJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2228cac4c232SBarry Smith PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 22293ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22308a381b04SJed Brown } 22318a381b04SJed Brown 223216353aafSBarry Smith /*@ 2233bcf0153eSBarry Smith TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 22344cc180ffSJed Brown 223520f4b53cSBarry Smith Logically Collective 22364cc180ffSJed Brown 2237d8d19677SJose E. Roman Input Parameters: 22384cc180ffSJed Brown + ts - timestepping context 2239bcf0153eSBarry Smith - flg - `PETSC_TRUE` for fully implicit 22404cc180ffSJed Brown 22414cc180ffSJed Brown Level: intermediate 22424cc180ffSJed Brown 22431cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 22444cc180ffSJed Brown @*/ 2245d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 2246d71ae5a4SJacob Faibussowitsch { 22474cc180ffSJed Brown PetscFunctionBegin; 22484cc180ffSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 22493a28c0e5SStefano Zampini PetscValidLogicalCollectiveBool(ts, flg, 2); 2250cac4c232SBarry Smith PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 22513ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22524cc180ffSJed Brown } 22534cc180ffSJed Brown 22543a28c0e5SStefano Zampini /*@ 22553a28c0e5SStefano Zampini TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 22563a28c0e5SStefano Zampini 225720f4b53cSBarry Smith Logically Collective 22583a28c0e5SStefano Zampini 22593a28c0e5SStefano Zampini Input Parameter: 22603a28c0e5SStefano Zampini . ts - timestepping context 22613a28c0e5SStefano Zampini 22627a7aea1fSJed Brown Output Parameter: 2263bcf0153eSBarry Smith . flg - `PETSC_TRUE` for fully implicit 22643a28c0e5SStefano Zampini 22653a28c0e5SStefano Zampini Level: intermediate 22663a28c0e5SStefano Zampini 22671cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 22683a28c0e5SStefano Zampini @*/ 2269d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 2270d71ae5a4SJacob Faibussowitsch { 22713a28c0e5SStefano Zampini PetscFunctionBegin; 22723a28c0e5SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 22734f572ea9SToby Isaac PetscAssertPointer(flg, 2); 2274cac4c232SBarry Smith PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 22753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22763a28c0e5SStefano Zampini } 22773a28c0e5SStefano Zampini 2278d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 2279d71ae5a4SJacob Faibussowitsch { 22808a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 22818a381b04SJed Brown 22828a381b04SJed Brown PetscFunctionBegin; 22838a381b04SJed Brown *arktype = ark->tableau->name; 22843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22858a381b04SJed Brown } 2286d27334e2SStefano Zampini 2287d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 2288d71ae5a4SJacob Faibussowitsch { 22898a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 22908a381b04SJed Brown PetscBool match; 22918a381b04SJed Brown ARKTableauLink link; 22928a381b04SJed Brown 22938a381b04SJed Brown PetscFunctionBegin; 22948a381b04SJed Brown if (ark->tableau) { 22959566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 22963ba16761SJacob Faibussowitsch if (match) PetscFunctionReturn(PETSC_SUCCESS); 22978a381b04SJed Brown } 22988a381b04SJed Brown for (link = ARKTableauList; link; link = link->next) { 22999566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 23008a381b04SJed Brown if (match) { 23019566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 23028a381b04SJed Brown ark->tableau = &link->tab; 23039566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 23042ffb9264SLisandro Dalcin ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 23053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 23068a381b04SJed Brown } 23078a381b04SJed Brown } 230898921bdaSJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 23098a381b04SJed Brown } 2310e0877f53SBarry Smith 2311d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 2312d71ae5a4SJacob Faibussowitsch { 23134cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 23144cc180ffSJed Brown 23154cc180ffSJed Brown PetscFunctionBegin; 23164cc180ffSJed Brown ark->imex = (PetscBool)!flg; 23173ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 23184cc180ffSJed Brown } 23198a381b04SJed Brown 2320d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 2321d71ae5a4SJacob Faibussowitsch { 23223a28c0e5SStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 23233a28c0e5SStefano Zampini 23243a28c0e5SStefano Zampini PetscFunctionBegin; 23253a28c0e5SStefano Zampini *flg = (PetscBool)!ark->imex; 23263ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 23273a28c0e5SStefano Zampini } 23283a28c0e5SStefano Zampini 2329d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 2330d71ae5a4SJacob Faibussowitsch { 2331b3a6b972SJed Brown PetscFunctionBegin; 23329566063dSJacob Faibussowitsch PetscCall(TSReset_ARKIMEX(ts)); 2333b3a6b972SJed Brown if (ts->dm) { 23349566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 23359566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2336b3a6b972SJed Brown } 23379566063dSJacob Faibussowitsch PetscCall(PetscFree(ts->data)); 2338d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL)); 2339d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL)); 23409566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 23419566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 23429566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 23432e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 2344*3a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", NULL)); 2345*3a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", NULL)); 2346*3a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", NULL)); 2347*3a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_ARKIMEX_FastSlowSplit_C", NULL)); 23483ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2349b3a6b972SJed Brown } 2350b3a6b972SJed Brown 23518a381b04SJed Brown /* ------------------------------------------------------------ */ 23528a381b04SJed Brown /*MC 2353c79dcfbdSBarry Smith TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 23548a381b04SJed Brown 2355fca742c7SJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2356fca742c7SJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2357bcf0153eSBarry Smith of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2358d0685a90SJed Brown 23598a381b04SJed Brown Level: beginner 23608a381b04SJed Brown 2361bcf0153eSBarry Smith Notes: 2362bcf0153eSBarry Smith The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 23638a381b04SJed Brown 2364bcf0153eSBarry Smith If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2365bcf0153eSBarry Smith 2366bcf0153eSBarry 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). 2367bcf0153eSBarry Smith 2368bcf0153eSBarry Smith Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2369bcf0153eSBarry Smith 23701cc06b55SBarry Smith .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2371bcf0153eSBarry Smith `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2372bcf0153eSBarry Smith `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 23738a381b04SJed Brown M*/ 2374d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2375d71ae5a4SJacob Faibussowitsch { 237680ab5e31SHong Zhang TS_ARKIMEX *ark; 2377d27334e2SStefano Zampini PetscBool dirk; 23788a381b04SJed Brown 23798a381b04SJed Brown PetscFunctionBegin; 23809566063dSJacob Faibussowitsch PetscCall(TSARKIMEXInitializePackage()); 2381d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 23828a381b04SJed Brown 23838a381b04SJed Brown ts->ops->reset = TSReset_ARKIMEX; 238480ab5e31SHong Zhang ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 23858a381b04SJed Brown ts->ops->destroy = TSDestroy_ARKIMEX; 23868a381b04SJed Brown ts->ops->view = TSView_ARKIMEX; 2387f2c2a1b9SBarry Smith ts->ops->load = TSLoad_ARKIMEX; 23888a381b04SJed Brown ts->ops->setup = TSSetUp_ARKIMEX; 238980ab5e31SHong Zhang ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 23908a381b04SJed Brown ts->ops->step = TSStep_ARKIMEX; 2391cd652676SJed Brown ts->ops->interpolate = TSInterpolate_ARKIMEX; 2392108c343cSJed Brown ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 23938a381b04SJed Brown ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 23948a381b04SJed Brown ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 23958a381b04SJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 239680ab5e31SHong Zhang ts->ops->getstages = TSGetStages_ARKIMEX; 239780ab5e31SHong Zhang ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 23988a381b04SJed Brown 2399efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 2400efd4aadfSBarry Smith 240180ab5e31SHong Zhang PetscCall(PetscNew(&ark)); 240280ab5e31SHong Zhang ts->data = (void *)ark; 2403d27334e2SStefano Zampini ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 240480ab5e31SHong Zhang 240580ab5e31SHong Zhang ark->VecsDeltaLam = NULL; 240680ab5e31SHong Zhang ark->VecsSensiTemp = NULL; 240780ab5e31SHong Zhang ark->VecsSensiPTemp = NULL; 24088a381b04SJed Brown 24099566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2410d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2411d27334e2SStefano Zampini if (!dirk) { 24129566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 24139566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 2414*3a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX)); 2415*3a2a065fSHong Zhang PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX)); 24169566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2417d27334e2SStefano Zampini } 2418d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2419d27334e2SStefano Zampini } 2420d27334e2SStefano Zampini 2421d27334e2SStefano Zampini /* ------------------------------------------------------------ */ 2422d27334e2SStefano Zampini 2423d27334e2SStefano Zampini static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2424d27334e2SStefano Zampini { 2425d27334e2SStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2426d27334e2SStefano Zampini 2427d27334e2SStefano Zampini PetscFunctionBegin; 2428d27334e2SStefano Zampini PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2429d27334e2SStefano Zampini PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2430d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2431d27334e2SStefano Zampini } 2432d27334e2SStefano Zampini 2433cc4c1da9SBarry Smith /*@ 2434d27334e2SStefano Zampini TSDIRKSetType - Set the type of `TSDIRK` scheme 2435d27334e2SStefano Zampini 2436d27334e2SStefano Zampini Logically Collective 2437d27334e2SStefano Zampini 2438d27334e2SStefano Zampini Input Parameters: 2439d27334e2SStefano Zampini + ts - timestepping context 2440d27334e2SStefano Zampini - dirktype - type of `TSDIRK` scheme 2441d27334e2SStefano Zampini 2442d27334e2SStefano Zampini Options Database Key: 2443d27334e2SStefano Zampini . -ts_dirkimex_type - set `TSDIRK` scheme type 2444d27334e2SStefano Zampini 2445d27334e2SStefano Zampini Level: intermediate 2446d27334e2SStefano Zampini 2447d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2448d27334e2SStefano Zampini @*/ 2449d27334e2SStefano Zampini PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2450d27334e2SStefano Zampini { 2451d27334e2SStefano Zampini PetscFunctionBegin; 2452d27334e2SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2453d27334e2SStefano Zampini PetscAssertPointer(dirktype, 2); 2454d27334e2SStefano Zampini PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2455d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2456d27334e2SStefano Zampini } 2457d27334e2SStefano Zampini 2458cc4c1da9SBarry Smith /*@ 2459d27334e2SStefano Zampini TSDIRKGetType - Get the type of `TSDIRK` scheme 2460d27334e2SStefano Zampini 2461d27334e2SStefano Zampini Logically Collective 2462d27334e2SStefano Zampini 2463d27334e2SStefano Zampini Input Parameter: 2464d27334e2SStefano Zampini . ts - timestepping context 2465d27334e2SStefano Zampini 2466d27334e2SStefano Zampini Output Parameter: 2467d27334e2SStefano Zampini . dirktype - type of `TSDIRK` scheme 2468d27334e2SStefano Zampini 2469d27334e2SStefano Zampini Level: intermediate 2470d27334e2SStefano Zampini 2471d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRKSetType()` 2472d27334e2SStefano Zampini @*/ 2473d27334e2SStefano Zampini PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2474d27334e2SStefano Zampini { 2475d27334e2SStefano Zampini PetscFunctionBegin; 2476d27334e2SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2477d27334e2SStefano Zampini PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2478d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2479d27334e2SStefano Zampini } 2480d27334e2SStefano Zampini 2481d27334e2SStefano Zampini /*MC 24823405e92cSStefano Zampini TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes. 2483d27334e2SStefano Zampini 2484d27334e2SStefano Zampini Level: beginner 2485d27334e2SStefano Zampini 2486d27334e2SStefano Zampini Notes: 24873405e92cSStefano Zampini The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type. 24883405e92cSStefano Zampini The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with: 24893405e92cSStefano Zampini + E - whether the method has an explicit first stage 24903405e92cSStefano Zampini . S - whether the method is single diagonal 24913405e92cSStefano Zampini . P - order of the advancing method 24923405e92cSStefano Zampini . Q - order of the embedded method 24933405e92cSStefano Zampini . S - number of stages 24943405e92cSStefano Zampini . SA - whether the method is stiffly accurate 24953405e92cSStefano Zampini . L - whether the method is L-stable 24963405e92cSStefano Zampini - A - whether the method is A-stable 2497d27334e2SStefano Zampini 2498d27334e2SStefano Zampini .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2499d27334e2SStefano Zampini M*/ 2500d27334e2SStefano Zampini PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2501d27334e2SStefano Zampini { 2502d27334e2SStefano Zampini PetscFunctionBegin; 2503d27334e2SStefano Zampini PetscCall(TSCreate_ARKIMEX(ts)); 2504d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2505d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2506d27334e2SStefano Zampini PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 25073ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 25088a381b04SJed Brown } 2509*3a2a065fSHong Zhang 2510*3a2a065fSHong Zhang /*@ 2511*3a2a065fSHong Zhang TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system 2512*3a2a065fSHong Zhang 2513*3a2a065fSHong Zhang Logically Collective 2514*3a2a065fSHong Zhang 2515*3a2a065fSHong Zhang Input Parameters: 2516*3a2a065fSHong Zhang + ts - timestepping context 2517*3a2a065fSHong Zhang - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise. 2518*3a2a065fSHong Zhang 2519*3a2a065fSHong Zhang Options Database Key: 2520*3a2a065fSHong Zhang . -ts_arkimex_fastslowsplit - <true,false> 2521*3a2a065fSHong Zhang 2522*3a2a065fSHong Zhang Level: intermediate 2523*3a2a065fSHong Zhang 2524*3a2a065fSHong Zhang .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()` 2525*3a2a065fSHong Zhang @*/ 2526*3a2a065fSHong Zhang PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow) 2527*3a2a065fSHong Zhang { 2528*3a2a065fSHong Zhang PetscFunctionBegin; 2529*3a2a065fSHong Zhang PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow)); 2530*3a2a065fSHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 2531*3a2a065fSHong Zhang } 2532*3a2a065fSHong Zhang 2533*3a2a065fSHong Zhang /*@ 2534*3a2a065fSHong Zhang TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system 2535*3a2a065fSHong Zhang 2536*3a2a065fSHong Zhang Not Collective 2537*3a2a065fSHong Zhang 2538*3a2a065fSHong Zhang Input Parameter: 2539*3a2a065fSHong Zhang . ts - timestepping context 2540*3a2a065fSHong Zhang 2541*3a2a065fSHong Zhang Output Parameter: 2542*3a2a065fSHong Zhang . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise 2543*3a2a065fSHong Zhang 2544*3a2a065fSHong Zhang Level: intermediate 2545*3a2a065fSHong Zhang 2546*3a2a065fSHong Zhang .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()` 2547*3a2a065fSHong Zhang @*/ 2548*3a2a065fSHong Zhang PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow) 2549*3a2a065fSHong Zhang { 2550*3a2a065fSHong Zhang PetscFunctionBegin; 2551*3a2a065fSHong Zhang PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow)); 2552*3a2a065fSHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 2553*3a2a065fSHong Zhang } 2554