18a381b04SJed Brown /* 2*d27334e2SStefano Zampini Code for timestepping with additive Runge-Kutta IMEX method or Diagonally Implicit Runge-Kutta methods. 38a381b04SJed Brown 48a381b04SJed Brown Notes: 5*d27334e2SStefano 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> 148a381b04SJed Brown 1519fd82e9SBarry Smith static TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3; 16*d27334e2SStefano Zampini static TSDIRKType TSDIRKDefault = TSDIRKS212; 178a381b04SJed Brown static PetscBool TSARKIMEXRegisterAllCalled; 188a381b04SJed Brown static PetscBool TSARKIMEXPackageInitialized; 1956dcabbaSDebojyoti Ghosh static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec); 208a381b04SJed Brown 218a381b04SJed Brown typedef struct _ARKTableau *ARKTableau; 228a381b04SJed Brown struct _ARKTableau { 238a381b04SJed Brown char *name; 24*d27334e2SStefano Zampini PetscBool additive; /* If False, it is a DIRK method */ 254f385281SJed Brown PetscInt order; /* Classical approximation order of the method */ 264f385281SJed Brown PetscInt s; /* Number of stages */ 27e817cc15SEmil Constantinescu PetscBool stiffly_accurate; /* The implicit part is stiffly accurate */ 28e817cc15SEmil Constantinescu PetscBool FSAL_implicit; /* The implicit part is FSAL */ 29e817cc15SEmil Constantinescu PetscBool explicit_first_stage; /* The implicit part has an explicit first stage */ 304f385281SJed Brown PetscInt pinterp; /* Interpolation order */ 314f385281SJed Brown PetscReal *At, *bt, *ct; /* Stiff tableau */ 328a381b04SJed Brown PetscReal *A, *b, *c; /* Non-stiff tableau */ 33108c343cSJed Brown PetscReal *bembedt, *bembed; /* Embedded formula of order one less (order-1) */ 34cd652676SJed Brown PetscReal *binterpt, *binterp; /* Dense output formula */ 35108c343cSJed Brown PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 368a381b04SJed Brown }; 378a381b04SJed Brown typedef struct _ARKTableauLink *ARKTableauLink; 388a381b04SJed Brown struct _ARKTableauLink { 398a381b04SJed Brown struct _ARKTableau tab; 408a381b04SJed Brown ARKTableauLink next; 418a381b04SJed Brown }; 428a381b04SJed Brown static ARKTableauLink ARKTableauList; 438a381b04SJed Brown 448a381b04SJed Brown typedef struct { 458a381b04SJed Brown ARKTableau tableau; 468a381b04SJed Brown Vec *Y; /* States computed during the step */ 478a381b04SJed Brown Vec *YdotI; /* Time derivatives for the stiff part */ 488a381b04SJed Brown Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 4956dcabbaSDebojyoti Ghosh Vec *Y_prev; /* States computed during the previous time step */ 5056dcabbaSDebojyoti Ghosh Vec *YdotI_prev; /* Time derivatives for the stiff part for the previous time step*/ 5156dcabbaSDebojyoti Ghosh Vec *YdotRHS_prev; /* Function evaluations for the non-stiff part for the previous time step*/ 52e817cc15SEmil Constantinescu Vec Ydot0; /* Holds the slope from the previous step in FSAL case */ 538a381b04SJed Brown Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 548a381b04SJed Brown Vec Z; /* Ydot = shift(Y-Z) */ 558a381b04SJed Brown PetscScalar *work; /* Scalar work */ 56b296d7d5SJed Brown PetscReal scoeff; /* shift = scoeff/dt */ 578a381b04SJed Brown PetscReal stage_time; 584cc180ffSJed Brown PetscBool imex; 5996400bd6SLisandro Dalcin PetscBool extrapolate; /* Extrapolate initial guess from previous time-step stage values */ 60108c343cSJed Brown TSStepStatus status; 6180ab5e31SHong Zhang 6280ab5e31SHong Zhang /* context for sensitivity analysis */ 6380ab5e31SHong Zhang Vec *VecsDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */ 6480ab5e31SHong Zhang Vec *VecsSensiTemp; /* Vectors to be multiplied with Jacobian transpose */ 6580ab5e31SHong Zhang Vec *VecsSensiPTemp; /* Temporary Vectors to store JacobianP-transpose-vector product */ 668a381b04SJed Brown } TS_ARKIMEX; 67*d27334e2SStefano Zampini 681f80e275SEmil Constantinescu /*MC 691f80e275SEmil Constantinescu TSARKIMEXARS122 - Second order ARK IMEX scheme. 708a381b04SJed Brown 711f80e275SEmil Constantinescu This method has one explicit stage and one implicit stage. 721f80e275SEmil Constantinescu 73bcf0153eSBarry Smith Options Database Key: 7467b8a455SSatish Balay . -ts_arkimex_type ars122 - set arkimex type to ars122 759bd3e852SBarry Smith 76bcf0153eSBarry Smith Level: advanced 77bcf0153eSBarry Smith 781f80e275SEmil Constantinescu References: 79606c0280SSatish Balay . * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 801f80e275SEmil Constantinescu 811cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 821f80e275SEmil Constantinescu M*/ 83*d27334e2SStefano Zampini 841f80e275SEmil Constantinescu /*MC 851f80e275SEmil Constantinescu TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part. 861f80e275SEmil Constantinescu 871f80e275SEmil 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. 881f80e275SEmil Constantinescu 89bcf0153eSBarry Smith Options Database Key: 9067b8a455SSatish Balay . -ts_arkimex_type a2 - set arkimex type to a2 919bd3e852SBarry Smith 921f80e275SEmil Constantinescu Level: advanced 931f80e275SEmil Constantinescu 941cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 951f80e275SEmil Constantinescu M*/ 96*d27334e2SStefano Zampini 971f80e275SEmil Constantinescu /*MC 981f80e275SEmil Constantinescu TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part. 991f80e275SEmil Constantinescu 1001f80e275SEmil Constantinescu This method has two implicit stages, and L-stable implicit scheme. 1011f80e275SEmil Constantinescu 102bcf0153eSBarry Smith Options Database Key: 10367b8a455SSatish Balay . -ts_arkimex_type l2 - set arkimex type to l2 1049bd3e852SBarry Smith 105bcf0153eSBarry Smith Level: advanced 106bcf0153eSBarry Smith 1071f80e275SEmil Constantinescu References: 108606c0280SSatish Balay . * - L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005. 1091f80e275SEmil Constantinescu 1101cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 1111f80e275SEmil Constantinescu M*/ 112*d27334e2SStefano Zampini 1131f80e275SEmil Constantinescu /*MC 114c79dcfbdSBarry Smith TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method. 115e817cc15SEmil Constantinescu 116e817cc15SEmil Constantinescu This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used. 117e817cc15SEmil Constantinescu 118bcf0153eSBarry Smith Options Database Key: 11967b8a455SSatish Balay . -ts_arkimex_type 1bee - set arkimex type to 1bee 1209bd3e852SBarry Smith 121e817cc15SEmil Constantinescu Level: advanced 122e817cc15SEmil Constantinescu 1231cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 124e817cc15SEmil Constantinescu M*/ 125*d27334e2SStefano Zampini 126e817cc15SEmil Constantinescu /*MC 1271f80e275SEmil Constantinescu TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part. 1281f80e275SEmil Constantinescu 1291f80e275SEmil 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. 1301f80e275SEmil Constantinescu 131bcf0153eSBarry Smith Options Database Key: 13267b8a455SSatish Balay . -ts_arkimex_type 2c - set arkimex type to 2c 1339bd3e852SBarry Smith 1341f80e275SEmil Constantinescu Level: advanced 1351f80e275SEmil Constantinescu 1361cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 1371f80e275SEmil Constantinescu M*/ 138*d27334e2SStefano Zampini 13964f491ddSJed Brown /*MC 14064f491ddSJed Brown TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 14164f491ddSJed Brown 142da81f932SPierre 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. 14364f491ddSJed Brown 144bcf0153eSBarry Smith Options Database Key: 14567b8a455SSatish Balay . -ts_arkimex_type 2d - set arkimex type to 2d 1469bd3e852SBarry Smith 147b330ce4dSSatish Balay Level: advanced 148b330ce4dSSatish Balay 1491cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 15064f491ddSJed Brown M*/ 151*d27334e2SStefano Zampini 15264f491ddSJed Brown /*MC 15364f491ddSJed Brown TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 15464f491ddSJed Brown 15564f491ddSJed Brown This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu. 15664f491ddSJed Brown 157bcf0153eSBarry Smith Options Database Key: 15867b8a455SSatish Balay . -ts_arkimex_type 2e - set arkimex type to 2e 1599bd3e852SBarry Smith 160b330ce4dSSatish Balay Level: advanced 161b330ce4dSSatish Balay 1621cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 16364f491ddSJed Brown M*/ 164*d27334e2SStefano Zampini 16564f491ddSJed Brown /*MC 1666cf0794eSJed Brown TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme. 1676cf0794eSJed Brown 1686cf0794eSJed Brown This method has three implicit stages. 1696cf0794eSJed Brown 1706cf0794eSJed Brown References: 171606c0280SSatish Balay . * - L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005. 1726cf0794eSJed Brown 173a8d69d7bSBarry Smith This method is referred to as SSP2-(3,3,2) in https://arxiv.org/abs/1110.4375 1746cf0794eSJed Brown 175bcf0153eSBarry Smith Options Database Key: 17667b8a455SSatish Balay . -ts_arkimex_type prssp2 - set arkimex type to prssp2 1779bd3e852SBarry Smith 1786cf0794eSJed Brown Level: advanced 1796cf0794eSJed Brown 1801cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 1816cf0794eSJed Brown M*/ 182*d27334e2SStefano Zampini 1836cf0794eSJed Brown /*MC 18464f491ddSJed Brown TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part. 18564f491ddSJed Brown 18664f491ddSJed Brown This method has one explicit stage and three implicit stages. 18764f491ddSJed Brown 188bcf0153eSBarry Smith Options Database Key: 18967b8a455SSatish Balay . -ts_arkimex_type 3 - set arkimex type to 3 1909bd3e852SBarry Smith 191bcf0153eSBarry Smith Level: advanced 192bcf0153eSBarry Smith 19364f491ddSJed Brown References: 194606c0280SSatish Balay . * - Kennedy and Carpenter 2003. 19564f491ddSJed Brown 1961cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 19764f491ddSJed Brown M*/ 198*d27334e2SStefano Zampini 19964f491ddSJed Brown /*MC 2006cf0794eSJed Brown TSARKIMEXARS443 - Third order ARK IMEX scheme. 2016cf0794eSJed Brown 2026cf0794eSJed Brown This method has one explicit stage and four implicit stages. 2036cf0794eSJed Brown 204bcf0153eSBarry Smith Options Database Key: 20567b8a455SSatish Balay . -ts_arkimex_type ars443 - set arkimex type to ars443 2069bd3e852SBarry Smith 207bcf0153eSBarry Smith Level: advanced 208bcf0153eSBarry Smith 2096cf0794eSJed Brown References: 210606c0280SSatish Balay + * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 211606c0280SSatish Balay - * - This method is referred to as ARS(4,4,3) in https://arxiv.org/abs/1110.4375 2126cf0794eSJed Brown 2131cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 2146cf0794eSJed Brown M*/ 215*d27334e2SStefano Zampini 2166cf0794eSJed Brown /*MC 2176cf0794eSJed Brown TSARKIMEXBPR3 - Third order ARK IMEX scheme. 2186cf0794eSJed Brown 2196cf0794eSJed Brown This method has one explicit stage and four implicit stages. 2206cf0794eSJed Brown 221bcf0153eSBarry Smith Options Database Key: 22267b8a455SSatish Balay . -ts_arkimex_type bpr3 - set arkimex type to bpr3 2239bd3e852SBarry Smith 224bcf0153eSBarry Smith Level: advanced 225bcf0153eSBarry Smith 2266cf0794eSJed Brown References: 227606c0280SSatish Balay . * - This method is referred to as ARK3 in https://arxiv.org/abs/1110.4375 2286cf0794eSJed Brown 2291cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 2306cf0794eSJed Brown M*/ 231*d27334e2SStefano Zampini 2326cf0794eSJed Brown /*MC 23364f491ddSJed Brown TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part. 23464f491ddSJed Brown 23564f491ddSJed Brown This method has one explicit stage and four implicit stages. 23664f491ddSJed Brown 237bcf0153eSBarry Smith Options Database Key: 23867b8a455SSatish Balay . -ts_arkimex_type 4 - set arkimex type to4 2399bd3e852SBarry Smith 240bcf0153eSBarry Smith Level: advanced 241bcf0153eSBarry Smith 24264f491ddSJed Brown References: 243606c0280SSatish Balay . * - Kennedy and Carpenter 2003. 24464f491ddSJed Brown 2451cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 24664f491ddSJed Brown M*/ 247*d27334e2SStefano Zampini 24864f491ddSJed Brown /*MC 24964f491ddSJed Brown TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part. 25064f491ddSJed Brown 25164f491ddSJed Brown This method has one explicit stage and five implicit stages. 25264f491ddSJed Brown 253bcf0153eSBarry Smith Options Database Key: 25467b8a455SSatish Balay . -ts_arkimex_type 5 - set arkimex type to 5 2559bd3e852SBarry Smith 256bcf0153eSBarry Smith Level: advanced 257bcf0153eSBarry Smith 25864f491ddSJed Brown References: 259606c0280SSatish Balay . * - Kennedy and Carpenter 2003. 26064f491ddSJed Brown 2611cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 26264f491ddSJed Brown M*/ 26364f491ddSJed Brown 264*d27334e2SStefano Zampini static PetscErrorCode TSHasRHSFunction(TS ts, PetscBool *has) 265*d27334e2SStefano Zampini { 266*d27334e2SStefano Zampini TSRHSFunction func; 267*d27334e2SStefano Zampini 268*d27334e2SStefano Zampini PetscFunctionBegin; 269*d27334e2SStefano Zampini *has = PETSC_FALSE; 270*d27334e2SStefano Zampini PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL)); 271*d27334e2SStefano Zampini if (func) *has = PETSC_TRUE; 272*d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 273*d27334e2SStefano Zampini } 274*d27334e2SStefano Zampini 2758a381b04SJed Brown /*@C 276bcf0153eSBarry Smith TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX` 2778a381b04SJed Brown 278fca742c7SJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 2798a381b04SJed Brown 2808a381b04SJed Brown Level: advanced 2818a381b04SJed Brown 2821cc06b55SBarry Smith .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()` 2838a381b04SJed Brown @*/ 284d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXRegisterAll(void) 285d71ae5a4SJacob Faibussowitsch { 2868a381b04SJed Brown PetscFunctionBegin; 2873ba16761SJacob Faibussowitsch if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 2888a381b04SJed Brown TSARKIMEXRegisterAllCalled = PETSC_TRUE; 289e817cc15SEmil Constantinescu 290*d27334e2SStefano Zampini #define RC PetscRealConstant 291*d27334e2SStefano Zampini #define s2 RC(1.414213562373095048802) /* PetscSqrtReal((PetscReal)2.0) */ 292*d27334e2SStefano Zampini #define us2 RC(0.2928932188134524755992) /* 1.0-1.0/PetscSqrtReal((PetscReal)2.0); */ 293*d27334e2SStefano Zampini 294*d27334e2SStefano Zampini /* Diagonally implicit methods */ 295e817cc15SEmil Constantinescu { 296*d27334e2SStefano Zampini /* DIRK212, default of SUNDIALS */ 297*d27334e2SStefano Zampini const PetscReal A[2][2] = { 298*d27334e2SStefano Zampini {RC(1.0), RC(0.0)}, 299*d27334e2SStefano Zampini {RC(-1.0), RC(1.0)} 300*d27334e2SStefano Zampini }; 301*d27334e2SStefano Zampini const PetscReal b[2] = {RC(0.5), RC(0.5)}; 302*d27334e2SStefano Zampini const PetscReal bembed[2] = {RC(1.0), RC(0.0)}; 303*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKS212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 304*d27334e2SStefano Zampini } 305*d27334e2SStefano Zampini 3069371c9d4SSatish Balay { 307*d27334e2SStefano Zampini /* ESDIRK12 from https://arxiv.org/pdf/1803.01613.pdf */ 308*d27334e2SStefano Zampini const PetscReal A[2][2] = { 309*d27334e2SStefano Zampini {RC(0.0), RC(0.0)}, 310*d27334e2SStefano Zampini {RC(0.0), RC(1.0)} 311*d27334e2SStefano Zampini }; 312*d27334e2SStefano Zampini const PetscReal b[2] = {RC(0.0), RC(1.0)}; 313*d27334e2SStefano Zampini const PetscReal bembed[2] = {RC(0.5), RC(0.5)}; 314*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 315*d27334e2SStefano Zampini } 316*d27334e2SStefano Zampini 317*d27334e2SStefano Zampini { 318*d27334e2SStefano Zampini /* ESDIRK213L[2]SA from KC16. 319*d27334e2SStefano Zampini TR-BDF2 from Osea and Shampine */ 320*d27334e2SStefano Zampini const PetscReal A[3][3] = { 321*d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0)}, 322*d27334e2SStefano Zampini {us2, us2, RC(0.0)}, 323*d27334e2SStefano Zampini {s2 / RC(4.0), s2 / RC(4.0), us2 }, 324*d27334e2SStefano Zampini }; 325*d27334e2SStefano Zampini const PetscReal b[3] = {s2 / RC(4.0), s2 / RC(4.0), us2}; 326*d27334e2SStefano 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)}; 327*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES213SAL, 2, 3, &A[0][0], b, NULL, bembed, 1, b)); 328*d27334e2SStefano Zampini } 329*d27334e2SStefano Zampini 330*d27334e2SStefano Zampini { 331*d27334e2SStefano Zampini #define g RC(0.43586652150845899941601945) 332*d27334e2SStefano Zampini #define g2 PetscSqr(g) 333*d27334e2SStefano Zampini #define g3 g *g2 334*d27334e2SStefano Zampini #define g4 PetscSqr(g2) 335*d27334e2SStefano Zampini #define g5 g *g4 336*d27334e2SStefano Zampini #define c3 RC(1.0) 337*d27334e2SStefano Zampini #define a32 c3 *(c3 - RC(2.0) * g) / (RC(4.0) * g) 338*d27334e2SStefano 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)) 339*d27334e2SStefano Zampini #define b3 (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) / (RC(3.0) * c3 * (c3 - RC(2.0) * g)) 340*d27334e2SStefano Zampini #if 0 341*d27334e2SStefano Zampini /* This is for c3 = 3/5 */ 342*d27334e2SStefano Zampini #define bh2 \ 343*d27334e2SStefano 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)) 344*d27334e2SStefano 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)) 345*d27334e2SStefano 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) 346*d27334e2SStefano Zampini #else 347*d27334e2SStefano Zampini /* This is for c3 = 1.0 */ 348*d27334e2SStefano Zampini #define bh2 a32 349*d27334e2SStefano Zampini #define bh3 g 350*d27334e2SStefano Zampini #define bh4 RC(0.0) 351*d27334e2SStefano Zampini #endif 352*d27334e2SStefano Zampini /* ESDIRK3(2I)4L[2]SA from KC16 with c3 = 1.0 */ 353*d27334e2SStefano Zampini /* Given by Kvaerno https://link.springer.com/article/10.1023/b:bitn.0000046811.70614.38 */ 354*d27334e2SStefano Zampini const PetscReal A[4][4] = { 355*d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0), RC(0.0)}, 356*d27334e2SStefano Zampini {g, g, RC(0.0), RC(0.0)}, 357*d27334e2SStefano Zampini {c3 - a32 - g, a32, g, RC(0.0)}, 358*d27334e2SStefano Zampini {RC(1.0) - b2 - b3 - g, b2, b3, g }, 359*d27334e2SStefano Zampini }; 360*d27334e2SStefano Zampini const PetscReal b[4] = {RC(1.0) - b2 - b3 - g, b2, b3, g}; 361*d27334e2SStefano Zampini const PetscReal bembed[4] = {RC(1.0) - bh2 - bh3 - bh4, bh2, bh3, bh4}; 362*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES324SAL, 3, 4, &A[0][0], b, NULL, bembed, 1, b)); 363*d27334e2SStefano Zampini #undef g 364*d27334e2SStefano Zampini #undef g2 365*d27334e2SStefano Zampini #undef g3 366*d27334e2SStefano Zampini #undef c3 367*d27334e2SStefano Zampini #undef a32 368*d27334e2SStefano Zampini #undef b2 369*d27334e2SStefano Zampini #undef b3 370*d27334e2SStefano Zampini #undef bh2 371*d27334e2SStefano Zampini #undef bh3 372*d27334e2SStefano Zampini #undef bh4 373*d27334e2SStefano Zampini } 374*d27334e2SStefano Zampini 375*d27334e2SStefano Zampini { 376*d27334e2SStefano Zampini /* ESDIRK3(2I)5L[2]SA from KC16 */ 377*d27334e2SStefano Zampini const PetscReal A[5][5] = { 378*d27334e2SStefano Zampini {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 379*d27334e2SStefano Zampini {RC(9.0) / RC(40.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0), RC(0.0) }, 380*d27334e2SStefano 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) }, 381*d27334e2SStefano 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) }, 382*d27334e2SStefano 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)}, 383*d27334e2SStefano Zampini }; 384*d27334e2SStefano 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)}; 385*d27334e2SStefano 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)}; 386*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES325SAL, 3, 5, &A[0][0], b, NULL, bembed, 1, b)); 387*d27334e2SStefano Zampini } 388*d27334e2SStefano Zampini 389*d27334e2SStefano Zampini { 390*d27334e2SStefano Zampini // DIRK(6,6)[1]A[(7,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 391*d27334e2SStefano Zampini const PetscReal A[7][7] = { 392*d27334e2SStefano Zampini {RC(0.303487844706747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 393*d27334e2SStefano Zampini {RC(-0.279756492709814), RC(0.500032236020747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 394*d27334e2SStefano Zampini {RC(0.280583215743895), RC(-0.438560061586751), RC(0.217250734515736), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 395*d27334e2SStefano Zampini {RC(-0.0677678738539846), RC(0.984312781232293), RC(-0.266720192540149), RC(0.2476680834526), RC(0.0), RC(0.0), RC(0.0) }, 396*d27334e2SStefano Zampini {RC(0.125671616147993), RC(-0.995401751002415), RC(0.761333109549059), RC(-0.210281837202208), RC(0.866743712636936), RC(0.0), RC(0.0) }, 397*d27334e2SStefano Zampini {RC(-0.368056238801488), RC(-0.999928082701516), RC(0.534734253232519), RC(-0.174856916279082), RC(0.615007160285509), RC(0.696549912132029), RC(0.0) }, 398*d27334e2SStefano Zampini {RC(-0.00570546839653984), RC(-0.113110431835656), RC(-0.000965563207671587), RC(-0.000130490084629567), RC(0.00111737736895673), RC(-0.279385587378871), RC(0.618455906845342)} 399*d27334e2SStefano Zampini }; 400*d27334e2SStefano 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)}; 401*d27334e2SStefano 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)}; 402*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK657A, 6, 7, &A[0][0], b, NULL, bembed, 1, b)); 403*d27334e2SStefano Zampini } 404*d27334e2SStefano Zampini { 405*d27334e2SStefano Zampini // ESDIRK(8,6)[2]SA[(8,4)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 406*d27334e2SStefano Zampini const PetscReal A[8][8] = { 407*d27334e2SStefano 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) }, 408*d27334e2SStefano 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) }, 409*d27334e2SStefano 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) }, 410*d27334e2SStefano 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) }, 411*d27334e2SStefano 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) }, 412*d27334e2SStefano 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) }, 413*d27334e2SStefano 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) }, 414*d27334e2SStefano 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)} 415*d27334e2SStefano Zampini }; 416*d27334e2SStefano 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)}; 417*d27334e2SStefano 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)}; 418*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES648SA, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 419*d27334e2SStefano Zampini } 420*d27334e2SStefano Zampini { 421*d27334e2SStefano Zampini // DIRK(8,6)[1]SAL[(8,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 422*d27334e2SStefano Zampini const PetscReal A[8][8] = { 423*d27334e2SStefano 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) }, 424*d27334e2SStefano 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) }, 425*d27334e2SStefano 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) }, 426*d27334e2SStefano 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) }, 427*d27334e2SStefano 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) }, 428*d27334e2SStefano 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) }, 429*d27334e2SStefano 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) }, 430*d27334e2SStefano 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)} 431*d27334e2SStefano Zampini }; 432*d27334e2SStefano 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)}; 433*d27334e2SStefano 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)}; 434*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK658A, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 435*d27334e2SStefano Zampini } 436*d27334e2SStefano Zampini { 437*d27334e2SStefano Zampini // SDIRK(9,6)[1]SAL[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 438*d27334e2SStefano Zampini const PetscReal A[9][9] = { 439*d27334e2SStefano 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) }, 440*d27334e2SStefano 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) }, 441*d27334e2SStefano 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) }, 442*d27334e2SStefano 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) }, 443*d27334e2SStefano 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) }, 444*d27334e2SStefano 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) }, 445*d27334e2SStefano 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) }, 446*d27334e2SStefano 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) }, 447*d27334e2SStefano 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)} 448*d27334e2SStefano Zampini }; 449*d27334e2SStefano 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)}; 450*d27334e2SStefano 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)}; 451*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKS659A, 6, 9, &A[0][0], b, NULL, bembed, 1, b)); 452*d27334e2SStefano Zampini } 453*d27334e2SStefano Zampini { 454*d27334e2SStefano Zampini // DIRK(10,7)[1]SAL[(10,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 455*d27334e2SStefano Zampini const PetscReal A[10][10] = { 456*d27334e2SStefano 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) }, 457*d27334e2SStefano 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) }, 458*d27334e2SStefano 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) }, 459*d27334e2SStefano 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) }, 460*d27334e2SStefano 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) }, 461*d27334e2SStefano 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) }, 462*d27334e2SStefano 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) }, 463*d27334e2SStefano 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) }, 464*d27334e2SStefano 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) }, 465*d27334e2SStefano 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)} 466*d27334e2SStefano Zampini }; 467*d27334e2SStefano 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)}; 468*d27334e2SStefano Zampini const PetscReal bembed[10] = 469*d27334e2SStefano 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)}; 470*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK7510SAL, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 471*d27334e2SStefano Zampini } 472*d27334e2SStefano Zampini { 473*d27334e2SStefano Zampini // ESDIRK(10,7)[2]SA[(10,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 474*d27334e2SStefano Zampini const PetscReal A[10][10] = { 475*d27334e2SStefano 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) }, 476*d27334e2SStefano 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) }, 477*d27334e2SStefano 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) }, 478*d27334e2SStefano 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) }, 479*d27334e2SStefano 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) }, 480*d27334e2SStefano 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) }, 481*d27334e2SStefano 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) }, 482*d27334e2SStefano 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) }, 483*d27334e2SStefano 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) }, 484*d27334e2SStefano 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)} 485*d27334e2SStefano Zampini }; 486*d27334e2SStefano Zampini const PetscReal b[10] = {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), 487*d27334e2SStefano Zampini RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)}; 488*d27334e2SStefano Zampini const PetscReal bembed[10] = {RC(-0.015494246543626), RC(0.167657963820093), RC(0.269858958144236), RC(-0.0443258997755156), RC(0.150049236875266), 489*d27334e2SStefano Zampini RC(0.259452082755846), RC(0.244624573502521), RC(-0.215528446920284), RC(0.0487601760292619), RC(0.134945602112201)}; 490*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES7510SA, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 491*d27334e2SStefano Zampini } 492*d27334e2SStefano Zampini { 493*d27334e2SStefano Zampini // DIRK(9,7)[1]A[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 494*d27334e2SStefano Zampini const PetscReal A[9][9] = { 495*d27334e2SStefano 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) }, 496*d27334e2SStefano 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) }, 497*d27334e2SStefano 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) }, 498*d27334e2SStefano 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) }, 499*d27334e2SStefano 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) }, 500*d27334e2SStefano 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) }, 501*d27334e2SStefano 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) }, 502*d27334e2SStefano 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) }, 503*d27334e2SStefano 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)} 504*d27334e2SStefano Zampini }; 505*d27334e2SStefano Zampini 506*d27334e2SStefano 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)}; 507*d27334e2SStefano 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)}; 508*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK759A, 7, 9, &A[0][0], b, NULL, bembed, 1, b)); 509*d27334e2SStefano Zampini } 510*d27334e2SStefano Zampini { 511*d27334e2SStefano Zampini // SDIRK(11,7)[1]SAL[(11,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 512*d27334e2SStefano Zampini const PetscReal A[11][11] = { 513*d27334e2SStefano 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) }, 514*d27334e2SStefano 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) }, 515*d27334e2SStefano 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) }, 516*d27334e2SStefano 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) }, 517*d27334e2SStefano 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) }, 518*d27334e2SStefano 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) }, 519*d27334e2SStefano 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) }, 520*d27334e2SStefano 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) }, 521*d27334e2SStefano 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) }, 522*d27334e2SStefano 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) }, 523*d27334e2SStefano 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)} 524*d27334e2SStefano Zampini }; 525*d27334e2SStefano Zampini const PetscReal b[11] = 526*d27334e2SStefano 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)}; 527*d27334e2SStefano Zampini const PetscReal bembed[11] = 528*d27334e2SStefano 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)}; 529*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKS7511SAL, 7, 11, &A[0][0], b, NULL, bembed, 1, b)); 530*d27334e2SStefano Zampini } 531*d27334e2SStefano Zampini { 532*d27334e2SStefano Zampini // DIRK(13,8)[1]A[(14,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 533*d27334e2SStefano Zampini const PetscReal A[14][14] = { 534*d27334e2SStefano 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) }, 535*d27334e2SStefano 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) }, 536*d27334e2SStefano 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) }, 537*d27334e2SStefano 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) }, 538*d27334e2SStefano 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) }, 539*d27334e2SStefano 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) }, 540*d27334e2SStefano 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) }, 541*d27334e2SStefano 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) }, 542*d27334e2SStefano 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) }, 543*d27334e2SStefano 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) }, 544*d27334e2SStefano 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) }, 545*d27334e2SStefano 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) }, 546*d27334e2SStefano 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) }, 547*d27334e2SStefano 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)} 548*d27334e2SStefano Zampini }; 549*d27334e2SStefano 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)}; 550*d27334e2SStefano 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), 551*d27334e2SStefano Zampini RC(0.0417664613347638), RC(0.223636394275293), RC(0.231569156867596), RC(0.240526201277663), RC(-0.222933582911926), RC(-0.0111479879597561), RC(0.19915314335888)}; 552*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK8614A, 8, 14, &A[0][0], b, NULL, bembed, 1, b)); 553*d27334e2SStefano Zampini } 554*d27334e2SStefano Zampini { 555*d27334e2SStefano Zampini // DIRK(15,8)[1]SAL[(16,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 556*d27334e2SStefano Zampini const PetscReal A[16][16] = { 557*d27334e2SStefano 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) }, 558*d27334e2SStefano 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) }, 559*d27334e2SStefano 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) }, 560*d27334e2SStefano 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) }, 561*d27334e2SStefano 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) }, 562*d27334e2SStefano 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) }, 563*d27334e2SStefano 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) }, 564*d27334e2SStefano 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) }, 565*d27334e2SStefano 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) }, 566*d27334e2SStefano 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) }, 567*d27334e2SStefano 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) }, 568*d27334e2SStefano 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) }, 569*d27334e2SStefano 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) }, 570*d27334e2SStefano 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) }, 571*d27334e2SStefano 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) }, 572*d27334e2SStefano 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)} 573*d27334e2SStefano Zampini }; 574*d27334e2SStefano 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)}; 575*d27334e2SStefano 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), 576*d27334e2SStefano 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)}; 577*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRK8616SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 578*d27334e2SStefano Zampini } 579*d27334e2SStefano Zampini { 580*d27334e2SStefano Zampini // ESDIRK(16,8)[2]SAL[(16,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 581*d27334e2SStefano Zampini const PetscReal A[16][16] = { 582*d27334e2SStefano 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) }, 583*d27334e2SStefano 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) }, 584*d27334e2SStefano 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) }, 585*d27334e2SStefano 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) }, 586*d27334e2SStefano 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) }, 587*d27334e2SStefano 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) }, 588*d27334e2SStefano 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) }, 589*d27334e2SStefano 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) }, 590*d27334e2SStefano 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) }, 591*d27334e2SStefano 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) }, 592*d27334e2SStefano 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) }, 593*d27334e2SStefano 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) }, 594*d27334e2SStefano 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) }, 595*d27334e2SStefano 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) }, 596*d27334e2SStefano 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) }, 597*d27334e2SStefano 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)} 598*d27334e2SStefano Zampini }; 599*d27334e2SStefano 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), 600*d27334e2SStefano 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)}; 601*d27334e2SStefano 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), 602*d27334e2SStefano 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)}; 603*d27334e2SStefano Zampini PetscCall(TSDIRKRegister(TSDIRKES8516SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 604*d27334e2SStefano Zampini } 605*d27334e2SStefano Zampini 606*d27334e2SStefano Zampini /* Additive methods */ 607*d27334e2SStefano Zampini { 608*d27334e2SStefano Zampini const PetscReal A[3][3] = { 609e817cc15SEmil Constantinescu {0.0, 0.0, 0.0}, 6109371c9d4SSatish Balay {0.0, 0.0, 0.0}, 6119371c9d4SSatish Balay {0.0, 0.5, 0.0} 612*d27334e2SStefano Zampini }; 613*d27334e2SStefano Zampini const PetscReal At[3][3] = { 614*d27334e2SStefano Zampini {1.0, 0.0, 0.0}, 615*d27334e2SStefano Zampini {0.0, 0.5, 0.0}, 616*d27334e2SStefano Zampini {0.0, 0.5, 0.5} 617*d27334e2SStefano Zampini }; 618*d27334e2SStefano Zampini const PetscReal b[3] = {0.0, 0.5, 0.5}; 619*d27334e2SStefano Zampini const PetscReal bembedt[3] = {1.0, 0.0, 0.0}; 6209566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 621e817cc15SEmil Constantinescu } 6228a381b04SJed Brown { 623*d27334e2SStefano Zampini const PetscReal A[2][2] = { 6249371c9d4SSatish Balay {0.0, 0.0}, 6259371c9d4SSatish Balay {0.5, 0.0} 626*d27334e2SStefano Zampini }; 627*d27334e2SStefano Zampini const PetscReal At[2][2] = { 628*d27334e2SStefano Zampini {0.0, 0.0}, 629*d27334e2SStefano Zampini {0.0, 0.5} 630*d27334e2SStefano Zampini }; 631*d27334e2SStefano Zampini const PetscReal b[2] = {0.0, 1.0}; 632*d27334e2SStefano Zampini const PetscReal bembedt[2] = {0.5, 0.5}; 6331f80e275SEmil 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 */ 6349566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 6351f80e275SEmil Constantinescu } 6361f80e275SEmil Constantinescu { 637*d27334e2SStefano Zampini const PetscReal A[2][2] = { 6389371c9d4SSatish Balay {0.0, 0.0}, 6399371c9d4SSatish Balay {1.0, 0.0} 640*d27334e2SStefano Zampini }; 641*d27334e2SStefano Zampini const PetscReal At[2][2] = { 642*d27334e2SStefano Zampini {0.0, 0.0}, 643*d27334e2SStefano Zampini {0.5, 0.5} 644*d27334e2SStefano Zampini }; 645*d27334e2SStefano Zampini const PetscReal b[2] = {0.5, 0.5}; 646*d27334e2SStefano Zampini const PetscReal bembedt[2] = {0.0, 1.0}; 6471f80e275SEmil 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 */ 6489566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 6491f80e275SEmil Constantinescu } 6501f80e275SEmil Constantinescu { 651*d27334e2SStefano Zampini const PetscReal A[2][2] = { 6529371c9d4SSatish Balay {0.0, 0.0}, 6539371c9d4SSatish Balay {1.0, 0.0} 654*d27334e2SStefano Zampini }; 655*d27334e2SStefano Zampini const PetscReal At[2][2] = { 656*d27334e2SStefano Zampini {us2, 0.0}, 657*d27334e2SStefano Zampini {1.0 - 2.0 * us2, us2} 658*d27334e2SStefano Zampini }; 659*d27334e2SStefano Zampini const PetscReal b[2] = {0.5, 0.5}; 660*d27334e2SStefano Zampini const PetscReal bembedt[2] = {0.0, 1.0}; 661*d27334e2SStefano Zampini const PetscReal binterpt[2][2] = { 662*d27334e2SStefano Zampini {(us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))}, 663*d27334e2SStefano Zampini {1 - (us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))} 664*d27334e2SStefano Zampini }; 665*d27334e2SStefano Zampini const PetscReal binterp[2][2] = { 666*d27334e2SStefano Zampini {1.0, -0.5}, 667*d27334e2SStefano Zampini {0.0, 0.5 } 668*d27334e2SStefano Zampini }; 6699566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 6701f80e275SEmil Constantinescu } 6711f80e275SEmil Constantinescu { 672*d27334e2SStefano Zampini const PetscReal A[3][3] = { 6739371c9d4SSatish Balay {0, 0, 0}, 674*d27334e2SStefano Zampini {2 - s2, 0, 0}, 6759371c9d4SSatish Balay {0.5, 0.5, 0} 676*d27334e2SStefano Zampini }; 677*d27334e2SStefano Zampini const PetscReal At[3][3] = { 678*d27334e2SStefano Zampini {0, 0, 0 }, 679*d27334e2SStefano Zampini {1 - 1 / s2, 1 - 1 / s2, 0 }, 680*d27334e2SStefano Zampini {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 681*d27334e2SStefano Zampini }; 682*d27334e2SStefano Zampini const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 683*d27334e2SStefano Zampini const PetscReal binterpt[3][2] = { 684*d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 685*d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 686*d27334e2SStefano Zampini {1.0 - s2, 1.0 / s2 } 687*d27334e2SStefano Zampini }; 6889566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 6891f80e275SEmil Constantinescu } 6901f80e275SEmil Constantinescu { 691*d27334e2SStefano Zampini const PetscReal A[3][3] = { 6929371c9d4SSatish Balay {0, 0, 0}, 693*d27334e2SStefano Zampini {2 - s2, 0, 0}, 6949371c9d4SSatish Balay {0.75, 0.25, 0} 695*d27334e2SStefano Zampini }; 696*d27334e2SStefano Zampini const PetscReal At[3][3] = { 697*d27334e2SStefano Zampini {0, 0, 0 }, 698*d27334e2SStefano Zampini {1 - 1 / s2, 1 - 1 / s2, 0 }, 699*d27334e2SStefano Zampini {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 700*d27334e2SStefano Zampini }; 701*d27334e2SStefano Zampini const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 702*d27334e2SStefano Zampini const PetscReal binterpt[3][2] = { 703*d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 704*d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 705*d27334e2SStefano Zampini {1.0 - s2, 1.0 / s2 } 706*d27334e2SStefano Zampini }; 7079566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 7088a381b04SJed Brown } 70906db7b1cSJed Brown { /* Optimal for linear implicit part */ 710*d27334e2SStefano Zampini const PetscReal A[3][3] = { 7119371c9d4SSatish Balay {0, 0, 0}, 712*d27334e2SStefano Zampini {2 - s2, 0, 0}, 713*d27334e2SStefano Zampini {(3 - 2 * s2) / 6, (3 + 2 * s2) / 6, 0} 714*d27334e2SStefano Zampini }; 715*d27334e2SStefano Zampini const PetscReal At[3][3] = { 716*d27334e2SStefano Zampini {0, 0, 0 }, 717*d27334e2SStefano Zampini {1 - 1 / s2, 1 - 1 / s2, 0 }, 718*d27334e2SStefano Zampini {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 719*d27334e2SStefano Zampini }; 720*d27334e2SStefano Zampini const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 721*d27334e2SStefano Zampini const PetscReal binterpt[3][2] = { 722*d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 723*d27334e2SStefano Zampini {1.0 / s2, -1.0 / (2.0 * s2)}, 724*d27334e2SStefano Zampini {1.0 - s2, 1.0 / s2 } 725*d27334e2SStefano Zampini }; 7269566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 727a3a57f36SJed Brown } 7286cf0794eSJed Brown { /* Optimal for linear implicit part */ 729*d27334e2SStefano Zampini const PetscReal A[3][3] = { 7309371c9d4SSatish Balay {0, 0, 0}, 7316cf0794eSJed Brown {0.5, 0, 0}, 7329371c9d4SSatish Balay {0.5, 0.5, 0} 733*d27334e2SStefano Zampini }; 734*d27334e2SStefano Zampini const PetscReal At[3][3] = { 735*d27334e2SStefano Zampini {0.25, 0, 0 }, 736*d27334e2SStefano Zampini {0, 0.25, 0 }, 737*d27334e2SStefano Zampini {1. / 3, 1. / 3, 1. / 3} 738*d27334e2SStefano Zampini }; 7399566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 7406cf0794eSJed Brown } 741a3a57f36SJed Brown { 742*d27334e2SStefano Zampini const PetscReal A[4][4] = { 7439371c9d4SSatish Balay {0, 0, 0, 0}, 7444040e9f2SJed Brown {1767732205903. / 2027836641118., 0, 0, 0}, 7454040e9f2SJed Brown {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 7469371c9d4SSatish Balay {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 747*d27334e2SStefano Zampini }; 748*d27334e2SStefano Zampini const PetscReal At[4][4] = { 749*d27334e2SStefano Zampini {0, 0, 0, 0 }, 750*d27334e2SStefano Zampini {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0 }, 751*d27334e2SStefano Zampini {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0 }, 752*d27334e2SStefano Zampini {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.} 753*d27334e2SStefano Zampini }; 754*d27334e2SStefano Zampini const PetscReal bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}; 755*d27334e2SStefano Zampini const PetscReal binterpt[4][2] = { 756*d27334e2SStefano Zampini {4655552711362. / 22874653954995., -215264564351. / 13552729205753. }, 757*d27334e2SStefano Zampini {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119. }, 758*d27334e2SStefano Zampini {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, 759*d27334e2SStefano Zampini {584795268549. / 6622622206610., 2508943948391. / 7218656332882. } 760*d27334e2SStefano Zampini }; 7619566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 762a3a57f36SJed Brown } 763a3a57f36SJed Brown { 764*d27334e2SStefano Zampini const PetscReal A[5][5] = { 7659371c9d4SSatish Balay {0, 0, 0, 0, 0}, 7666cf0794eSJed Brown {1. / 2, 0, 0, 0, 0}, 7676cf0794eSJed Brown {11. / 18, 1. / 18, 0, 0, 0}, 7686cf0794eSJed Brown {5. / 6, -5. / 6, .5, 0, 0}, 7699371c9d4SSatish Balay {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 770*d27334e2SStefano Zampini }; 771*d27334e2SStefano Zampini const PetscReal At[5][5] = { 772*d27334e2SStefano Zampini {0, 0, 0, 0, 0 }, 773*d27334e2SStefano Zampini {0, 1. / 2, 0, 0, 0 }, 774*d27334e2SStefano Zampini {0, 1. / 6, 1. / 2, 0, 0 }, 775*d27334e2SStefano Zampini {0, -1. / 2, 1. / 2, 1. / 2, 0 }, 776*d27334e2SStefano Zampini {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2} 777*d27334e2SStefano Zampini }; 778*d27334e2SStefano Zampini PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 7796cf0794eSJed Brown } 7806cf0794eSJed Brown { 781*d27334e2SStefano Zampini const PetscReal A[5][5] = { 7829371c9d4SSatish Balay {0, 0, 0, 0, 0}, 7836cf0794eSJed Brown {1, 0, 0, 0, 0}, 7846cf0794eSJed Brown {4. / 9, 2. / 9, 0, 0, 0}, 7856cf0794eSJed Brown {1. / 4, 0, 3. / 4, 0, 0}, 7869371c9d4SSatish Balay {1. / 4, 0, 3. / 5, 0, 0} 787*d27334e2SStefano Zampini }; 788*d27334e2SStefano Zampini const PetscReal At[5][5] = { 789*d27334e2SStefano Zampini {0, 0, 0, 0, 0 }, 790*d27334e2SStefano Zampini {.5, .5, 0, 0, 0 }, 791*d27334e2SStefano Zampini {5. / 18, -1. / 9, .5, 0, 0 }, 792*d27334e2SStefano Zampini {.5, 0, 0, .5, 0 }, 793*d27334e2SStefano Zampini {.25, 0, .75, -.5, .5} 794*d27334e2SStefano Zampini }; 795*d27334e2SStefano Zampini PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 7966cf0794eSJed Brown } 7976cf0794eSJed Brown { 798*d27334e2SStefano Zampini const PetscReal A[6][6] = { 7999371c9d4SSatish Balay {0, 0, 0, 0, 0, 0}, 800a3a57f36SJed Brown {1. / 2, 0, 0, 0, 0, 0}, 8014040e9f2SJed Brown {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 8024040e9f2SJed Brown {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 8034040e9f2SJed Brown {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 8049371c9d4SSatish Balay {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 805*d27334e2SStefano Zampini }; 806*d27334e2SStefano Zampini const PetscReal At[6][6] = { 807*d27334e2SStefano Zampini {0, 0, 0, 0, 0, 0 }, 808*d27334e2SStefano Zampini {1. / 4, 1. / 4, 0, 0, 0, 0 }, 809*d27334e2SStefano Zampini {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0 }, 810*d27334e2SStefano Zampini {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0 }, 811*d27334e2SStefano Zampini {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0 }, 812*d27334e2SStefano Zampini {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4} 813*d27334e2SStefano Zampini }; 814*d27334e2SStefano Zampini const PetscReal bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}; 815*d27334e2SStefano Zampini const PetscReal binterpt[6][3] = { 816*d27334e2SStefano Zampini {6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025. }, 817*d27334e2SStefano Zampini {0, 0, 0 }, 818*d27334e2SStefano Zampini {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035. }, 819*d27334e2SStefano Zampini {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 820*d27334e2SStefano Zampini {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469. }, 821*d27334e2SStefano Zampini {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.} 822*d27334e2SStefano Zampini }; 8239566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 824a3a57f36SJed Brown } 825a3a57f36SJed Brown { 826*d27334e2SStefano Zampini const PetscReal A[8][8] = { 8279371c9d4SSatish Balay {0, 0, 0, 0, 0, 0, 0, 0}, 828a3a57f36SJed Brown {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 8294040e9f2SJed Brown {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 8304040e9f2SJed Brown {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 8314040e9f2SJed Brown {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 8324040e9f2SJed Brown {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 8334040e9f2SJed Brown {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 8349371c9d4SSatish Balay {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 835*d27334e2SStefano Zampini }; 836*d27334e2SStefano Zampini const PetscReal At[8][8] = { 837*d27334e2SStefano Zampini {0, 0, 0, 0, 0, 0, 0, 0 }, 838*d27334e2SStefano Zampini {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0 }, 839*d27334e2SStefano Zampini {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0 }, 840*d27334e2SStefano Zampini {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0 }, 841*d27334e2SStefano Zampini {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0 }, 842*d27334e2SStefano Zampini {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0 }, 843*d27334e2SStefano Zampini {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0 }, 844*d27334e2SStefano Zampini {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.} 845*d27334e2SStefano Zampini }; 846*d27334e2SStefano Zampini const PetscReal bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}; 847*d27334e2SStefano Zampini const PetscReal binterpt[8][3] = { 848*d27334e2SStefano Zampini {-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439. }, 849*d27334e2SStefano Zampini {0, 0, 0 }, 850*d27334e2SStefano Zampini {0, 0, 0 }, 851*d27334e2SStefano Zampini {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, 852*d27334e2SStefano Zampini {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, 853*d27334e2SStefano Zampini {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, 854*d27334e2SStefano Zampini {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, 855*d27334e2SStefano Zampini {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460. } 856*d27334e2SStefano Zampini }; 8579566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 858a3a57f36SJed Brown } 859*d27334e2SStefano Zampini #undef RC 860*d27334e2SStefano Zampini #undef us2 861*d27334e2SStefano Zampini #undef s2 8623ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 8638a381b04SJed Brown } 8648a381b04SJed Brown 8658a381b04SJed Brown /*@C 866bcf0153eSBarry Smith TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`. 8678a381b04SJed Brown 8688a381b04SJed Brown Not Collective 8698a381b04SJed Brown 8708a381b04SJed Brown Level: advanced 8718a381b04SJed Brown 8721cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 8738a381b04SJed Brown @*/ 874d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXRegisterDestroy(void) 875d71ae5a4SJacob Faibussowitsch { 8768a381b04SJed Brown ARKTableauLink link; 8778a381b04SJed Brown 8788a381b04SJed Brown PetscFunctionBegin; 8798a381b04SJed Brown while ((link = ARKTableauList)) { 8808a381b04SJed Brown ARKTableau t = &link->tab; 8818a381b04SJed Brown ARKTableauList = link->next; 8829566063dSJacob Faibussowitsch PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 8839566063dSJacob Faibussowitsch PetscCall(PetscFree2(t->bembedt, t->bembed)); 8849566063dSJacob Faibussowitsch PetscCall(PetscFree2(t->binterpt, t->binterp)); 8859566063dSJacob Faibussowitsch PetscCall(PetscFree(t->name)); 8869566063dSJacob Faibussowitsch PetscCall(PetscFree(link)); 8878a381b04SJed Brown } 8888a381b04SJed Brown TSARKIMEXRegisterAllCalled = PETSC_FALSE; 8893ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 8908a381b04SJed Brown } 8918a381b04SJed Brown 8928a381b04SJed Brown /*@C 893bcf0153eSBarry Smith TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called 894bcf0153eSBarry Smith from `TSInitializePackage()`. 8958a381b04SJed Brown 8968a381b04SJed Brown Level: developer 8978a381b04SJed Brown 8981cc06b55SBarry Smith .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()` 8998a381b04SJed Brown @*/ 900d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXInitializePackage(void) 901d71ae5a4SJacob Faibussowitsch { 9028a381b04SJed Brown PetscFunctionBegin; 9033ba16761SJacob Faibussowitsch if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 9048a381b04SJed Brown TSARKIMEXPackageInitialized = PETSC_TRUE; 9059566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegisterAll()); 9069566063dSJacob Faibussowitsch PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 9073ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 9088a381b04SJed Brown } 9098a381b04SJed Brown 9108a381b04SJed Brown /*@C 911bcf0153eSBarry Smith TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is 912bcf0153eSBarry Smith called from `PetscFinalize()`. 9138a381b04SJed Brown 9148a381b04SJed Brown Level: developer 9158a381b04SJed Brown 9161cc06b55SBarry Smith .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()` 9178a381b04SJed Brown @*/ 918d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXFinalizePackage(void) 919d71ae5a4SJacob Faibussowitsch { 9208a381b04SJed Brown PetscFunctionBegin; 9218a381b04SJed Brown TSARKIMEXPackageInitialized = PETSC_FALSE; 9229566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRegisterDestroy()); 9233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 9248a381b04SJed Brown } 9258a381b04SJed Brown 926cd652676SJed Brown /*@C 927bcf0153eSBarry Smith TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 928cd652676SJed Brown 929*d27334e2SStefano Zampini Logically Collective. 930cd652676SJed Brown 931cd652676SJed Brown Input Parameters: 932cd652676SJed Brown + name - identifier for method 933cd652676SJed Brown . order - approximation order of method 934cd652676SJed Brown . s - number of stages, this is the dimension of the matrices below 935cd652676SJed Brown . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 9360298fd71SBarry Smith . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 9370298fd71SBarry Smith . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 938cd652676SJed Brown . A - Non-stiff stage coefficients (dimension s*s, row-major) 9390298fd71SBarry Smith . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 9400298fd71SBarry Smith . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 9410298fd71SBarry Smith . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 9420298fd71SBarry Smith . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 943cd652676SJed Brown . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 944cd652676SJed Brown . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 9450298fd71SBarry Smith - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 946cd652676SJed Brown 947cd652676SJed Brown Level: advanced 948cd652676SJed Brown 949bcf0153eSBarry Smith Note: 950bcf0153eSBarry Smith Several `TSARKIMEX` methods are provided, this function is only needed to create new methods. 951bcf0153eSBarry Smith 9521cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS` 953cd652676SJed Brown @*/ 954d71ae5a4SJacob 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[]) 955d71ae5a4SJacob Faibussowitsch { 9568a381b04SJed Brown ARKTableauLink link; 9578a381b04SJed Brown ARKTableau t; 9588a381b04SJed Brown PetscInt i, j; 9598a381b04SJed Brown 9608a381b04SJed Brown PetscFunctionBegin; 9619566063dSJacob Faibussowitsch PetscCall(TSARKIMEXInitializePackage()); 962*d27334e2SStefano Zampini for (link = ARKTableauList; link; link = link->next) { 963*d27334e2SStefano Zampini PetscBool match; 964*d27334e2SStefano Zampini 965*d27334e2SStefano Zampini PetscCall(PetscStrcmp(link->tab.name, name, &match)); 966*d27334e2SStefano Zampini PetscCheck(!match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Method with name \"%s\" already registered", name); 967*d27334e2SStefano Zampini } 9689566063dSJacob Faibussowitsch PetscCall(PetscNew(&link)); 9698a381b04SJed Brown t = &link->tab; 9709566063dSJacob Faibussowitsch PetscCall(PetscStrallocpy(name, &t->name)); 9718a381b04SJed Brown t->order = order; 9728a381b04SJed Brown t->s = s; 9739566063dSJacob Faibussowitsch PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 9749566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->At, At, s * s)); 975*d27334e2SStefano Zampini if (A) { 9769566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->A, A, s * s)); 977*d27334e2SStefano Zampini t->additive = PETSC_TRUE; 978*d27334e2SStefano Zampini } 979*d27334e2SStefano Zampini 9809566063dSJacob Faibussowitsch if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 9819371c9d4SSatish Balay else 9829371c9d4SSatish Balay for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 983*d27334e2SStefano Zampini 984*d27334e2SStefano Zampini if (t->additive) { 9859566063dSJacob Faibussowitsch if (b) PetscCall(PetscArraycpy(t->b, b, s)); 9869371c9d4SSatish Balay else 9879371c9d4SSatish Balay for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 988*d27334e2SStefano Zampini } else PetscCall(PetscArrayzero(t->b, s)); 989*d27334e2SStefano Zampini 9909566063dSJacob Faibussowitsch if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 9919371c9d4SSatish Balay else 9929371c9d4SSatish Balay for (i = 0; i < s; i++) 9939371c9d4SSatish Balay for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 994*d27334e2SStefano Zampini 995*d27334e2SStefano Zampini if (t->additive) { 9969566063dSJacob Faibussowitsch if (c) PetscCall(PetscArraycpy(t->c, c, s)); 9979371c9d4SSatish Balay else 9989371c9d4SSatish Balay for (i = 0; i < s; i++) 9999371c9d4SSatish Balay for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 1000*d27334e2SStefano Zampini } else PetscCall(PetscArrayzero(t->c, s)); 1001*d27334e2SStefano Zampini 1002e817cc15SEmil Constantinescu t->stiffly_accurate = PETSC_TRUE; 10039371c9d4SSatish Balay for (i = 0; i < s; i++) 10049371c9d4SSatish Balay if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 1005*d27334e2SStefano Zampini 1006e817cc15SEmil Constantinescu t->explicit_first_stage = PETSC_TRUE; 10079371c9d4SSatish Balay for (i = 0; i < s; i++) 10089371c9d4SSatish Balay if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 1009*d27334e2SStefano Zampini 1010e817cc15SEmil Constantinescu /* def of FSAL can be made more precise */ 10114e9d4bf5SJed Brown t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 1012*d27334e2SStefano Zampini 1013108c343cSJed Brown if (bembedt) { 10149566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 10159566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 10169566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 1017108c343cSJed Brown } 1018108c343cSJed Brown 10194f385281SJed Brown t->pinterp = pinterp; 10209566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 10219566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 10229566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 1023*d27334e2SStefano Zampini 10248a381b04SJed Brown link->next = ARKTableauList; 10258a381b04SJed Brown ARKTableauList = link; 10263ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10278a381b04SJed Brown } 10288a381b04SJed Brown 1029*d27334e2SStefano Zampini /*@C 1030*d27334e2SStefano Zampini TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation 1031*d27334e2SStefano Zampini 1032*d27334e2SStefano Zampini Logically Collective. 1033*d27334e2SStefano Zampini 1034*d27334e2SStefano Zampini Input Parameters: 1035*d27334e2SStefano Zampini + name - identifier for method 1036*d27334e2SStefano Zampini . order - approximation order of method 1037*d27334e2SStefano Zampini . s - number of stages, this is the dimension of the matrices below 1038*d27334e2SStefano Zampini . At - Butcher table of stage coefficients (dimension `s`*`s`, row-major order) 1039*d27334e2SStefano Zampini . bt - Butcher table for completing the step (dimension `s`; pass `NULL` to use the last row of `At`) 1040*d27334e2SStefano Zampini . ct - Abscissa of each stage (dimension s, NULL to use row sums of At) 1041*d27334e2SStefano Zampini . bembedt - Stiff part of completion table for embedded method (dimension s; `NULL` if not available) 1042*d27334e2SStefano Zampini . pinterp - Order of the interpolation scheme, equal to the number of columns of `binterpt` and `binterp` 1043*d27334e2SStefano Zampini - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 1044*d27334e2SStefano Zampini 1045*d27334e2SStefano Zampini Level: advanced 1046*d27334e2SStefano Zampini 1047*d27334e2SStefano Zampini Note: 1048*d27334e2SStefano Zampini Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods. 1049*d27334e2SStefano Zampini 1050*d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS` 1051*d27334e2SStefano Zampini @*/ 1052*d27334e2SStefano 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[]) 1053*d27334e2SStefano Zampini { 1054*d27334e2SStefano Zampini PetscFunctionBegin; 1055*d27334e2SStefano Zampini PetscCall(TSARKIMEXRegister(name, order, s, At, bt, ct, NULL, NULL, NULL, bembedt, NULL, pinterp, binterpt, NULL)); 1056*d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 1057*d27334e2SStefano Zampini } 1058*d27334e2SStefano Zampini 1059108c343cSJed Brown /* 1060108c343cSJed Brown The step completion formula is 1061108c343cSJed Brown 1062108c343cSJed Brown x1 = x0 - h bt^T YdotI + h b^T YdotRHS 1063108c343cSJed Brown 1064108c343cSJed Brown This function can be called before or after ts->vec_sol has been updated. 1065108c343cSJed Brown Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 1066108c343cSJed Brown We can write 1067108c343cSJed Brown 1068108c343cSJed Brown x1e = x0 - h bet^T YdotI + h be^T YdotRHS 1069108c343cSJed Brown = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 1070108c343cSJed Brown = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 1071108c343cSJed Brown 1072108c343cSJed Brown so we can evaluate the method with different order even after the step has been optimistically completed. 1073108c343cSJed Brown */ 1074d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) 1075d71ae5a4SJacob Faibussowitsch { 1076108c343cSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1077108c343cSJed Brown ARKTableau tab = ark->tableau; 1078108c343cSJed Brown PetscScalar *w = ark->work; 1079108c343cSJed Brown PetscReal h; 1080108c343cSJed Brown PetscInt s = tab->s, j; 1081*d27334e2SStefano Zampini PetscBool hasE; 1082108c343cSJed Brown 1083108c343cSJed Brown PetscFunctionBegin; 1084108c343cSJed Brown switch (ark->status) { 1085108c343cSJed Brown case TS_STEP_INCOMPLETE: 1086d71ae5a4SJacob Faibussowitsch case TS_STEP_PENDING: 1087d71ae5a4SJacob Faibussowitsch h = ts->time_step; 1088d71ae5a4SJacob Faibussowitsch break; 1089d71ae5a4SJacob Faibussowitsch case TS_STEP_COMPLETE: 1090d71ae5a4SJacob Faibussowitsch h = ts->ptime - ts->ptime_prev; 1091d71ae5a4SJacob Faibussowitsch break; 1092d71ae5a4SJacob Faibussowitsch default: 1093d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1094108c343cSJed Brown } 1095108c343cSJed Brown if (order == tab->order) { 1096e817cc15SEmil Constantinescu if (ark->status == TS_STEP_INCOMPLETE) { 1097740132f1SEmil Constantinescu if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 10989566063dSJacob Faibussowitsch PetscCall(VecCopy(ark->Y[s - 1], X)); 1099e817cc15SEmil Constantinescu } else { /* Use the standard completion formula (bt,b) */ 11009566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, X)); 1101e817cc15SEmil Constantinescu for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 11029566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1103*d27334e2SStefano Zampini if (tab->additive && ark->imex) { /* Method is IMEX, complete the explicit formula */ 1104*d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1105*d27334e2SStefano Zampini if (hasE) { 1106108c343cSJed Brown for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 11079566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1108e817cc15SEmil Constantinescu } 1109e817cc15SEmil Constantinescu } 1110*d27334e2SStefano Zampini } 11119566063dSJacob Faibussowitsch } else PetscCall(VecCopy(ts->vec_sol, X)); 1112108c343cSJed Brown if (done) *done = PETSC_TRUE; 11133ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1114108c343cSJed Brown } else if (order == tab->order - 1) { 1115108c343cSJed Brown if (!tab->bembedt) goto unavailable; 1116108c343cSJed Brown if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 11179566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, X)); 1118e817cc15SEmil Constantinescu for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 11199566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1120*d27334e2SStefano Zampini if (tab->additive) { 1121*d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1122*d27334e2SStefano Zampini if (hasE) { 1123108c343cSJed Brown for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 11249566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1125*d27334e2SStefano Zampini } 1126*d27334e2SStefano Zampini } 1127108c343cSJed Brown } else { /* Rollback and re-complete using (bet-be,be-b) */ 11289566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, X)); 1129e817cc15SEmil Constantinescu for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 11309566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 1131*d27334e2SStefano Zampini if (tab->additive) { 1132*d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1133*d27334e2SStefano Zampini if (hasE) { 1134108c343cSJed Brown for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 11359566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1136108c343cSJed Brown } 1137*d27334e2SStefano Zampini } 1138*d27334e2SStefano Zampini } 1139108c343cSJed Brown if (done) *done = PETSC_TRUE; 11403ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1141108c343cSJed Brown } 1142108c343cSJed Brown unavailable: 1143108c343cSJed Brown if (done) *done = PETSC_FALSE; 11449371c9d4SSatish Balay else 11459371c9d4SSatish 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, 11469371c9d4SSatish Balay tab->order, order); 11473ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1148108c343cSJed Brown } 1149108c343cSJed Brown 1150d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 1151d71ae5a4SJacob Faibussowitsch { 1152c79dcfbdSBarry Smith Vec Udot, Y1, Y2; 1153c79dcfbdSBarry Smith TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1154c79dcfbdSBarry Smith PetscReal norm; 1155c79dcfbdSBarry Smith 1156c79dcfbdSBarry Smith PetscFunctionBegin; 11579566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 11589566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 11599566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 11609566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 11619566063dSJacob Faibussowitsch PetscCall(VecSetRandom(Udot, NULL)); 11629566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 11639566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y2, -1.0, Y1)); 11649566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y2, -1.0, Udot)); 11659566063dSJacob Faibussowitsch PetscCall(VecNorm(Y2, NORM_2, &norm)); 1166c79dcfbdSBarry Smith if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 1167c79dcfbdSBarry Smith *id = PETSC_TRUE; 1168c79dcfbdSBarry Smith } else { 1169c79dcfbdSBarry Smith *id = PETSC_FALSE; 11709566063dSJacob 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)); 1171c79dcfbdSBarry Smith } 11729566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Udot)); 11739566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Y1)); 11749566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Y2)); 11753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1176c79dcfbdSBarry Smith } 1177c79dcfbdSBarry Smith 1178d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRollBack_ARKIMEX(TS ts) 1179d71ae5a4SJacob Faibussowitsch { 118024655328SShri TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 118124655328SShri ARKTableau tab = ark->tableau; 118224655328SShri const PetscInt s = tab->s; 118324655328SShri const PetscReal *bt = tab->bt, *b = tab->b; 118424655328SShri PetscScalar *w = ark->work; 118524655328SShri Vec *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS; 118624655328SShri PetscInt j; 1187be5899b3SLisandro Dalcin PetscReal h; 118824655328SShri 118924655328SShri PetscFunctionBegin; 1190be5899b3SLisandro Dalcin switch (ark->status) { 1191be5899b3SLisandro Dalcin case TS_STEP_INCOMPLETE: 1192d71ae5a4SJacob Faibussowitsch case TS_STEP_PENDING: 1193d71ae5a4SJacob Faibussowitsch h = ts->time_step; 1194d71ae5a4SJacob Faibussowitsch break; 1195d71ae5a4SJacob Faibussowitsch case TS_STEP_COMPLETE: 1196d71ae5a4SJacob Faibussowitsch h = ts->ptime - ts->ptime_prev; 1197d71ae5a4SJacob Faibussowitsch break; 1198d71ae5a4SJacob Faibussowitsch default: 1199d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1200be5899b3SLisandro Dalcin } 120124655328SShri for (j = 0; j < s; j++) w[j] = -h * bt[j]; 12029566063dSJacob Faibussowitsch PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotI)); 1203*d27334e2SStefano Zampini if (tab->additive) { 1204*d27334e2SStefano Zampini PetscBool hasE; 1205*d27334e2SStefano Zampini 1206*d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1207*d27334e2SStefano Zampini if (hasE) { 120824655328SShri for (j = 0; j < s; j++) w[j] = -h * b[j]; 12099566063dSJacob Faibussowitsch PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS)); 1210*d27334e2SStefano Zampini } 1211*d27334e2SStefano Zampini } 12123ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 121324655328SShri } 121424655328SShri 1215d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSStep_ARKIMEX(TS ts) 1216d71ae5a4SJacob Faibussowitsch { 12178a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 12188a381b04SJed Brown ARKTableau tab = ark->tableau; 12198a381b04SJed Brown const PetscInt s = tab->s; 122024655328SShri const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 1221406d0ec2SJed Brown PetscScalar *w = ark->work; 12221297b224SEmil Constantinescu Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 122396400bd6SLisandro Dalcin PetscBool extrapolate = ark->extrapolate; 1224108c343cSJed Brown TSAdapt adapt; 12258a381b04SJed Brown SNES snes; 1226fecfb714SLisandro Dalcin PetscInt i, j, its, lits; 1227be5899b3SLisandro Dalcin PetscInt rejections = 0; 1228*d27334e2SStefano Zampini PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE; 122996400bd6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 12308a381b04SJed Brown 12318a381b04SJed Brown PetscFunctionBegin; 123296400bd6SLisandro Dalcin if (ark->extrapolate && !ark->Y_prev) { 12339566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 12349566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1235*d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 123696400bd6SLisandro Dalcin } 123796400bd6SLisandro Dalcin 1238*d27334e2SStefano Zampini if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE)); 1239*d27334e2SStefano Zampini if (!hasE) dirk = PETSC_TRUE; 1240*d27334e2SStefano Zampini 124196400bd6SLisandro Dalcin if (!ts->steprollback) { 1242*d27334e2SStefano Zampini if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 12439566063dSJacob Faibussowitsch PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 124496400bd6SLisandro Dalcin } 1245fecfb714SLisandro Dalcin if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 124696400bd6SLisandro Dalcin for (i = 0; i < s; i++) { 12479566063dSJacob Faibussowitsch PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 12489566063dSJacob Faibussowitsch PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 1249*d27334e2SStefano Zampini if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 125096400bd6SLisandro Dalcin } 125196400bd6SLisandro Dalcin } 125296400bd6SLisandro Dalcin } 125396400bd6SLisandro Dalcin 1254*d27334e2SStefano Zampini /* For fully implicit formulations we can solve the equations 1255*d27334e2SStefano Zampini F(tn,xn,xdot) = 0 1256*d27334e2SStefano Zampini for the explicit first stage */ 1257*d27334e2SStefano Zampini if (dirk && tab->explicit_first_stage && ts->steprestart) { 1258*d27334e2SStefano Zampini ark->scoeff = 0.0; 1259*d27334e2SStefano Zampini PetscCall(VecCopy(ts->vec_sol, Z)); 1260*d27334e2SStefano Zampini PetscCall(TSGetSNES(ts, &snes)); 1261*d27334e2SStefano Zampini PetscCall(SNESSolve(snes, NULL, Ydot0)); 1262*d27334e2SStefano Zampini } 1263*d27334e2SStefano Zampini 1264*d27334e2SStefano Zampini /* For IMEX we compute a step */ 1265*d27334e2SStefano Zampini if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 126696400bd6SLisandro Dalcin TS ts_start; 1267*d27334e2SStefano Zampini if (PetscDefined(USE_DEBUG) && hasE) { 1268c79dcfbdSBarry Smith PetscBool id = PETSC_FALSE; 12699566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1270*d27334e2SStefano Zampini PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "This scheme requires an identity mass matrix, however the TSIFunction you provided does not utilize an identity mass matrix"); 1271c79dcfbdSBarry Smith } 12729566063dSJacob Faibussowitsch PetscCall(TSClone(ts, &ts_start)); 12739566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 12749566063dSJacob Faibussowitsch PetscCall(TSSetTime(ts_start, ts->ptime)); 12759566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 12769566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 12779566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 12789566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 12799566063dSJacob Faibussowitsch PetscCall(TSSetType(ts_start, TSARKIMEX)); 12809566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 12819566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 128234561852SEmil Constantinescu 12839566063dSJacob Faibussowitsch PetscCall(TSRestartStep(ts_start)); 12849566063dSJacob Faibussowitsch PetscCall(TSSolve(ts_start, ts->vec_sol)); 12859566063dSJacob Faibussowitsch PetscCall(TSGetTime(ts_start, &ts->ptime)); 12869566063dSJacob Faibussowitsch PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 1287bbd56ea5SKarl Rupp 128885fc7851SLisandro Dalcin { /* Save the initial slope for the next step */ 128985fc7851SLisandro Dalcin TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 12909566063dSJacob Faibussowitsch PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 129185fc7851SLisandro Dalcin } 129296400bd6SLisandro Dalcin ts->steps++; 129334561852SEmil Constantinescu 1294d15a3a53SEmil Constantinescu /* Set the correct TS in SNES */ 1295d15a3a53SEmil Constantinescu /* We'll try to bypass this by changing the method on the fly */ 129696400bd6SLisandro Dalcin { 12979566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 12989566063dSJacob Faibussowitsch PetscCall(TSSetSNES(ts, snes)); 129996400bd6SLisandro Dalcin } 13009566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts_start)); 1301e817cc15SEmil Constantinescu } 1302e817cc15SEmil Constantinescu 1303108c343cSJed Brown ark->status = TS_STEP_INCOMPLETE; 130496400bd6SLisandro Dalcin while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 130596400bd6SLisandro Dalcin PetscReal t = ts->ptime; 1306108c343cSJed Brown PetscReal h = ts->time_step; 13078a381b04SJed Brown for (i = 0; i < s; i++) { 13089be3e283SDebojyoti Ghosh ark->stage_time = t + h * ct[i]; 13099566063dSJacob Faibussowitsch PetscCall(TSPreStage(ts, ark->stage_time)); 13108a381b04SJed Brown if (At[i * s + i] == 0) { /* This stage is explicit */ 13113c633725SBarry 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"); 13129566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, Y[i])); 1313e817cc15SEmil Constantinescu for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 13149566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 1315*d27334e2SStefano Zampini if (tab->additive && hasE) { 13168a381b04SJed Brown for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 13179566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 1318*d27334e2SStefano Zampini } 13198a381b04SJed Brown } else { 1320b296d7d5SJed Brown ark->scoeff = 1. / At[i * s + i]; 13218a381b04SJed Brown /* Ydot = shift*(Y-Z) */ 13229566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, Z)); 1323e817cc15SEmil Constantinescu for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 13249566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Z, i, w, YdotI)); 1325*d27334e2SStefano Zampini if (tab->additive && hasE) { 1326c58d1302SEmil Constantinescu for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 13279566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 1328*d27334e2SStefano Zampini } 1329fecfb714SLisandro Dalcin if (extrapolate && !ts->steprestart) { 133056dcabbaSDebojyoti Ghosh /* Initial guess extrapolated from previous time step stage values */ 13319566063dSJacob Faibussowitsch PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 133256dcabbaSDebojyoti Ghosh } else { 13338a381b04SJed Brown /* Initial guess taken from last stage */ 13349566063dSJacob Faibussowitsch PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 133556dcabbaSDebojyoti Ghosh } 13369566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 13379566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, Y[i])); 13389566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its)); 13399566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 13409371c9d4SSatish Balay ts->snes_its += its; 13419371c9d4SSatish Balay ts->ksp_its += lits; 13429566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 13439566063dSJacob Faibussowitsch PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 134496400bd6SLisandro Dalcin if (!stageok) { 13451be93e3eSJed Brown /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 13461be93e3eSJed Brown * use extrapolation to initialize the solves on the next attempt. */ 134796400bd6SLisandro Dalcin extrapolate = PETSC_FALSE; 13481be93e3eSJed Brown goto reject_step; 13491be93e3eSJed Brown } 13508a381b04SJed Brown } 1351*d27334e2SStefano Zampini if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { 1352e817cc15SEmil Constantinescu if (i == 0 && tab->explicit_first_stage) { 1353*d27334e2SStefano 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", 1354*d27334e2SStefano Zampini ((PetscObject)ts)->type_name, ark->tableau->name); 13559566063dSJacob Faibussowitsch PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 1356e817cc15SEmil Constantinescu } else { 13579566063dSJacob Faibussowitsch PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1358e817cc15SEmil Constantinescu } 1359e817cc15SEmil Constantinescu } else { 13605eca1a21SEmil Constantinescu if (i == 0 && tab->explicit_first_stage) { 13619566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Ydot)); 13629566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 13639566063dSJacob Faibussowitsch PetscCall(VecScale(YdotI[i], -1.0)); 13645eca1a21SEmil Constantinescu } else { 13659566063dSJacob Faibussowitsch PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 13665eca1a21SEmil Constantinescu } 1367*d27334e2SStefano Zampini if (hasE) { 13684cc180ffSJed Brown if (ark->imex) { 13699566063dSJacob Faibussowitsch PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 13704cc180ffSJed Brown } else { 13719566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(YdotRHS[i])); 13724cc180ffSJed Brown } 13738a381b04SJed Brown } 1374*d27334e2SStefano Zampini } 13759566063dSJacob Faibussowitsch PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 1376e817cc15SEmil Constantinescu } 137796400bd6SLisandro Dalcin 1378be5899b3SLisandro Dalcin ark->status = TS_STEP_INCOMPLETE; 13799566063dSJacob Faibussowitsch PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 1380108c343cSJed Brown ark->status = TS_STEP_PENDING; 13819566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 13829566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidatesClear(adapt)); 13839566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 13849566063dSJacob Faibussowitsch PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 138596400bd6SLisandro Dalcin ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 138696400bd6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 13879566063dSJacob Faibussowitsch PetscCall(TSRollBack_ARKIMEX(ts)); 1388be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1389be5899b3SLisandro Dalcin goto reject_step; 139096400bd6SLisandro Dalcin } 139196400bd6SLisandro Dalcin 13928a381b04SJed Brown ts->ptime += ts->time_step; 1393cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 1394108c343cSJed Brown break; 139596400bd6SLisandro Dalcin 139696400bd6SLisandro Dalcin reject_step: 13979371c9d4SSatish Balay ts->reject++; 13989371c9d4SSatish Balay accept = PETSC_FALSE; 1399be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 140096400bd6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 140163a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1402108c343cSJed Brown } 1403f85781f1SEmil Constantinescu } 14043ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 14058a381b04SJed Brown } 1406*d27334e2SStefano Zampini 140780ab5e31SHong Zhang /* 140880ab5e31SHong Zhang This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 140980ab5e31SHong Zhang Udot = H(t,U) + G(t,U) 141080ab5e31SHong Zhang This corresponds to F(t,U,Udot) = Udot-H(t,U). 141180ab5e31SHong Zhang 141280ab5e31SHong Zhang The complete adjoint equations are 141380ab5e31SHong Zhang (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 141480ab5e31SHong Zhang (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 141580ab5e31SHong Zhang + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0 141680ab5e31SHong Zhang lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 141780ab5e31SHong Zhang mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 141880ab5e31SHong Zhang + (b_i dGdP + bt[i] dHdP) lambda_{n+1} + dGdP \sum_{j=i+1}^s a[j][i] lambda_s[j] 141980ab5e31SHong Zhang + dHdP \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0 142080ab5e31SHong Zhang mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 142180ab5e31SHong Zhang where shift = 1/(at[i][i]*h) 142280ab5e31SHong Zhang 142380ab5e31SHong Zhang If at[i][i] is 0, the first equation falls back to 142480ab5e31SHong Zhang lambda_s[i] = h ( 142580ab5e31SHong Zhang (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 142680ab5e31SHong Zhang + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 142780ab5e31SHong Zhang 142880ab5e31SHong Zhang */ 142980ab5e31SHong Zhang static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 143080ab5e31SHong Zhang { 143180ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 143280ab5e31SHong Zhang ARKTableau tab = ark->tableau; 143380ab5e31SHong Zhang const PetscInt s = tab->s; 143480ab5e31SHong Zhang const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 143580ab5e31SHong Zhang PetscScalar *w = ark->work; 143680ab5e31SHong Zhang Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 143780ab5e31SHong Zhang Mat Jex, Jim, Jimpre; 143880ab5e31SHong Zhang PetscInt i, j, nadj; 143980ab5e31SHong Zhang PetscReal t = ts->ptime, stage_time_ex; 144080ab5e31SHong Zhang PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 144180ab5e31SHong Zhang KSP ksp; 144280ab5e31SHong Zhang 144380ab5e31SHong Zhang PetscFunctionBegin; 144480ab5e31SHong Zhang ark->status = TS_STEP_INCOMPLETE; 144580ab5e31SHong Zhang PetscCall(SNESGetKSP(ts->snes, &ksp)); 144680ab5e31SHong Zhang PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 144780ab5e31SHong Zhang PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 144880ab5e31SHong Zhang 144980ab5e31SHong Zhang for (i = s - 1; i >= 0; i--) { 145080ab5e31SHong Zhang ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 145180ab5e31SHong Zhang stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 145280ab5e31SHong Zhang if (At[i * s + i] == 0) { // This stage is explicit 145380ab5e31SHong Zhang ark->scoeff = 0.; 145480ab5e31SHong Zhang } else { 145580ab5e31SHong Zhang ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 145680ab5e31SHong Zhang } 145780ab5e31SHong Zhang PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 145880ab5e31SHong Zhang PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 145980ab5e31SHong Zhang if (ts->vecs_sensip) { 146080ab5e31SHong 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 146180ab5e31SHong Zhang PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 146280ab5e31SHong Zhang } 146380ab5e31SHong Zhang for (nadj = 0; nadj < ts->numcost; nadj++) { 146480ab5e31SHong Zhang /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 146580ab5e31SHong Zhang if (s - i - 1 > 0) { 146680ab5e31SHong Zhang /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 146780ab5e31SHong Zhang for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 146880ab5e31SHong Zhang PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 146980ab5e31SHong Zhang PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 147080ab5e31SHong Zhang /* (shift I - dHdU) Temp */ 147180ab5e31SHong Zhang PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 147280ab5e31SHong Zhang /* cancel out shift Temp where shift=-scoeff/h */ 147380ab5e31SHong Zhang PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 147480ab5e31SHong Zhang if (ts->vecs_sensip) { 147580ab5e31SHong Zhang /* - dHdP Temp */ 147680ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 147780ab5e31SHong Zhang /* mu_n += h dHdP Temp */ 147880ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsSensiPTemp[nadj])); 147980ab5e31SHong Zhang } 148080ab5e31SHong Zhang /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 148180ab5e31SHong Zhang for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 148280ab5e31SHong Zhang PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 148380ab5e31SHong Zhang PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 148480ab5e31SHong Zhang /* dGdU Temp */ 148580ab5e31SHong Zhang PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 148680ab5e31SHong Zhang if (ts->vecs_sensip) { 148780ab5e31SHong Zhang /* dGdP Temp */ 148880ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 148980ab5e31SHong Zhang /* mu_n += h dGdP Temp */ 149080ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 149180ab5e31SHong Zhang } 149280ab5e31SHong Zhang } else { 149380ab5e31SHong Zhang PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 149480ab5e31SHong Zhang } 149580ab5e31SHong Zhang if (bt[i]) { // bt[i] dHdU lambda_{n+1} 149680ab5e31SHong Zhang /* (shift I - dHdU)^T lambda_{n+1} */ 149780ab5e31SHong Zhang PetscCall(MatMultTranspose(Jim, ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); 149880ab5e31SHong Zhang /* -bt[i] (shift I - dHdU)^T lambda_{n+1} */ 149980ab5e31SHong Zhang PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i], VecsSensiTemp[nadj])); 150080ab5e31SHong Zhang /* cancel out -bt[i] shift lambda_{n+1} */ 150180ab5e31SHong Zhang PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i] * ark->scoeff / adjoint_time_step, ts->vecs_sensi[nadj])); 150280ab5e31SHong Zhang if (ts->vecs_sensip) { 150380ab5e31SHong Zhang /* -dHdP lambda_{n+1} */ 150480ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacp, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj])); 150580ab5e31SHong Zhang /* mu_n += h bt[i] dHdP lambda_{n+1} */ 150680ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], -bt[i] * adjoint_time_step, VecsSensiPTemp[nadj])); 150780ab5e31SHong Zhang } 150880ab5e31SHong Zhang } 150980ab5e31SHong Zhang if (b[i]) { // b[i] dGdU lambda_{n+1} 151080ab5e31SHong Zhang /* dGdU lambda_{n+1} */ 151180ab5e31SHong Zhang PetscCall(MatMultTranspose(Jex, ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); 151280ab5e31SHong Zhang /* b[i] dGdU lambda_{n+1} */ 151380ab5e31SHong Zhang PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], b[i], VecsSensiTemp[nadj])); 151480ab5e31SHong Zhang if (ts->vecs_sensip) { 151580ab5e31SHong Zhang /* dGdP lambda_{n+1} */ 151680ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacprhs, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj])); 151780ab5e31SHong Zhang /* mu_n += h b[i] dGdP lambda_{n+1} */ 151880ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], b[i] * adjoint_time_step, VecsSensiPTemp[nadj])); 151980ab5e31SHong Zhang } 152080ab5e31SHong Zhang } 152180ab5e31SHong Zhang /* Build LHS for first-order adjoint */ 152280ab5e31SHong Zhang if (At[i * s + i] == 0) { // This stage is explicit 152380ab5e31SHong Zhang PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 152480ab5e31SHong Zhang } else { 152580ab5e31SHong Zhang KSP ksp; 152680ab5e31SHong Zhang KSPConvergedReason kspreason; 152780ab5e31SHong Zhang PetscCall(SNESGetKSP(ts->snes, &ksp)); 152880ab5e31SHong Zhang PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 152980ab5e31SHong Zhang PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 153080ab5e31SHong Zhang PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 153180ab5e31SHong Zhang PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 153280ab5e31SHong Zhang if (kspreason < 0) { 153380ab5e31SHong Zhang ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 153480ab5e31SHong 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)); 153580ab5e31SHong Zhang } 153680ab5e31SHong Zhang if (ts->vecs_sensip) { 153780ab5e31SHong Zhang /* -dHdP lambda_s[i] */ 153880ab5e31SHong Zhang PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 153980ab5e31SHong Zhang /* mu_n += h at[i][i] dHdP lambda_s[i] */ 154080ab5e31SHong Zhang PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 154180ab5e31SHong Zhang } 154280ab5e31SHong Zhang } 154380ab5e31SHong Zhang } 154480ab5e31SHong Zhang } 154580ab5e31SHong Zhang for (j = 0; j < s; j++) w[j] = 1.0; 154680ab5e31SHong Zhang for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 154780ab5e31SHong Zhang PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 154880ab5e31SHong Zhang ark->status = TS_STEP_COMPLETE; 154980ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 155080ab5e31SHong Zhang } 15518a381b04SJed Brown 1552d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1553d71ae5a4SJacob Faibussowitsch { 1554cd652676SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1555*d27334e2SStefano Zampini ARKTableau tab = ark->tableau; 1556*d27334e2SStefano Zampini PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1557108c343cSJed Brown PetscReal h; 1558108c343cSJed Brown PetscReal tt, t; 1559*d27334e2SStefano Zampini PetscScalar *bt = ark->work, *b = ark->work + s; 1560*d27334e2SStefano Zampini const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1561cd652676SJed Brown 1562cd652676SJed Brown PetscFunctionBegin; 1563*d27334e2SStefano 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); 1564108c343cSJed Brown switch (ark->status) { 1565108c343cSJed Brown case TS_STEP_INCOMPLETE: 1566108c343cSJed Brown case TS_STEP_PENDING: 1567108c343cSJed Brown h = ts->time_step; 1568108c343cSJed Brown t = (itime - ts->ptime) / h; 1569108c343cSJed Brown break; 1570108c343cSJed Brown case TS_STEP_COMPLETE: 1571be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1572108c343cSJed Brown t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1573108c343cSJed Brown break; 1574d71ae5a4SJacob Faibussowitsch default: 1575d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1576108c343cSJed Brown } 1577cd652676SJed Brown for (i = 0; i < s; i++) bt[i] = b[i] = 0; 15784f385281SJed Brown for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1579cd652676SJed Brown for (i = 0; i < s; i++) { 1580c1758d98SDebojyoti Ghosh bt[i] += h * Bt[i * pinterp + j] * tt; 1581108c343cSJed Brown b[i] += h * B[i * pinterp + j] * tt; 1582cd652676SJed Brown } 1583cd652676SJed Brown } 15849566063dSJacob Faibussowitsch PetscCall(VecCopy(ark->Y[0], X)); 15859566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1586*d27334e2SStefano Zampini if (tab->additive) { 1587*d27334e2SStefano Zampini PetscBool hasE; 1588*d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1589*d27334e2SStefano Zampini if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1590*d27334e2SStefano Zampini } 15913ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1592cd652676SJed Brown } 1593cd652676SJed Brown 1594d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1595d71ae5a4SJacob Faibussowitsch { 159656dcabbaSDebojyoti Ghosh TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1597*d27334e2SStefano Zampini ARKTableau tab = ark->tableau; 1598*d27334e2SStefano Zampini PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1599be5899b3SLisandro Dalcin PetscReal h, h_prev, t, tt; 1600*d27334e2SStefano Zampini PetscScalar *bt = ark->work, *b = ark->work + s; 1601*d27334e2SStefano Zampini const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 160256dcabbaSDebojyoti Ghosh 160356dcabbaSDebojyoti Ghosh PetscFunctionBegin; 16043c633725SBarry Smith PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 160581d12688SDebojyoti Ghosh h = ts->time_step; 1606be5899b3SLisandro Dalcin h_prev = ts->ptime - ts->ptime_prev; 1607be5899b3SLisandro Dalcin t = 1 + h / h_prev * c; 1608*d27334e2SStefano Zampini for (i = 0; i < s; i++) bt[i] = b[i] = 0; 160956dcabbaSDebojyoti Ghosh for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 161056dcabbaSDebojyoti Ghosh for (i = 0; i < s; i++) { 161181d12688SDebojyoti Ghosh bt[i] += h * Bt[i * pinterp + j] * tt; 161256dcabbaSDebojyoti Ghosh b[i] += h * B[i * pinterp + j] * tt; 161356dcabbaSDebojyoti Ghosh } 161456dcabbaSDebojyoti Ghosh } 16153c633725SBarry Smith PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 16169566063dSJacob Faibussowitsch PetscCall(VecCopy(ark->Y_prev[0], X)); 16179566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1618*d27334e2SStefano Zampini if (tab->additive) { 1619*d27334e2SStefano Zampini PetscBool hasE; 1620*d27334e2SStefano Zampini PetscCall(TSHasRHSFunction(ts, &hasE)); 1621*d27334e2SStefano Zampini if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1622*d27334e2SStefano Zampini } 16233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 162456dcabbaSDebojyoti Ghosh } 162556dcabbaSDebojyoti Ghosh 1626d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1627d71ae5a4SJacob Faibussowitsch { 162896400bd6SLisandro Dalcin TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 162996400bd6SLisandro Dalcin ARKTableau tab = ark->tableau; 163096400bd6SLisandro Dalcin 163196400bd6SLisandro Dalcin PetscFunctionBegin; 16323ba16761SJacob Faibussowitsch if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 16339566063dSJacob Faibussowitsch PetscCall(PetscFree(ark->work)); 16349566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 16359566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 16369566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 16379566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 16389566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 16399566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 16403ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 164196400bd6SLisandro Dalcin } 164296400bd6SLisandro Dalcin 1643d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSReset_ARKIMEX(TS ts) 1644d71ae5a4SJacob Faibussowitsch { 16458a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 16468a381b04SJed Brown 16478a381b04SJed Brown PetscFunctionBegin; 16489566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTableauReset(ts)); 16499566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ark->Ydot)); 16509566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ark->Ydot0)); 16519566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ark->Z)); 16523ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 16538a381b04SJed Brown } 16548a381b04SJed Brown 165580ab5e31SHong Zhang static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 165680ab5e31SHong Zhang { 165780ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 165880ab5e31SHong Zhang ARKTableau tab = ark->tableau; 165980ab5e31SHong Zhang 166080ab5e31SHong Zhang PetscFunctionBegin; 166180ab5e31SHong Zhang PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 166280ab5e31SHong Zhang PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 166380ab5e31SHong Zhang PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 166480ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 166580ab5e31SHong Zhang } 166680ab5e31SHong Zhang 1667d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1668d71ae5a4SJacob Faibussowitsch { 1669d5e6173cSPeter Brune TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1670d5e6173cSPeter Brune 1671d5e6173cSPeter Brune PetscFunctionBegin; 1672d5e6173cSPeter Brune if (Z) { 1673d5e6173cSPeter Brune if (dm && dm != ts->dm) { 16749566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1675d5e6173cSPeter Brune } else *Z = ax->Z; 1676d5e6173cSPeter Brune } 1677d5e6173cSPeter Brune if (Ydot) { 1678d5e6173cSPeter Brune if (dm && dm != ts->dm) { 16799566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1680d5e6173cSPeter Brune } else *Ydot = ax->Ydot; 1681d5e6173cSPeter Brune } 16823ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1683d5e6173cSPeter Brune } 1684d5e6173cSPeter Brune 1685d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1686d71ae5a4SJacob Faibussowitsch { 1687d5e6173cSPeter Brune PetscFunctionBegin; 1688d5e6173cSPeter Brune if (Z) { 168948a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1690d5e6173cSPeter Brune } 1691d5e6173cSPeter Brune if (Ydot) { 169248a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1693d5e6173cSPeter Brune } 16943ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1695d5e6173cSPeter Brune } 1696d5e6173cSPeter Brune 16978a381b04SJed Brown /* 16988a381b04SJed Brown This defines the nonlinear equation that is to be solved with SNES 16998a381b04SJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 17008a381b04SJed Brown */ 1701d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1702d71ae5a4SJacob Faibussowitsch { 17038a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1704d5e6173cSPeter Brune DM dm, dmsave; 1705d5e6173cSPeter Brune Vec Z, Ydot; 1706b296d7d5SJed Brown PetscReal shift = ark->scoeff / ts->time_step; 17078a381b04SJed Brown 17088a381b04SJed Brown PetscFunctionBegin; 17099566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 17109566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1711d5e6173cSPeter Brune dmsave = ts->dm; 1712d5e6173cSPeter Brune ts->dm = dm; 1713740132f1SEmil Constantinescu 1714*d27334e2SStefano Zampini if (ark->scoeff == 0.0) { 1715*d27334e2SStefano Zampini /* We are solving F(t,x_n,xdot) = 0 to start the method */ 1716*d27334e2SStefano Zampini PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 1717*d27334e2SStefano Zampini } else { 1718*d27334e2SStefano Zampini PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 17199566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1720*d27334e2SStefano Zampini } 1721e817cc15SEmil Constantinescu 1722d5e6173cSPeter Brune ts->dm = dmsave; 17239566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 17243ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 17258a381b04SJed Brown } 17268a381b04SJed Brown 1727d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1728d71ae5a4SJacob Faibussowitsch { 17298a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1730d5e6173cSPeter Brune DM dm, dmsave; 1731*d27334e2SStefano Zampini Vec Ydot, Z; 1732b296d7d5SJed Brown PetscReal shift = ark->scoeff / ts->time_step; 17338a381b04SJed Brown 17348a381b04SJed Brown PetscFunctionBegin; 17359566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 1736*d27334e2SStefano Zampini PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 17378a381b04SJed Brown /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1738d5e6173cSPeter Brune dmsave = ts->dm; 1739d5e6173cSPeter Brune ts->dm = dm; 1740740132f1SEmil Constantinescu 1741*d27334e2SStefano Zampini if (ark->scoeff == 0.0) { 1742*d27334e2SStefano Zampini /* We are solving F(t,x_n,xdot) = 0 to start the method, we only only dF/dXdot 1743*d27334e2SStefano Zampini Jed's proposal is to compute with a very large shift and scale back the matrix */ 1744*d27334e2SStefano Zampini shift = 1.0 / PETSC_MACHINE_EPSILON; 1745*d27334e2SStefano Zampini PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1746*d27334e2SStefano Zampini PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 1747*d27334e2SStefano Zampini if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1748*d27334e2SStefano Zampini } else { 17499566063dSJacob Faibussowitsch PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1750*d27334e2SStefano Zampini } 1751d5e6173cSPeter Brune ts->dm = dmsave; 1752*d27334e2SStefano Zampini PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 17533ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1754d5e6173cSPeter Brune } 1755d5e6173cSPeter Brune 175680ab5e31SHong Zhang static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 175780ab5e31SHong Zhang { 175880ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 175980ab5e31SHong Zhang 176080ab5e31SHong Zhang PetscFunctionBegin; 176180ab5e31SHong Zhang if (ns) *ns = ark->tableau->s; 176280ab5e31SHong Zhang if (Y) *Y = ark->Y; 176380ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 176480ab5e31SHong Zhang } 176580ab5e31SHong Zhang 1766d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1767d71ae5a4SJacob Faibussowitsch { 1768d5e6173cSPeter Brune PetscFunctionBegin; 17693ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1770d5e6173cSPeter Brune } 1771d5e6173cSPeter Brune 1772d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1773d71ae5a4SJacob Faibussowitsch { 1774d5e6173cSPeter Brune TS ts = (TS)ctx; 1775d5e6173cSPeter Brune Vec Z, Z_c; 1776d5e6173cSPeter Brune 1777d5e6173cSPeter Brune PetscFunctionBegin; 17789566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 17799566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 17809566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Z, Z_c)); 17819566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 17829566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 17839566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 17843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 17858a381b04SJed Brown } 17868a381b04SJed Brown 1787d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1788d71ae5a4SJacob Faibussowitsch { 1789cdb298fcSPeter Brune PetscFunctionBegin; 17903ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1791cdb298fcSPeter Brune } 1792cdb298fcSPeter Brune 1793d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1794d71ae5a4SJacob Faibussowitsch { 1795cdb298fcSPeter Brune TS ts = (TS)ctx; 1796cdb298fcSPeter Brune Vec Z, Z_c; 1797cdb298fcSPeter Brune 1798cdb298fcSPeter Brune PetscFunctionBegin; 17999566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 18009566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 1801cdb298fcSPeter Brune 18029566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 18039566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1804cdb298fcSPeter Brune 18059566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 18069566063dSJacob Faibussowitsch PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 18073ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1808cdb298fcSPeter Brune } 1809cdb298fcSPeter Brune 1810d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 1811d71ae5a4SJacob Faibussowitsch { 181296400bd6SLisandro Dalcin TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 181396400bd6SLisandro Dalcin ARKTableau tab = ark->tableau; 181496400bd6SLisandro Dalcin 181596400bd6SLisandro Dalcin PetscFunctionBegin; 1816*d27334e2SStefano Zampini PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 18179566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 18189566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 1819*d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 182096400bd6SLisandro Dalcin if (ark->extrapolate) { 18219566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 18229566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1823*d27334e2SStefano Zampini if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 182496400bd6SLisandro Dalcin } 18253ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 182696400bd6SLisandro Dalcin } 182796400bd6SLisandro Dalcin 1828d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 1829d71ae5a4SJacob Faibussowitsch { 18308a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1831d5e6173cSPeter Brune DM dm; 183296400bd6SLisandro Dalcin SNES snes; 1833f9c1d6abSBarry Smith 18348a381b04SJed Brown PetscFunctionBegin; 18359566063dSJacob Faibussowitsch PetscCall(TSARKIMEXTableauSetUp(ts)); 18369566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 18379566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 18389566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 18399566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm)); 18409566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 18419566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 18429566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 18433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 18448a381b04SJed Brown } 18458a381b04SJed Brown 184680ab5e31SHong Zhang static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 184780ab5e31SHong Zhang { 184880ab5e31SHong Zhang TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 184980ab5e31SHong Zhang ARKTableau tab = ark->tableau; 185080ab5e31SHong Zhang 185180ab5e31SHong Zhang PetscFunctionBegin; 185280ab5e31SHong Zhang PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 185380ab5e31SHong Zhang PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 185480ab5e31SHong Zhang if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); } 185580ab5e31SHong Zhang if (PetscDefined(USE_DEBUG)) { 185680ab5e31SHong Zhang PetscBool id = PETSC_FALSE; 185780ab5e31SHong Zhang PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1858*d27334e2SStefano Zampini PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Adjoint ARKIMEX requires an identity mass matrix, however the TSIFunction you provided does not utilize an identity mass matrix"); 185980ab5e31SHong Zhang } 186080ab5e31SHong Zhang PetscFunctionReturn(PETSC_SUCCESS); 186180ab5e31SHong Zhang } 186280ab5e31SHong Zhang 1863d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject) 1864d71ae5a4SJacob Faibussowitsch { 18654cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1866*d27334e2SStefano Zampini PetscBool dirk; 18678a381b04SJed Brown 18688a381b04SJed Brown PetscFunctionBegin; 1869*d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 1870*d27334e2SStefano Zampini PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 18718a381b04SJed Brown { 18728a381b04SJed Brown ARKTableauLink link; 18738a381b04SJed Brown PetscInt count, choice; 18748a381b04SJed Brown PetscBool flg; 18758a381b04SJed Brown const char **namelist; 1876*d27334e2SStefano Zampini for (link = ARKTableauList, count = 0; link; link = link->next) { 1877*d27334e2SStefano Zampini if (!dirk && link->tab.additive) count++; 1878*d27334e2SStefano Zampini if (dirk && !link->tab.additive) count++; 1879*d27334e2SStefano Zampini } 18809566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1881*d27334e2SStefano Zampini for (link = ARKTableauList, count = 0; link; link = link->next) { 1882*d27334e2SStefano Zampini if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 1883*d27334e2SStefano Zampini if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 1884*d27334e2SStefano Zampini } 1885*d27334e2SStefano Zampini if (dirk) { 1886*d27334e2SStefano Zampini PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 1887*d27334e2SStefano Zampini if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 1888*d27334e2SStefano Zampini } else { 18899566063dSJacob Faibussowitsch PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 18909566063dSJacob Faibussowitsch if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 18914cc180ffSJed Brown flg = (PetscBool)!ark->imex; 18929566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 18934cc180ffSJed Brown ark->imex = (PetscBool)!flg; 1894*d27334e2SStefano Zampini } 1895*d27334e2SStefano Zampini PetscCall(PetscFree(namelist)); 18969566063dSJacob 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)); 18978a381b04SJed Brown } 1898d0609cedSBarry Smith PetscOptionsHeadEnd(); 18993ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19008a381b04SJed Brown } 19018a381b04SJed Brown 1902d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 1903d71ae5a4SJacob Faibussowitsch { 19048a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1905*d27334e2SStefano Zampini PetscBool iascii, dirk; 19068a381b04SJed Brown 19078a381b04SJed Brown PetscFunctionBegin; 1908*d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 19099566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 19108a381b04SJed Brown if (iascii) { 1911*d27334e2SStefano Zampini PetscViewerFormat format; 19129c334d8fSLisandro Dalcin ARKTableau tab = ark->tableau; 191319fd82e9SBarry Smith TSARKIMEXType arktype; 1914*d27334e2SStefano Zampini char buf[2048]; 19153a28c0e5SStefano Zampini PetscBool flg; 19163a28c0e5SStefano Zampini 19179566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetType(ts, &arktype)); 19189566063dSJacob Faibussowitsch PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 1919*d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 19209566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 1921*d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 1922*d27334e2SStefano Zampini PetscCall(PetscViewerGetFormat(viewer, &format)); 1923*d27334e2SStefano Zampini if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1924*d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 1925*d27334e2SStefano Zampini for (PetscInt i = 0; i < tab->s; i++) { 1926*d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 1927*d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 1928*d27334e2SStefano Zampini } 1929*d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 1930*d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 1931*d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 1932*d27334e2SStefano Zampini PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 1933*d27334e2SStefano Zampini } 19349566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 19359566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 19369566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 19379566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 1938*d27334e2SStefano Zampini if (!dirk) { 1939*d27334e2SStefano Zampini PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 19409566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 19418a381b04SJed Brown } 1942*d27334e2SStefano Zampini } 19433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19448a381b04SJed Brown } 19458a381b04SJed Brown 1946d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 1947d71ae5a4SJacob Faibussowitsch { 1948f2c2a1b9SBarry Smith SNES snes; 19499c334d8fSLisandro Dalcin TSAdapt adapt; 1950f2c2a1b9SBarry Smith 1951f2c2a1b9SBarry Smith PetscFunctionBegin; 19529566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 19539566063dSJacob Faibussowitsch PetscCall(TSAdaptLoad(adapt, viewer)); 19549566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 19559566063dSJacob Faibussowitsch PetscCall(SNESLoad(snes, viewer)); 1956ad6bc421SBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 19579566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 19589566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 19593ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1960f2c2a1b9SBarry Smith } 1961f2c2a1b9SBarry Smith 19628a381b04SJed Brown /*@C 1963bcf0153eSBarry Smith TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 19648a381b04SJed Brown 196520f4b53cSBarry Smith Logically Collective 19668a381b04SJed Brown 1967d8d19677SJose E. Roman Input Parameters: 19688a381b04SJed Brown + ts - timestepping context 1969bcf0153eSBarry Smith - arktype - type of `TSARKIMEX` scheme 19708a381b04SJed Brown 1971bcf0153eSBarry Smith Options Database Key: 1972bcf0153eSBarry Smith . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 19739bd3e852SBarry Smith 19748a381b04SJed Brown Level: intermediate 19758a381b04SJed Brown 19761cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 1977db781477SPatrick Sanan `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 19788a381b04SJed Brown @*/ 1979d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 1980d71ae5a4SJacob Faibussowitsch { 19818a381b04SJed Brown PetscFunctionBegin; 19828a381b04SJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 19834f572ea9SToby Isaac PetscAssertPointer(arktype, 2); 1984cac4c232SBarry Smith PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 19853ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 19868a381b04SJed Brown } 19878a381b04SJed Brown 19888a381b04SJed Brown /*@C 1989bcf0153eSBarry Smith TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 19908a381b04SJed Brown 199120f4b53cSBarry Smith Logically Collective 19928a381b04SJed Brown 19938a381b04SJed Brown Input Parameter: 19948a381b04SJed Brown . ts - timestepping context 19958a381b04SJed Brown 19968a381b04SJed Brown Output Parameter: 1997bcf0153eSBarry Smith . arktype - type of `TSARKIMEX` scheme 19988a381b04SJed Brown 19998a381b04SJed Brown Level: intermediate 20008a381b04SJed Brown 200142747ad1SJacob Faibussowitsch .seealso: [](ch_ts), `TSARKIMEXc` 20028a381b04SJed Brown @*/ 2003d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2004d71ae5a4SJacob Faibussowitsch { 20058a381b04SJed Brown PetscFunctionBegin; 20068a381b04SJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2007cac4c232SBarry Smith PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 20083ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20098a381b04SJed Brown } 20108a381b04SJed Brown 201116353aafSBarry Smith /*@ 2012bcf0153eSBarry Smith TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 20134cc180ffSJed Brown 201420f4b53cSBarry Smith Logically Collective 20154cc180ffSJed Brown 2016d8d19677SJose E. Roman Input Parameters: 20174cc180ffSJed Brown + ts - timestepping context 2018bcf0153eSBarry Smith - flg - `PETSC_TRUE` for fully implicit 20194cc180ffSJed Brown 20204cc180ffSJed Brown Level: intermediate 20214cc180ffSJed Brown 20221cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 20234cc180ffSJed Brown @*/ 2024d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 2025d71ae5a4SJacob Faibussowitsch { 20264cc180ffSJed Brown PetscFunctionBegin; 20274cc180ffSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 20283a28c0e5SStefano Zampini PetscValidLogicalCollectiveBool(ts, flg, 2); 2029cac4c232SBarry Smith PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 20303ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20314cc180ffSJed Brown } 20324cc180ffSJed Brown 20333a28c0e5SStefano Zampini /*@ 20343a28c0e5SStefano Zampini TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 20353a28c0e5SStefano Zampini 203620f4b53cSBarry Smith Logically Collective 20373a28c0e5SStefano Zampini 20383a28c0e5SStefano Zampini Input Parameter: 20393a28c0e5SStefano Zampini . ts - timestepping context 20403a28c0e5SStefano Zampini 20417a7aea1fSJed Brown Output Parameter: 2042bcf0153eSBarry Smith . flg - `PETSC_TRUE` for fully implicit 20433a28c0e5SStefano Zampini 20443a28c0e5SStefano Zampini Level: intermediate 20453a28c0e5SStefano Zampini 20461cc06b55SBarry Smith .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 20473a28c0e5SStefano Zampini @*/ 2048d71ae5a4SJacob Faibussowitsch PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 2049d71ae5a4SJacob Faibussowitsch { 20503a28c0e5SStefano Zampini PetscFunctionBegin; 20513a28c0e5SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 20524f572ea9SToby Isaac PetscAssertPointer(flg, 2); 2053cac4c232SBarry Smith PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 20543ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20553a28c0e5SStefano Zampini } 20563a28c0e5SStefano Zampini 2057d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 2058d71ae5a4SJacob Faibussowitsch { 20598a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 20608a381b04SJed Brown 20618a381b04SJed Brown PetscFunctionBegin; 20628a381b04SJed Brown *arktype = ark->tableau->name; 20633ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20648a381b04SJed Brown } 2065*d27334e2SStefano Zampini 2066d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 2067d71ae5a4SJacob Faibussowitsch { 20688a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 20698a381b04SJed Brown PetscBool match; 20708a381b04SJed Brown ARKTableauLink link; 20718a381b04SJed Brown 20728a381b04SJed Brown PetscFunctionBegin; 20738a381b04SJed Brown if (ark->tableau) { 20749566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 20753ba16761SJacob Faibussowitsch if (match) PetscFunctionReturn(PETSC_SUCCESS); 20768a381b04SJed Brown } 20778a381b04SJed Brown for (link = ARKTableauList; link; link = link->next) { 20789566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 20798a381b04SJed Brown if (match) { 20809566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 20818a381b04SJed Brown ark->tableau = &link->tab; 20829566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 20832ffb9264SLisandro Dalcin ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 20843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20858a381b04SJed Brown } 20868a381b04SJed Brown } 208798921bdaSJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 20888a381b04SJed Brown } 2089e0877f53SBarry Smith 2090d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 2091d71ae5a4SJacob Faibussowitsch { 20924cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 20934cc180ffSJed Brown 20944cc180ffSJed Brown PetscFunctionBegin; 20954cc180ffSJed Brown ark->imex = (PetscBool)!flg; 20963ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20974cc180ffSJed Brown } 20988a381b04SJed Brown 2099d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 2100d71ae5a4SJacob Faibussowitsch { 21013a28c0e5SStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 21023a28c0e5SStefano Zampini 21033a28c0e5SStefano Zampini PetscFunctionBegin; 21043a28c0e5SStefano Zampini *flg = (PetscBool)!ark->imex; 21053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 21063a28c0e5SStefano Zampini } 21073a28c0e5SStefano Zampini 2108d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 2109d71ae5a4SJacob Faibussowitsch { 2110b3a6b972SJed Brown PetscFunctionBegin; 21119566063dSJacob Faibussowitsch PetscCall(TSReset_ARKIMEX(ts)); 2112b3a6b972SJed Brown if (ts->dm) { 21139566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 21149566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2115b3a6b972SJed Brown } 21169566063dSJacob Faibussowitsch PetscCall(PetscFree(ts->data)); 2117*d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL)); 2118*d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL)); 21199566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 21209566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 21219566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 21222e956fe4SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 21233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2124b3a6b972SJed Brown } 2125b3a6b972SJed Brown 21268a381b04SJed Brown /* ------------------------------------------------------------ */ 21278a381b04SJed Brown /*MC 2128c79dcfbdSBarry Smith TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 21298a381b04SJed Brown 2130fca742c7SJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2131fca742c7SJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2132bcf0153eSBarry Smith of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2133d0685a90SJed Brown 21348a381b04SJed Brown Level: beginner 21358a381b04SJed Brown 2136bcf0153eSBarry Smith Notes: 2137bcf0153eSBarry Smith The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 21388a381b04SJed Brown 2139bcf0153eSBarry Smith If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2140bcf0153eSBarry Smith 2141bcf0153eSBarry 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). 2142bcf0153eSBarry Smith 2143bcf0153eSBarry Smith Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2144bcf0153eSBarry Smith 21451cc06b55SBarry Smith .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2146bcf0153eSBarry Smith `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2147bcf0153eSBarry Smith `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 21488a381b04SJed Brown M*/ 2149d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2150d71ae5a4SJacob Faibussowitsch { 215180ab5e31SHong Zhang TS_ARKIMEX *ark; 2152*d27334e2SStefano Zampini PetscBool dirk; 21538a381b04SJed Brown 21548a381b04SJed Brown PetscFunctionBegin; 21559566063dSJacob Faibussowitsch PetscCall(TSARKIMEXInitializePackage()); 2156*d27334e2SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 21578a381b04SJed Brown 21588a381b04SJed Brown ts->ops->reset = TSReset_ARKIMEX; 215980ab5e31SHong Zhang ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 21608a381b04SJed Brown ts->ops->destroy = TSDestroy_ARKIMEX; 21618a381b04SJed Brown ts->ops->view = TSView_ARKIMEX; 2162f2c2a1b9SBarry Smith ts->ops->load = TSLoad_ARKIMEX; 21638a381b04SJed Brown ts->ops->setup = TSSetUp_ARKIMEX; 216480ab5e31SHong Zhang ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 21658a381b04SJed Brown ts->ops->step = TSStep_ARKIMEX; 2166cd652676SJed Brown ts->ops->interpolate = TSInterpolate_ARKIMEX; 2167108c343cSJed Brown ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 216824655328SShri ts->ops->rollback = TSRollBack_ARKIMEX; 21698a381b04SJed Brown ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 21708a381b04SJed Brown ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 21718a381b04SJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 217280ab5e31SHong Zhang ts->ops->getstages = TSGetStages_ARKIMEX; 217380ab5e31SHong Zhang ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 21748a381b04SJed Brown 2175efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 2176efd4aadfSBarry Smith 217780ab5e31SHong Zhang PetscCall(PetscNew(&ark)); 217880ab5e31SHong Zhang ts->data = (void *)ark; 2179*d27334e2SStefano Zampini ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 218080ab5e31SHong Zhang 218180ab5e31SHong Zhang ark->VecsDeltaLam = NULL; 218280ab5e31SHong Zhang ark->VecsSensiTemp = NULL; 218380ab5e31SHong Zhang ark->VecsSensiPTemp = NULL; 21848a381b04SJed Brown 21859566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2186*d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2187*d27334e2SStefano Zampini if (!dirk) { 21889566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 21899566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 21909566063dSJacob Faibussowitsch PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2191*d27334e2SStefano Zampini } 2192*d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2193*d27334e2SStefano Zampini } 2194*d27334e2SStefano Zampini 2195*d27334e2SStefano Zampini /* ------------------------------------------------------------ */ 2196*d27334e2SStefano Zampini 2197*d27334e2SStefano Zampini static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2198*d27334e2SStefano Zampini { 2199*d27334e2SStefano Zampini TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2200*d27334e2SStefano Zampini 2201*d27334e2SStefano Zampini PetscFunctionBegin; 2202*d27334e2SStefano Zampini PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2203*d27334e2SStefano Zampini PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2204*d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2205*d27334e2SStefano Zampini } 2206*d27334e2SStefano Zampini 2207*d27334e2SStefano Zampini /*@C 2208*d27334e2SStefano Zampini TSDIRKSetType - Set the type of `TSDIRK` scheme 2209*d27334e2SStefano Zampini 2210*d27334e2SStefano Zampini Logically Collective 2211*d27334e2SStefano Zampini 2212*d27334e2SStefano Zampini Input Parameters: 2213*d27334e2SStefano Zampini + ts - timestepping context 2214*d27334e2SStefano Zampini - dirktype - type of `TSDIRK` scheme 2215*d27334e2SStefano Zampini 2216*d27334e2SStefano Zampini Options Database Key: 2217*d27334e2SStefano Zampini . -ts_dirkimex_type - set `TSDIRK` scheme type 2218*d27334e2SStefano Zampini 2219*d27334e2SStefano Zampini Level: intermediate 2220*d27334e2SStefano Zampini 2221*d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2222*d27334e2SStefano Zampini @*/ 2223*d27334e2SStefano Zampini PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2224*d27334e2SStefano Zampini { 2225*d27334e2SStefano Zampini PetscFunctionBegin; 2226*d27334e2SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2227*d27334e2SStefano Zampini PetscAssertPointer(dirktype, 2); 2228*d27334e2SStefano Zampini PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2229*d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2230*d27334e2SStefano Zampini } 2231*d27334e2SStefano Zampini 2232*d27334e2SStefano Zampini /*@C 2233*d27334e2SStefano Zampini TSDIRKGetType - Get the type of `TSDIRK` scheme 2234*d27334e2SStefano Zampini 2235*d27334e2SStefano Zampini Logically Collective 2236*d27334e2SStefano Zampini 2237*d27334e2SStefano Zampini Input Parameter: 2238*d27334e2SStefano Zampini . ts - timestepping context 2239*d27334e2SStefano Zampini 2240*d27334e2SStefano Zampini Output Parameter: 2241*d27334e2SStefano Zampini . dirktype - type of `TSDIRK` scheme 2242*d27334e2SStefano Zampini 2243*d27334e2SStefano Zampini Level: intermediate 2244*d27334e2SStefano Zampini 2245*d27334e2SStefano Zampini .seealso: [](ch_ts), `TSDIRKSetType()` 2246*d27334e2SStefano Zampini @*/ 2247*d27334e2SStefano Zampini PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2248*d27334e2SStefano Zampini { 2249*d27334e2SStefano Zampini PetscFunctionBegin; 2250*d27334e2SStefano Zampini PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2251*d27334e2SStefano Zampini PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2252*d27334e2SStefano Zampini PetscFunctionReturn(PETSC_SUCCESS); 2253*d27334e2SStefano Zampini } 2254*d27334e2SStefano Zampini 2255*d27334e2SStefano Zampini /*MC 2256*d27334e2SStefano Zampini TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes 2257*d27334e2SStefano Zampini 2258*d27334e2SStefano Zampini Level: beginner 2259*d27334e2SStefano Zampini 2260*d27334e2SStefano Zampini Notes: 2261*d27334e2SStefano Zampini The default is `TSDIRKS212`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type 2262*d27334e2SStefano Zampini 2263*d27334e2SStefano Zampini .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2264*d27334e2SStefano Zampini M*/ 2265*d27334e2SStefano Zampini PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2266*d27334e2SStefano Zampini { 2267*d27334e2SStefano Zampini PetscFunctionBegin; 2268*d27334e2SStefano Zampini PetscCall(TSCreate_ARKIMEX(ts)); 2269*d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2270*d27334e2SStefano Zampini PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2271*d27334e2SStefano Zampini PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 22723ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 22738a381b04SJed Brown } 2274