1 /* 2 Code for timestepping with additive Runge-Kutta IMEX method or Diagonally Implicit Runge-Kutta methods. 3 4 Notes: 5 For ARK, the general system is written as 6 7 F(t,U,Udot) = G(t,U) 8 9 where F represents the stiff part of the physics and G represents the non-stiff part. 10 11 */ 12 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 13 #include <petscdm.h> 14 15 static TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3; 16 static TSDIRKType TSDIRKDefault = TSDIRKS212; 17 static PetscBool TSARKIMEXRegisterAllCalled; 18 static PetscBool TSARKIMEXPackageInitialized; 19 static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec); 20 21 typedef struct _ARKTableau *ARKTableau; 22 struct _ARKTableau { 23 char *name; 24 PetscBool additive; /* If False, it is a DIRK method */ 25 PetscInt order; /* Classical approximation order of the method */ 26 PetscInt s; /* Number of stages */ 27 PetscBool stiffly_accurate; /* The implicit part is stiffly accurate */ 28 PetscBool FSAL_implicit; /* The implicit part is FSAL */ 29 PetscBool explicit_first_stage; /* The implicit part has an explicit first stage */ 30 PetscInt pinterp; /* Interpolation order */ 31 PetscReal *At, *bt, *ct; /* Stiff tableau */ 32 PetscReal *A, *b, *c; /* Non-stiff tableau */ 33 PetscReal *bembedt, *bembed; /* Embedded formula of order one less (order-1) */ 34 PetscReal *binterpt, *binterp; /* Dense output formula */ 35 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 36 }; 37 typedef struct _ARKTableauLink *ARKTableauLink; 38 struct _ARKTableauLink { 39 struct _ARKTableau tab; 40 ARKTableauLink next; 41 }; 42 static ARKTableauLink ARKTableauList; 43 44 typedef struct { 45 ARKTableau tableau; 46 Vec *Y; /* States computed during the step */ 47 Vec *YdotI; /* Time derivatives for the stiff part */ 48 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 49 Vec *Y_prev; /* States computed during the previous time step */ 50 Vec *YdotI_prev; /* Time derivatives for the stiff part for the previous time step*/ 51 Vec *YdotRHS_prev; /* Function evaluations for the non-stiff part for the previous time step*/ 52 Vec Ydot0; /* Holds the slope from the previous step in FSAL case */ 53 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 54 Vec Z; /* Ydot = shift(Y-Z) */ 55 PetscScalar *work; /* Scalar work */ 56 PetscReal scoeff; /* shift = scoeff/dt */ 57 PetscReal stage_time; 58 PetscBool imex; 59 PetscBool extrapolate; /* Extrapolate initial guess from previous time-step stage values */ 60 TSStepStatus status; 61 62 /* context for sensitivity analysis */ 63 Vec *VecsDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */ 64 Vec *VecsSensiTemp; /* Vectors to be multiplied with Jacobian transpose */ 65 Vec *VecsSensiPTemp; /* Temporary Vectors to store JacobianP-transpose-vector product */ 66 } TS_ARKIMEX; 67 68 /*MC 69 TSARKIMEXARS122 - Second order ARK IMEX scheme. 70 71 This method has one explicit stage and one implicit stage. 72 73 Options Database Key: 74 . -ts_arkimex_type ars122 - set arkimex type to ars122 75 76 Level: advanced 77 78 References: 79 . * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 80 81 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 82 M*/ 83 84 /*MC 85 TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part. 86 87 This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu. 88 89 Options Database Key: 90 . -ts_arkimex_type a2 - set arkimex type to a2 91 92 Level: advanced 93 94 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 95 M*/ 96 97 /*MC 98 TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part. 99 100 This method has two implicit stages, and L-stable implicit scheme. 101 102 Options Database Key: 103 . -ts_arkimex_type l2 - set arkimex type to l2 104 105 Level: advanced 106 107 References: 108 . * - 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. 109 110 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 111 M*/ 112 113 /*MC 114 TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method. 115 116 This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used. 117 118 Options Database Key: 119 . -ts_arkimex_type 1bee - set arkimex type to 1bee 120 121 Level: advanced 122 123 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 124 M*/ 125 126 /*MC 127 TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part. 128 129 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. 130 131 Options Database Key: 132 . -ts_arkimex_type 2c - set arkimex type to 2c 133 134 Level: advanced 135 136 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 137 M*/ 138 139 /*MC 140 TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 141 142 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. 143 144 Options Database Key: 145 . -ts_arkimex_type 2d - set arkimex type to 2d 146 147 Level: advanced 148 149 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 150 M*/ 151 152 /*MC 153 TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 154 155 This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu. 156 157 Options Database Key: 158 . -ts_arkimex_type 2e - set arkimex type to 2e 159 160 Level: advanced 161 162 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 163 M*/ 164 165 /*MC 166 TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme. 167 168 This method has three implicit stages. 169 170 References: 171 . * - 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. 172 173 This method is referred to as SSP2-(3,3,2) in https://arxiv.org/abs/1110.4375 174 175 Options Database Key: 176 . -ts_arkimex_type prssp2 - set arkimex type to prssp2 177 178 Level: advanced 179 180 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 181 M*/ 182 183 /*MC 184 TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part. 185 186 This method has one explicit stage and three implicit stages. 187 188 Options Database Key: 189 . -ts_arkimex_type 3 - set arkimex type to 3 190 191 Level: advanced 192 193 References: 194 . * - Kennedy and Carpenter 2003. 195 196 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 197 M*/ 198 199 /*MC 200 TSARKIMEXARS443 - Third order ARK IMEX scheme. 201 202 This method has one explicit stage and four implicit stages. 203 204 Options Database Key: 205 . -ts_arkimex_type ars443 - set arkimex type to ars443 206 207 Level: advanced 208 209 References: 210 + * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 211 - * - This method is referred to as ARS(4,4,3) in https://arxiv.org/abs/1110.4375 212 213 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 214 M*/ 215 216 /*MC 217 TSARKIMEXBPR3 - Third order ARK IMEX scheme. 218 219 This method has one explicit stage and four implicit stages. 220 221 Options Database Key: 222 . -ts_arkimex_type bpr3 - set arkimex type to bpr3 223 224 Level: advanced 225 226 References: 227 . * - This method is referred to as ARK3 in https://arxiv.org/abs/1110.4375 228 229 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 230 M*/ 231 232 /*MC 233 TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part. 234 235 This method has one explicit stage and four implicit stages. 236 237 Options Database Key: 238 . -ts_arkimex_type 4 - set arkimex type to4 239 240 Level: advanced 241 242 References: 243 . * - Kennedy and Carpenter 2003. 244 245 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 246 M*/ 247 248 /*MC 249 TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part. 250 251 This method has one explicit stage and five implicit stages. 252 253 Options Database Key: 254 . -ts_arkimex_type 5 - set arkimex type to 5 255 256 Level: advanced 257 258 References: 259 . * - Kennedy and Carpenter 2003. 260 261 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 262 M*/ 263 264 static PetscErrorCode TSHasRHSFunction(TS ts, PetscBool *has) 265 { 266 TSRHSFunction func; 267 268 PetscFunctionBegin; 269 *has = PETSC_FALSE; 270 PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL)); 271 if (func) *has = PETSC_TRUE; 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 /*@C 276 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX` 277 278 Not Collective, but should be called by all processes which will need the schemes to be registered 279 280 Level: advanced 281 282 .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()` 283 @*/ 284 PetscErrorCode TSARKIMEXRegisterAll(void) 285 { 286 PetscFunctionBegin; 287 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 288 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 289 290 #define RC PetscRealConstant 291 #define s2 RC(1.414213562373095048802) /* PetscSqrtReal((PetscReal)2.0) */ 292 #define us2 RC(0.2928932188134524755992) /* 1.0-1.0/PetscSqrtReal((PetscReal)2.0); */ 293 294 /* Diagonally implicit methods */ 295 { 296 /* DIRK212, default of SUNDIALS */ 297 const PetscReal A[2][2] = { 298 {RC(1.0), RC(0.0)}, 299 {RC(-1.0), RC(1.0)} 300 }; 301 const PetscReal b[2] = {RC(0.5), RC(0.5)}; 302 const PetscReal bembed[2] = {RC(1.0), RC(0.0)}; 303 PetscCall(TSDIRKRegister(TSDIRKS212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 304 } 305 306 { 307 /* ESDIRK12 from https://arxiv.org/pdf/1803.01613.pdf */ 308 const PetscReal A[2][2] = { 309 {RC(0.0), RC(0.0)}, 310 {RC(0.0), RC(1.0)} 311 }; 312 const PetscReal b[2] = {RC(0.0), RC(1.0)}; 313 const PetscReal bembed[2] = {RC(0.5), RC(0.5)}; 314 PetscCall(TSDIRKRegister(TSDIRKES212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 315 } 316 317 { 318 /* ESDIRK213L[2]SA from KC16. 319 TR-BDF2 from Osea and Shampine */ 320 const PetscReal A[3][3] = { 321 {RC(0.0), RC(0.0), RC(0.0)}, 322 {us2, us2, RC(0.0)}, 323 {s2 / RC(4.0), s2 / RC(4.0), us2 }, 324 }; 325 const PetscReal b[3] = {s2 / RC(4.0), s2 / RC(4.0), us2}; 326 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 PetscCall(TSDIRKRegister(TSDIRKES213SAL, 2, 3, &A[0][0], b, NULL, bembed, 1, b)); 328 } 329 330 { 331 #define g RC(0.43586652150845899941601945) 332 #define g2 PetscSqr(g) 333 #define g3 g *g2 334 #define g4 PetscSqr(g2) 335 #define g5 g *g4 336 #define c3 RC(1.0) 337 #define a32 c3 *(c3 - RC(2.0) * g) / (RC(4.0) * g) 338 #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 #define b3 (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) / (RC(3.0) * c3 * (c3 - RC(2.0) * g)) 340 #if 0 341 /* This is for c3 = 3/5 */ 342 #define bh2 \ 343 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 #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 #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 #else 347 /* This is for c3 = 1.0 */ 348 #define bh2 a32 349 #define bh3 g 350 #define bh4 RC(0.0) 351 #endif 352 /* ESDIRK3(2I)4L[2]SA from KC16 with c3 = 1.0 */ 353 /* Given by Kvaerno https://link.springer.com/article/10.1023/b:bitn.0000046811.70614.38 */ 354 const PetscReal A[4][4] = { 355 {RC(0.0), RC(0.0), RC(0.0), RC(0.0)}, 356 {g, g, RC(0.0), RC(0.0)}, 357 {c3 - a32 - g, a32, g, RC(0.0)}, 358 {RC(1.0) - b2 - b3 - g, b2, b3, g }, 359 }; 360 const PetscReal b[4] = {RC(1.0) - b2 - b3 - g, b2, b3, g}; 361 const PetscReal bembed[4] = {RC(1.0) - bh2 - bh3 - bh4, bh2, bh3, bh4}; 362 PetscCall(TSDIRKRegister(TSDIRKES324SAL, 3, 4, &A[0][0], b, NULL, bembed, 1, b)); 363 #undef g 364 #undef g2 365 #undef g3 366 #undef c3 367 #undef a32 368 #undef b2 369 #undef b3 370 #undef bh2 371 #undef bh3 372 #undef bh4 373 } 374 375 { 376 /* ESDIRK3(2I)5L[2]SA from KC16 */ 377 const PetscReal A[5][5] = { 378 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 379 {RC(9.0) / RC(40.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0), RC(0.0) }, 380 {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 {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 {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 }; 384 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 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 PetscCall(TSDIRKRegister(TSDIRKES325SAL, 3, 5, &A[0][0], b, NULL, bembed, 1, b)); 387 } 388 389 { 390 // DIRK(6,6)[1]A[(7,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 391 const PetscReal A[7][7] = { 392 {RC(0.303487844706747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 393 {RC(-0.279756492709814), RC(0.500032236020747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 394 {RC(0.280583215743895), RC(-0.438560061586751), RC(0.217250734515736), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 395 {RC(-0.0677678738539846), RC(0.984312781232293), RC(-0.266720192540149), RC(0.2476680834526), RC(0.0), RC(0.0), RC(0.0) }, 396 {RC(0.125671616147993), RC(-0.995401751002415), RC(0.761333109549059), RC(-0.210281837202208), RC(0.866743712636936), RC(0.0), RC(0.0) }, 397 {RC(-0.368056238801488), RC(-0.999928082701516), RC(0.534734253232519), RC(-0.174856916279082), RC(0.615007160285509), RC(0.696549912132029), RC(0.0) }, 398 {RC(-0.00570546839653984), RC(-0.113110431835656), RC(-0.000965563207671587), RC(-0.000130490084629567), RC(0.00111737736895673), RC(-0.279385587378871), RC(0.618455906845342)} 399 }; 400 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 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 PetscCall(TSDIRKRegister(TSDIRK657A, 6, 7, &A[0][0], b, NULL, bembed, 1, b)); 403 } 404 { 405 // ESDIRK(8,6)[2]SA[(8,4)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 406 const PetscReal A[8][8] = { 407 {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 {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 {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 {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 {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 {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 {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 {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 }; 416 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 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 PetscCall(TSDIRKRegister(TSDIRKES648SA, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 419 } 420 { 421 // DIRK(8,6)[1]SAL[(8,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 422 const PetscReal A[8][8] = { 423 {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 {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 {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 {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 {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 {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 {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 {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 }; 432 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 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 PetscCall(TSDIRKRegister(TSDIRK658A, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 435 } 436 { 437 // SDIRK(9,6)[1]SAL[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 438 const PetscReal A[9][9] = { 439 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 449 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 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 PetscCall(TSDIRKRegister(TSDIRKS659A, 6, 9, &A[0][0], b, NULL, bembed, 1, b)); 452 } 453 { 454 // DIRK(10,7)[1]SAL[(10,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 455 const PetscReal A[10][10] = { 456 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 467 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 const PetscReal bembed[10] = 469 {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 PetscCall(TSDIRKRegister(TSDIRK7510SAL, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 471 } 472 { 473 // ESDIRK(10,7)[2]SA[(10,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 474 const PetscReal A[10][10] = { 475 {RC(0.0), RC(0.0), 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 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 486 const PetscReal b[10] = {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), 487 RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)}; 488 const PetscReal bembed[10] = {RC(-0.015494246543626), RC(0.167657963820093), RC(0.269858958144236), RC(-0.0443258997755156), RC(0.150049236875266), 489 RC(0.259452082755846), RC(0.244624573502521), RC(-0.215528446920284), RC(0.0487601760292619), RC(0.134945602112201)}; 490 PetscCall(TSDIRKRegister(TSDIRKES7510SA, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 491 } 492 { 493 // DIRK(9,7)[1]A[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 494 const PetscReal A[9][9] = { 495 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 505 506 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 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 PetscCall(TSDIRKRegister(TSDIRK759A, 7, 9, &A[0][0], b, NULL, bembed, 1, b)); 509 } 510 { 511 // SDIRK(11,7)[1]SAL[(11,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 512 const PetscReal A[11][11] = { 513 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 525 const PetscReal b[11] = 526 {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 const PetscReal bembed[11] = 528 {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 PetscCall(TSDIRKRegister(TSDIRKS7511SAL, 7, 11, &A[0][0], b, NULL, bembed, 1, b)); 530 } 531 { 532 // DIRK(13,8)[1]A[(14,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 533 const PetscReal A[14][14] = { 534 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 549 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 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 RC(0.0417664613347638), RC(0.223636394275293), RC(0.231569156867596), RC(0.240526201277663), RC(-0.222933582911926), RC(-0.0111479879597561), RC(0.19915314335888)}; 552 PetscCall(TSDIRKRegister(TSDIRK8614A, 8, 14, &A[0][0], b, NULL, bembed, 1, b)); 553 } 554 { 555 // DIRK(15,8)[1]SAL[(16,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 556 const PetscReal A[16][16] = { 557 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 574 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 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 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 PetscCall(TSDIRKRegister(TSDIRK8616SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 578 } 579 { 580 // ESDIRK(16,8)[2]SAL[(16,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 581 const PetscReal A[16][16] = { 582 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), 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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 {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 }; 599 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 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 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 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 PetscCall(TSDIRKRegister(TSDIRKES8516SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 604 } 605 606 /* Additive methods */ 607 { 608 const PetscReal A[3][3] = { 609 {0.0, 0.0, 0.0}, 610 {0.0, 0.0, 0.0}, 611 {0.0, 0.5, 0.0} 612 }; 613 const PetscReal At[3][3] = { 614 {1.0, 0.0, 0.0}, 615 {0.0, 0.5, 0.0}, 616 {0.0, 0.5, 0.5} 617 }; 618 const PetscReal b[3] = {0.0, 0.5, 0.5}; 619 const PetscReal bembedt[3] = {1.0, 0.0, 0.0}; 620 PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 621 } 622 { 623 const PetscReal A[2][2] = { 624 {0.0, 0.0}, 625 {0.5, 0.0} 626 }; 627 const PetscReal At[2][2] = { 628 {0.0, 0.0}, 629 {0.0, 0.5} 630 }; 631 const PetscReal b[2] = {0.0, 1.0}; 632 const PetscReal bembedt[2] = {0.5, 0.5}; 633 /* 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 */ 634 PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 635 } 636 { 637 const PetscReal A[2][2] = { 638 {0.0, 0.0}, 639 {1.0, 0.0} 640 }; 641 const PetscReal At[2][2] = { 642 {0.0, 0.0}, 643 {0.5, 0.5} 644 }; 645 const PetscReal b[2] = {0.5, 0.5}; 646 const PetscReal bembedt[2] = {0.0, 1.0}; 647 /* 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 */ 648 PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 649 } 650 { 651 const PetscReal A[2][2] = { 652 {0.0, 0.0}, 653 {1.0, 0.0} 654 }; 655 const PetscReal At[2][2] = { 656 {us2, 0.0}, 657 {1.0 - 2.0 * us2, us2} 658 }; 659 const PetscReal b[2] = {0.5, 0.5}; 660 const PetscReal bembedt[2] = {0.0, 1.0}; 661 const PetscReal binterpt[2][2] = { 662 {(us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))}, 663 {1 - (us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))} 664 }; 665 const PetscReal binterp[2][2] = { 666 {1.0, -0.5}, 667 {0.0, 0.5 } 668 }; 669 PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 670 } 671 { 672 const PetscReal A[3][3] = { 673 {0, 0, 0}, 674 {2 - s2, 0, 0}, 675 {0.5, 0.5, 0} 676 }; 677 const PetscReal At[3][3] = { 678 {0, 0, 0 }, 679 {1 - 1 / s2, 1 - 1 / s2, 0 }, 680 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 681 }; 682 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 683 const PetscReal binterpt[3][2] = { 684 {1.0 / s2, -1.0 / (2.0 * s2)}, 685 {1.0 / s2, -1.0 / (2.0 * s2)}, 686 {1.0 - s2, 1.0 / s2 } 687 }; 688 PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 689 } 690 { 691 const PetscReal A[3][3] = { 692 {0, 0, 0}, 693 {2 - s2, 0, 0}, 694 {0.75, 0.25, 0} 695 }; 696 const PetscReal At[3][3] = { 697 {0, 0, 0 }, 698 {1 - 1 / s2, 1 - 1 / s2, 0 }, 699 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 700 }; 701 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 702 const PetscReal binterpt[3][2] = { 703 {1.0 / s2, -1.0 / (2.0 * s2)}, 704 {1.0 / s2, -1.0 / (2.0 * s2)}, 705 {1.0 - s2, 1.0 / s2 } 706 }; 707 PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 708 } 709 { /* Optimal for linear implicit part */ 710 const PetscReal A[3][3] = { 711 {0, 0, 0}, 712 {2 - s2, 0, 0}, 713 {(3 - 2 * s2) / 6, (3 + 2 * s2) / 6, 0} 714 }; 715 const PetscReal At[3][3] = { 716 {0, 0, 0 }, 717 {1 - 1 / s2, 1 - 1 / s2, 0 }, 718 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 719 }; 720 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 721 const PetscReal binterpt[3][2] = { 722 {1.0 / s2, -1.0 / (2.0 * s2)}, 723 {1.0 / s2, -1.0 / (2.0 * s2)}, 724 {1.0 - s2, 1.0 / s2 } 725 }; 726 PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 727 } 728 { /* Optimal for linear implicit part */ 729 const PetscReal A[3][3] = { 730 {0, 0, 0}, 731 {0.5, 0, 0}, 732 {0.5, 0.5, 0} 733 }; 734 const PetscReal At[3][3] = { 735 {0.25, 0, 0 }, 736 {0, 0.25, 0 }, 737 {1. / 3, 1. / 3, 1. / 3} 738 }; 739 PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 740 } 741 { 742 const PetscReal A[4][4] = { 743 {0, 0, 0, 0}, 744 {1767732205903. / 2027836641118., 0, 0, 0}, 745 {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 746 {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 747 }; 748 const PetscReal At[4][4] = { 749 {0, 0, 0, 0 }, 750 {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0 }, 751 {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0 }, 752 {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.} 753 }; 754 const PetscReal bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}; 755 const PetscReal binterpt[4][2] = { 756 {4655552711362. / 22874653954995., -215264564351. / 13552729205753. }, 757 {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119. }, 758 {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, 759 {584795268549. / 6622622206610., 2508943948391. / 7218656332882. } 760 }; 761 PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 762 } 763 { 764 const PetscReal A[5][5] = { 765 {0, 0, 0, 0, 0}, 766 {1. / 2, 0, 0, 0, 0}, 767 {11. / 18, 1. / 18, 0, 0, 0}, 768 {5. / 6, -5. / 6, .5, 0, 0}, 769 {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 770 }; 771 const PetscReal At[5][5] = { 772 {0, 0, 0, 0, 0 }, 773 {0, 1. / 2, 0, 0, 0 }, 774 {0, 1. / 6, 1. / 2, 0, 0 }, 775 {0, -1. / 2, 1. / 2, 1. / 2, 0 }, 776 {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2} 777 }; 778 PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 779 } 780 { 781 const PetscReal A[5][5] = { 782 {0, 0, 0, 0, 0}, 783 {1, 0, 0, 0, 0}, 784 {4. / 9, 2. / 9, 0, 0, 0}, 785 {1. / 4, 0, 3. / 4, 0, 0}, 786 {1. / 4, 0, 3. / 5, 0, 0} 787 }; 788 const PetscReal At[5][5] = { 789 {0, 0, 0, 0, 0 }, 790 {.5, .5, 0, 0, 0 }, 791 {5. / 18, -1. / 9, .5, 0, 0 }, 792 {.5, 0, 0, .5, 0 }, 793 {.25, 0, .75, -.5, .5} 794 }; 795 PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 796 } 797 { 798 const PetscReal A[6][6] = { 799 {0, 0, 0, 0, 0, 0}, 800 {1. / 2, 0, 0, 0, 0, 0}, 801 {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 802 {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 803 {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 804 {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 805 }; 806 const PetscReal At[6][6] = { 807 {0, 0, 0, 0, 0, 0 }, 808 {1. / 4, 1. / 4, 0, 0, 0, 0 }, 809 {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0 }, 810 {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0 }, 811 {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0 }, 812 {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4} 813 }; 814 const PetscReal bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}; 815 const PetscReal binterpt[6][3] = { 816 {6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025. }, 817 {0, 0, 0 }, 818 {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035. }, 819 {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 820 {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469. }, 821 {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.} 822 }; 823 PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 824 } 825 { 826 const PetscReal A[8][8] = { 827 {0, 0, 0, 0, 0, 0, 0, 0}, 828 {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 829 {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 830 {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 831 {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 832 {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 833 {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 834 {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 835 }; 836 const PetscReal At[8][8] = { 837 {0, 0, 0, 0, 0, 0, 0, 0 }, 838 {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0 }, 839 {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0 }, 840 {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0 }, 841 {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0 }, 842 {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0 }, 843 {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0 }, 844 {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.} 845 }; 846 const PetscReal bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}; 847 const PetscReal binterpt[8][3] = { 848 {-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439. }, 849 {0, 0, 0 }, 850 {0, 0, 0 }, 851 {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, 852 {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, 853 {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, 854 {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, 855 {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460. } 856 }; 857 PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 858 } 859 #undef RC 860 #undef us2 861 #undef s2 862 PetscFunctionReturn(PETSC_SUCCESS); 863 } 864 865 /*@C 866 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`. 867 868 Not Collective 869 870 Level: advanced 871 872 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 873 @*/ 874 PetscErrorCode TSARKIMEXRegisterDestroy(void) 875 { 876 ARKTableauLink link; 877 878 PetscFunctionBegin; 879 while ((link = ARKTableauList)) { 880 ARKTableau t = &link->tab; 881 ARKTableauList = link->next; 882 PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 883 PetscCall(PetscFree2(t->bembedt, t->bembed)); 884 PetscCall(PetscFree2(t->binterpt, t->binterp)); 885 PetscCall(PetscFree(t->name)); 886 PetscCall(PetscFree(link)); 887 } 888 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 889 PetscFunctionReturn(PETSC_SUCCESS); 890 } 891 892 /*@C 893 TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called 894 from `TSInitializePackage()`. 895 896 Level: developer 897 898 .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()` 899 @*/ 900 PetscErrorCode TSARKIMEXInitializePackage(void) 901 { 902 PetscFunctionBegin; 903 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 904 TSARKIMEXPackageInitialized = PETSC_TRUE; 905 PetscCall(TSARKIMEXRegisterAll()); 906 PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 907 PetscFunctionReturn(PETSC_SUCCESS); 908 } 909 910 /*@C 911 TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is 912 called from `PetscFinalize()`. 913 914 Level: developer 915 916 .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()` 917 @*/ 918 PetscErrorCode TSARKIMEXFinalizePackage(void) 919 { 920 PetscFunctionBegin; 921 TSARKIMEXPackageInitialized = PETSC_FALSE; 922 PetscCall(TSARKIMEXRegisterDestroy()); 923 PetscFunctionReturn(PETSC_SUCCESS); 924 } 925 926 /*@C 927 TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 928 929 Logically Collective. 930 931 Input Parameters: 932 + name - identifier for method 933 . order - approximation order of method 934 . s - number of stages, this is the dimension of the matrices below 935 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 936 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 937 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 938 . A - Non-stiff stage coefficients (dimension s*s, row-major) 939 . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 940 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 941 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 942 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 943 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 944 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 945 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 946 947 Level: advanced 948 949 Note: 950 Several `TSARKIMEX` methods are provided, this function is only needed to create new methods. 951 952 .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS` 953 @*/ 954 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[]) 955 { 956 ARKTableauLink link; 957 ARKTableau t; 958 PetscInt i, j; 959 960 PetscFunctionBegin; 961 PetscCall(TSARKIMEXInitializePackage()); 962 for (link = ARKTableauList; link; link = link->next) { 963 PetscBool match; 964 965 PetscCall(PetscStrcmp(link->tab.name, name, &match)); 966 PetscCheck(!match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Method with name \"%s\" already registered", name); 967 } 968 PetscCall(PetscNew(&link)); 969 t = &link->tab; 970 PetscCall(PetscStrallocpy(name, &t->name)); 971 t->order = order; 972 t->s = s; 973 PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 974 PetscCall(PetscArraycpy(t->At, At, s * s)); 975 if (A) { 976 PetscCall(PetscArraycpy(t->A, A, s * s)); 977 t->additive = PETSC_TRUE; 978 } 979 980 if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 981 else 982 for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 983 984 if (t->additive) { 985 if (b) PetscCall(PetscArraycpy(t->b, b, s)); 986 else 987 for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 988 } else PetscCall(PetscArrayzero(t->b, s)); 989 990 if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 991 else 992 for (i = 0; i < s; i++) 993 for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 994 995 if (t->additive) { 996 if (c) PetscCall(PetscArraycpy(t->c, c, s)); 997 else 998 for (i = 0; i < s; i++) 999 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 1000 } else PetscCall(PetscArrayzero(t->c, s)); 1001 1002 t->stiffly_accurate = PETSC_TRUE; 1003 for (i = 0; i < s; i++) 1004 if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 1005 1006 t->explicit_first_stage = PETSC_TRUE; 1007 for (i = 0; i < s; i++) 1008 if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 1009 1010 /* def of FSAL can be made more precise */ 1011 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 1012 1013 if (bembedt) { 1014 PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 1015 PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 1016 PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 1017 } 1018 1019 t->pinterp = pinterp; 1020 PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 1021 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 1022 PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 1023 1024 link->next = ARKTableauList; 1025 ARKTableauList = link; 1026 PetscFunctionReturn(PETSC_SUCCESS); 1027 } 1028 1029 /*@C 1030 TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation 1031 1032 Logically Collective. 1033 1034 Input Parameters: 1035 + name - identifier for method 1036 . order - approximation order of method 1037 . s - number of stages, this is the dimension of the matrices below 1038 . At - Butcher table of stage coefficients (dimension `s`*`s`, row-major order) 1039 . bt - Butcher table for completing the step (dimension `s`; pass `NULL` to use the last row of `At`) 1040 . ct - Abscissa of each stage (dimension s, NULL to use row sums of At) 1041 . bembedt - Stiff part of completion table for embedded method (dimension s; `NULL` if not available) 1042 . pinterp - Order of the interpolation scheme, equal to the number of columns of `binterpt` and `binterp` 1043 - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 1044 1045 Level: advanced 1046 1047 Note: 1048 Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods. 1049 1050 .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS` 1051 @*/ 1052 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 { 1054 PetscFunctionBegin; 1055 PetscCall(TSARKIMEXRegister(name, order, s, At, bt, ct, NULL, NULL, NULL, bembedt, NULL, pinterp, binterpt, NULL)); 1056 PetscFunctionReturn(PETSC_SUCCESS); 1057 } 1058 1059 /* 1060 The step completion formula is 1061 1062 x1 = x0 - h bt^T YdotI + h b^T YdotRHS 1063 1064 This function can be called before or after ts->vec_sol has been updated. 1065 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 1066 We can write 1067 1068 x1e = x0 - h bet^T YdotI + h be^T YdotRHS 1069 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 1070 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 1071 1072 so we can evaluate the method with different order even after the step has been optimistically completed. 1073 */ 1074 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) 1075 { 1076 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1077 ARKTableau tab = ark->tableau; 1078 PetscScalar *w = ark->work; 1079 PetscReal h; 1080 PetscInt s = tab->s, j; 1081 PetscBool hasE; 1082 1083 PetscFunctionBegin; 1084 switch (ark->status) { 1085 case TS_STEP_INCOMPLETE: 1086 case TS_STEP_PENDING: 1087 h = ts->time_step; 1088 break; 1089 case TS_STEP_COMPLETE: 1090 h = ts->ptime - ts->ptime_prev; 1091 break; 1092 default: 1093 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1094 } 1095 if (order == tab->order) { 1096 if (ark->status == TS_STEP_INCOMPLETE) { 1097 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 1098 PetscCall(VecCopy(ark->Y[s - 1], X)); 1099 } else { /* Use the standard completion formula (bt,b) */ 1100 PetscCall(VecCopy(ts->vec_sol, X)); 1101 for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 1102 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1103 if (tab->additive && ark->imex) { /* Method is IMEX, complete the explicit formula */ 1104 PetscCall(TSHasRHSFunction(ts, &hasE)); 1105 if (hasE) { 1106 for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 1107 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1108 } 1109 } 1110 } 1111 } else PetscCall(VecCopy(ts->vec_sol, X)); 1112 if (done) *done = PETSC_TRUE; 1113 PetscFunctionReturn(PETSC_SUCCESS); 1114 } else if (order == tab->order - 1) { 1115 if (!tab->bembedt) goto unavailable; 1116 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 1117 PetscCall(VecCopy(ts->vec_sol, X)); 1118 for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 1119 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1120 if (tab->additive) { 1121 PetscCall(TSHasRHSFunction(ts, &hasE)); 1122 if (hasE) { 1123 for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 1124 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1125 } 1126 } 1127 } else { /* Rollback and re-complete using (bet-be,be-b) */ 1128 PetscCall(VecCopy(ts->vec_sol, X)); 1129 for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 1130 PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 1131 if (tab->additive) { 1132 PetscCall(TSHasRHSFunction(ts, &hasE)); 1133 if (hasE) { 1134 for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 1135 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1136 } 1137 } 1138 } 1139 if (done) *done = PETSC_TRUE; 1140 PetscFunctionReturn(PETSC_SUCCESS); 1141 } 1142 unavailable: 1143 if (done) *done = PETSC_FALSE; 1144 else 1145 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, 1146 tab->order, order); 1147 PetscFunctionReturn(PETSC_SUCCESS); 1148 } 1149 1150 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 1151 { 1152 Vec Udot, Y1, Y2; 1153 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1154 PetscReal norm; 1155 1156 PetscFunctionBegin; 1157 PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 1158 PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 1159 PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 1160 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 1161 PetscCall(VecSetRandom(Udot, NULL)); 1162 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 1163 PetscCall(VecAXPY(Y2, -1.0, Y1)); 1164 PetscCall(VecAXPY(Y2, -1.0, Udot)); 1165 PetscCall(VecNorm(Y2, NORM_2, &norm)); 1166 if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 1167 *id = PETSC_TRUE; 1168 } else { 1169 *id = PETSC_FALSE; 1170 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)); 1171 } 1172 PetscCall(VecDestroy(&Udot)); 1173 PetscCall(VecDestroy(&Y1)); 1174 PetscCall(VecDestroy(&Y2)); 1175 PetscFunctionReturn(PETSC_SUCCESS); 1176 } 1177 1178 static PetscErrorCode TSRollBack_ARKIMEX(TS ts) 1179 { 1180 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1181 ARKTableau tab = ark->tableau; 1182 const PetscInt s = tab->s; 1183 const PetscReal *bt = tab->bt, *b = tab->b; 1184 PetscScalar *w = ark->work; 1185 Vec *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS; 1186 PetscInt j; 1187 PetscReal h; 1188 1189 PetscFunctionBegin; 1190 switch (ark->status) { 1191 case TS_STEP_INCOMPLETE: 1192 case TS_STEP_PENDING: 1193 h = ts->time_step; 1194 break; 1195 case TS_STEP_COMPLETE: 1196 h = ts->ptime - ts->ptime_prev; 1197 break; 1198 default: 1199 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1200 } 1201 for (j = 0; j < s; j++) w[j] = -h * bt[j]; 1202 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotI)); 1203 if (tab->additive) { 1204 PetscBool hasE; 1205 1206 PetscCall(TSHasRHSFunction(ts, &hasE)); 1207 if (hasE) { 1208 for (j = 0; j < s; j++) w[j] = -h * b[j]; 1209 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS)); 1210 } 1211 } 1212 PetscFunctionReturn(PETSC_SUCCESS); 1213 } 1214 1215 static PetscErrorCode TSStep_ARKIMEX(TS ts) 1216 { 1217 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1218 ARKTableau tab = ark->tableau; 1219 const PetscInt s = tab->s; 1220 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 1221 PetscScalar *w = ark->work; 1222 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 1223 PetscBool extrapolate = ark->extrapolate; 1224 TSAdapt adapt; 1225 SNES snes; 1226 PetscInt i, j, its, lits; 1227 PetscInt rejections = 0; 1228 PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE; 1229 PetscReal next_time_step = ts->time_step; 1230 1231 PetscFunctionBegin; 1232 if (ark->extrapolate && !ark->Y_prev) { 1233 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1234 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1235 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1236 } 1237 1238 if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE)); 1239 if (!hasE) dirk = PETSC_TRUE; 1240 1241 if (!ts->steprollback) { 1242 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 1243 PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 1244 } 1245 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 1246 for (i = 0; i < s; i++) { 1247 PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 1248 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 1249 if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 1250 } 1251 } 1252 } 1253 1254 /* For fully implicit formulations we can solve the equations 1255 F(tn,xn,xdot) = 0 1256 for the explicit first stage */ 1257 if (dirk && tab->explicit_first_stage && ts->steprestart) { 1258 ark->scoeff = 0.0; 1259 PetscCall(VecCopy(ts->vec_sol, Z)); 1260 PetscCall(TSGetSNES(ts, &snes)); 1261 PetscCall(SNESSolve(snes, NULL, Ydot0)); 1262 } 1263 1264 /* For IMEX we compute a step */ 1265 if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 1266 TS ts_start; 1267 if (PetscDefined(USE_DEBUG) && hasE) { 1268 PetscBool id = PETSC_FALSE; 1269 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1270 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"); 1271 } 1272 PetscCall(TSClone(ts, &ts_start)); 1273 PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 1274 PetscCall(TSSetTime(ts_start, ts->ptime)); 1275 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 1276 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 1277 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 1278 PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 1279 PetscCall(TSSetType(ts_start, TSARKIMEX)); 1280 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 1281 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 1282 1283 PetscCall(TSRestartStep(ts_start)); 1284 PetscCall(TSSolve(ts_start, ts->vec_sol)); 1285 PetscCall(TSGetTime(ts_start, &ts->ptime)); 1286 PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 1287 1288 { /* Save the initial slope for the next step */ 1289 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 1290 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 1291 } 1292 ts->steps++; 1293 1294 /* Set the correct TS in SNES */ 1295 /* We'll try to bypass this by changing the method on the fly */ 1296 { 1297 PetscCall(TSGetSNES(ts, &snes)); 1298 PetscCall(TSSetSNES(ts, snes)); 1299 } 1300 PetscCall(TSDestroy(&ts_start)); 1301 } 1302 1303 ark->status = TS_STEP_INCOMPLETE; 1304 while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 1305 PetscReal t = ts->ptime; 1306 PetscReal h = ts->time_step; 1307 for (i = 0; i < s; i++) { 1308 ark->stage_time = t + h * ct[i]; 1309 PetscCall(TSPreStage(ts, ark->stage_time)); 1310 if (At[i * s + i] == 0) { /* This stage is explicit */ 1311 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"); 1312 PetscCall(VecCopy(ts->vec_sol, Y[i])); 1313 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1314 PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 1315 if (tab->additive && hasE) { 1316 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1317 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 1318 } 1319 } else { 1320 ark->scoeff = 1. / At[i * s + i]; 1321 /* Ydot = shift*(Y-Z) */ 1322 PetscCall(VecCopy(ts->vec_sol, Z)); 1323 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1324 PetscCall(VecMAXPY(Z, i, w, YdotI)); 1325 if (tab->additive && hasE) { 1326 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1327 PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 1328 } 1329 if (extrapolate && !ts->steprestart) { 1330 /* Initial guess extrapolated from previous time step stage values */ 1331 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 1332 } else { 1333 /* Initial guess taken from last stage */ 1334 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 1335 } 1336 PetscCall(TSGetSNES(ts, &snes)); 1337 PetscCall(SNESSolve(snes, NULL, Y[i])); 1338 PetscCall(SNESGetIterationNumber(snes, &its)); 1339 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 1340 ts->snes_its += its; 1341 ts->ksp_its += lits; 1342 PetscCall(TSGetAdapt(ts, &adapt)); 1343 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 1344 if (!stageok) { 1345 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 1346 * use extrapolation to initialize the solves on the next attempt. */ 1347 extrapolate = PETSC_FALSE; 1348 goto reject_step; 1349 } 1350 } 1351 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { 1352 if (i == 0 && tab->explicit_first_stage) { 1353 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 ((PetscObject)ts)->type_name, ark->tableau->name); 1355 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 1356 } else { 1357 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1358 } 1359 } else { 1360 if (i == 0 && tab->explicit_first_stage) { 1361 PetscCall(VecZeroEntries(Ydot)); 1362 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 1363 PetscCall(VecScale(YdotI[i], -1.0)); 1364 } else { 1365 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1366 } 1367 if (hasE) { 1368 if (ark->imex) { 1369 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 1370 } else { 1371 PetscCall(VecZeroEntries(YdotRHS[i])); 1372 } 1373 } 1374 } 1375 PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 1376 } 1377 1378 ark->status = TS_STEP_INCOMPLETE; 1379 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 1380 ark->status = TS_STEP_PENDING; 1381 PetscCall(TSGetAdapt(ts, &adapt)); 1382 PetscCall(TSAdaptCandidatesClear(adapt)); 1383 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 1384 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1385 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1386 if (!accept) { /* Roll back the current step */ 1387 PetscCall(TSRollBack_ARKIMEX(ts)); 1388 ts->time_step = next_time_step; 1389 goto reject_step; 1390 } 1391 1392 ts->ptime += ts->time_step; 1393 ts->time_step = next_time_step; 1394 break; 1395 1396 reject_step: 1397 ts->reject++; 1398 accept = PETSC_FALSE; 1399 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1400 ts->reason = TS_DIVERGED_STEP_REJECTED; 1401 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1402 } 1403 } 1404 PetscFunctionReturn(PETSC_SUCCESS); 1405 } 1406 1407 /* 1408 This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 1409 Udot = H(t,U) + G(t,U) 1410 This corresponds to F(t,U,Udot) = Udot-H(t,U). 1411 1412 The complete adjoint equations are 1413 (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 1414 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 1415 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0 1416 lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 1417 mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 1418 + (b_i dGdP + bt[i] dHdP) lambda_{n+1} + dGdP \sum_{j=i+1}^s a[j][i] lambda_s[j] 1419 + dHdP \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0 1420 mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 1421 where shift = 1/(at[i][i]*h) 1422 1423 If at[i][i] is 0, the first equation falls back to 1424 lambda_s[i] = h ( 1425 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 1426 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 1427 1428 */ 1429 static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 1430 { 1431 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1432 ARKTableau tab = ark->tableau; 1433 const PetscInt s = tab->s; 1434 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 1435 PetscScalar *w = ark->work; 1436 Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 1437 Mat Jex, Jim, Jimpre; 1438 PetscInt i, j, nadj; 1439 PetscReal t = ts->ptime, stage_time_ex; 1440 PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 1441 KSP ksp; 1442 1443 PetscFunctionBegin; 1444 ark->status = TS_STEP_INCOMPLETE; 1445 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1446 PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 1447 PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 1448 1449 for (i = s - 1; i >= 0; i--) { 1450 ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 1451 stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 1452 if (At[i * s + i] == 0) { // This stage is explicit 1453 ark->scoeff = 0.; 1454 } else { 1455 ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 1456 } 1457 PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 1458 PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 1459 if (ts->vecs_sensip) { 1460 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 1461 PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 1462 } 1463 for (nadj = 0; nadj < ts->numcost; nadj++) { 1464 /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 1465 if (s - i - 1 > 0) { 1466 /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 1467 for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 1468 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1469 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1470 /* (shift I - dHdU) Temp */ 1471 PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 1472 /* cancel out shift Temp where shift=-scoeff/h */ 1473 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 1474 if (ts->vecs_sensip) { 1475 /* - dHdP Temp */ 1476 PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1477 /* mu_n += h dHdP Temp */ 1478 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsSensiPTemp[nadj])); 1479 } 1480 /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 1481 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 1482 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1483 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1484 /* dGdU Temp */ 1485 PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1486 if (ts->vecs_sensip) { 1487 /* dGdP Temp */ 1488 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1489 /* mu_n += h dGdP Temp */ 1490 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1491 } 1492 } else { 1493 PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 1494 } 1495 if (bt[i]) { // bt[i] dHdU lambda_{n+1} 1496 /* (shift I - dHdU)^T lambda_{n+1} */ 1497 PetscCall(MatMultTranspose(Jim, ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); 1498 /* -bt[i] (shift I - dHdU)^T lambda_{n+1} */ 1499 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i], VecsSensiTemp[nadj])); 1500 /* cancel out -bt[i] shift lambda_{n+1} */ 1501 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i] * ark->scoeff / adjoint_time_step, ts->vecs_sensi[nadj])); 1502 if (ts->vecs_sensip) { 1503 /* -dHdP lambda_{n+1} */ 1504 PetscCall(MatMultTranspose(ts->Jacp, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj])); 1505 /* mu_n += h bt[i] dHdP lambda_{n+1} */ 1506 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -bt[i] * adjoint_time_step, VecsSensiPTemp[nadj])); 1507 } 1508 } 1509 if (b[i]) { // b[i] dGdU lambda_{n+1} 1510 /* dGdU lambda_{n+1} */ 1511 PetscCall(MatMultTranspose(Jex, ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); 1512 /* b[i] dGdU lambda_{n+1} */ 1513 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], b[i], VecsSensiTemp[nadj])); 1514 if (ts->vecs_sensip) { 1515 /* dGdP lambda_{n+1} */ 1516 PetscCall(MatMultTranspose(ts->Jacprhs, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj])); 1517 /* mu_n += h b[i] dGdP lambda_{n+1} */ 1518 PetscCall(VecAXPY(ts->vecs_sensip[nadj], b[i] * adjoint_time_step, VecsSensiPTemp[nadj])); 1519 } 1520 } 1521 /* Build LHS for first-order adjoint */ 1522 if (At[i * s + i] == 0) { // This stage is explicit 1523 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 1524 } else { 1525 KSP ksp; 1526 KSPConvergedReason kspreason; 1527 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1528 PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 1529 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 1530 PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1531 PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 1532 if (kspreason < 0) { 1533 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 1534 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj)); 1535 } 1536 if (ts->vecs_sensip) { 1537 /* -dHdP lambda_s[i] */ 1538 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 1539 /* mu_n += h at[i][i] dHdP lambda_s[i] */ 1540 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 1541 } 1542 } 1543 } 1544 } 1545 for (j = 0; j < s; j++) w[j] = 1.0; 1546 for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 1547 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 1548 ark->status = TS_STEP_COMPLETE; 1549 PetscFunctionReturn(PETSC_SUCCESS); 1550 } 1551 1552 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1553 { 1554 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1555 ARKTableau tab = ark->tableau; 1556 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1557 PetscReal h; 1558 PetscReal tt, t; 1559 PetscScalar *bt = ark->work, *b = ark->work + s; 1560 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1561 1562 PetscFunctionBegin; 1563 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name); 1564 switch (ark->status) { 1565 case TS_STEP_INCOMPLETE: 1566 case TS_STEP_PENDING: 1567 h = ts->time_step; 1568 t = (itime - ts->ptime) / h; 1569 break; 1570 case TS_STEP_COMPLETE: 1571 h = ts->ptime - ts->ptime_prev; 1572 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1573 break; 1574 default: 1575 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1576 } 1577 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1578 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1579 for (i = 0; i < s; i++) { 1580 bt[i] += h * Bt[i * pinterp + j] * tt; 1581 b[i] += h * B[i * pinterp + j] * tt; 1582 } 1583 } 1584 PetscCall(VecCopy(ark->Y[0], X)); 1585 PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1586 if (tab->additive) { 1587 PetscBool hasE; 1588 PetscCall(TSHasRHSFunction(ts, &hasE)); 1589 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1590 } 1591 PetscFunctionReturn(PETSC_SUCCESS); 1592 } 1593 1594 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1595 { 1596 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1597 ARKTableau tab = ark->tableau; 1598 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1599 PetscReal h, h_prev, t, tt; 1600 PetscScalar *bt = ark->work, *b = ark->work + s; 1601 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1602 1603 PetscFunctionBegin; 1604 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 1605 h = ts->time_step; 1606 h_prev = ts->ptime - ts->ptime_prev; 1607 t = 1 + h / h_prev * c; 1608 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1609 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1610 for (i = 0; i < s; i++) { 1611 bt[i] += h * Bt[i * pinterp + j] * tt; 1612 b[i] += h * B[i * pinterp + j] * tt; 1613 } 1614 } 1615 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 1616 PetscCall(VecCopy(ark->Y_prev[0], X)); 1617 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1618 if (tab->additive) { 1619 PetscBool hasE; 1620 PetscCall(TSHasRHSFunction(ts, &hasE)); 1621 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1622 } 1623 PetscFunctionReturn(PETSC_SUCCESS); 1624 } 1625 1626 static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1627 { 1628 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1629 ARKTableau tab = ark->tableau; 1630 1631 PetscFunctionBegin; 1632 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1633 PetscCall(PetscFree(ark->work)); 1634 PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 1635 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 1636 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 1637 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 1638 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 1639 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 1640 PetscFunctionReturn(PETSC_SUCCESS); 1641 } 1642 1643 static PetscErrorCode TSReset_ARKIMEX(TS ts) 1644 { 1645 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1646 1647 PetscFunctionBegin; 1648 PetscCall(TSARKIMEXTableauReset(ts)); 1649 PetscCall(VecDestroy(&ark->Ydot)); 1650 PetscCall(VecDestroy(&ark->Ydot0)); 1651 PetscCall(VecDestroy(&ark->Z)); 1652 PetscFunctionReturn(PETSC_SUCCESS); 1653 } 1654 1655 static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 1656 { 1657 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1658 ARKTableau tab = ark->tableau; 1659 1660 PetscFunctionBegin; 1661 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 1662 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 1663 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 1664 PetscFunctionReturn(PETSC_SUCCESS); 1665 } 1666 1667 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1668 { 1669 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1670 1671 PetscFunctionBegin; 1672 if (Z) { 1673 if (dm && dm != ts->dm) { 1674 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1675 } else *Z = ax->Z; 1676 } 1677 if (Ydot) { 1678 if (dm && dm != ts->dm) { 1679 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1680 } else *Ydot = ax->Ydot; 1681 } 1682 PetscFunctionReturn(PETSC_SUCCESS); 1683 } 1684 1685 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1686 { 1687 PetscFunctionBegin; 1688 if (Z) { 1689 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1690 } 1691 if (Ydot) { 1692 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1693 } 1694 PetscFunctionReturn(PETSC_SUCCESS); 1695 } 1696 1697 /* 1698 This defines the nonlinear equation that is to be solved with SNES 1699 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1700 */ 1701 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1702 { 1703 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1704 DM dm, dmsave; 1705 Vec Z, Ydot; 1706 PetscReal shift = ark->scoeff / ts->time_step; 1707 1708 PetscFunctionBegin; 1709 PetscCall(SNESGetDM(snes, &dm)); 1710 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1711 dmsave = ts->dm; 1712 ts->dm = dm; 1713 1714 if (ark->scoeff == 0.0) { 1715 /* We are solving F(t,x_n,xdot) = 0 to start the method */ 1716 PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 1717 } else { 1718 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1719 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1720 } 1721 1722 ts->dm = dmsave; 1723 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1724 PetscFunctionReturn(PETSC_SUCCESS); 1725 } 1726 1727 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1728 { 1729 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1730 DM dm, dmsave; 1731 Vec Ydot, Z; 1732 PetscReal shift = ark->scoeff / ts->time_step; 1733 1734 PetscFunctionBegin; 1735 PetscCall(SNESGetDM(snes, &dm)); 1736 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1737 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1738 dmsave = ts->dm; 1739 ts->dm = dm; 1740 1741 if (ark->scoeff == 0.0) { 1742 /* We are solving F(t,x_n,xdot) = 0 to start the method, we only only dF/dXdot 1743 Jed's proposal is to compute with a very large shift and scale back the matrix */ 1744 shift = 1.0 / PETSC_MACHINE_EPSILON; 1745 PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1746 PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 1747 if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1748 } else { 1749 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1750 } 1751 ts->dm = dmsave; 1752 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1753 PetscFunctionReturn(PETSC_SUCCESS); 1754 } 1755 1756 static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 1757 { 1758 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1759 1760 PetscFunctionBegin; 1761 if (ns) *ns = ark->tableau->s; 1762 if (Y) *Y = ark->Y; 1763 PetscFunctionReturn(PETSC_SUCCESS); 1764 } 1765 1766 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1767 { 1768 PetscFunctionBegin; 1769 PetscFunctionReturn(PETSC_SUCCESS); 1770 } 1771 1772 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1773 { 1774 TS ts = (TS)ctx; 1775 Vec Z, Z_c; 1776 1777 PetscFunctionBegin; 1778 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 1779 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 1780 PetscCall(MatRestrict(restrct, Z, Z_c)); 1781 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 1782 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 1783 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 1784 PetscFunctionReturn(PETSC_SUCCESS); 1785 } 1786 1787 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1788 { 1789 PetscFunctionBegin; 1790 PetscFunctionReturn(PETSC_SUCCESS); 1791 } 1792 1793 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1794 { 1795 TS ts = (TS)ctx; 1796 Vec Z, Z_c; 1797 1798 PetscFunctionBegin; 1799 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 1800 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 1801 1802 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1803 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1804 1805 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 1806 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 1807 PetscFunctionReturn(PETSC_SUCCESS); 1808 } 1809 1810 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 1811 { 1812 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1813 ARKTableau tab = ark->tableau; 1814 1815 PetscFunctionBegin; 1816 PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 1817 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 1818 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 1819 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 1820 if (ark->extrapolate) { 1821 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1822 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1823 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1824 } 1825 PetscFunctionReturn(PETSC_SUCCESS); 1826 } 1827 1828 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 1829 { 1830 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1831 DM dm; 1832 SNES snes; 1833 1834 PetscFunctionBegin; 1835 PetscCall(TSARKIMEXTableauSetUp(ts)); 1836 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 1837 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 1838 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 1839 PetscCall(TSGetDM(ts, &dm)); 1840 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 1841 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 1842 PetscCall(TSGetSNES(ts, &snes)); 1843 PetscFunctionReturn(PETSC_SUCCESS); 1844 } 1845 1846 static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 1847 { 1848 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1849 ARKTableau tab = ark->tableau; 1850 1851 PetscFunctionBegin; 1852 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 1853 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 1854 if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); } 1855 if (PetscDefined(USE_DEBUG)) { 1856 PetscBool id = PETSC_FALSE; 1857 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1858 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"); 1859 } 1860 PetscFunctionReturn(PETSC_SUCCESS); 1861 } 1862 1863 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject) 1864 { 1865 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1866 PetscBool dirk; 1867 1868 PetscFunctionBegin; 1869 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 1870 PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 1871 { 1872 ARKTableauLink link; 1873 PetscInt count, choice; 1874 PetscBool flg; 1875 const char **namelist; 1876 for (link = ARKTableauList, count = 0; link; link = link->next) { 1877 if (!dirk && link->tab.additive) count++; 1878 if (dirk && !link->tab.additive) count++; 1879 } 1880 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1881 for (link = ARKTableauList, count = 0; link; link = link->next) { 1882 if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 1883 if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 1884 } 1885 if (dirk) { 1886 PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 1887 if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 1888 } else { 1889 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 1890 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 1891 flg = (PetscBool)!ark->imex; 1892 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 1893 ark->imex = (PetscBool)!flg; 1894 } 1895 PetscCall(PetscFree(namelist)); 1896 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)); 1897 } 1898 PetscOptionsHeadEnd(); 1899 PetscFunctionReturn(PETSC_SUCCESS); 1900 } 1901 1902 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 1903 { 1904 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1905 PetscBool iascii, dirk; 1906 1907 PetscFunctionBegin; 1908 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 1909 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1910 if (iascii) { 1911 PetscViewerFormat format; 1912 ARKTableau tab = ark->tableau; 1913 TSARKIMEXType arktype; 1914 char buf[2048]; 1915 PetscBool flg; 1916 1917 PetscCall(TSARKIMEXGetType(ts, &arktype)); 1918 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 1919 PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 1920 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 1921 PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 1922 PetscCall(PetscViewerGetFormat(viewer, &format)); 1923 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1924 PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 1925 for (PetscInt i = 0; i < tab->s; i++) { 1926 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 1927 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 1928 } 1929 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 1930 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 1931 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 1932 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 1933 } 1934 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 1935 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 1936 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 1937 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 1938 if (!dirk) { 1939 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 1940 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 1941 } 1942 } 1943 PetscFunctionReturn(PETSC_SUCCESS); 1944 } 1945 1946 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 1947 { 1948 SNES snes; 1949 TSAdapt adapt; 1950 1951 PetscFunctionBegin; 1952 PetscCall(TSGetAdapt(ts, &adapt)); 1953 PetscCall(TSAdaptLoad(adapt, viewer)); 1954 PetscCall(TSGetSNES(ts, &snes)); 1955 PetscCall(SNESLoad(snes, viewer)); 1956 /* function and Jacobian context for SNES when used with TS is always ts object */ 1957 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 1958 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 1959 PetscFunctionReturn(PETSC_SUCCESS); 1960 } 1961 1962 /*@C 1963 TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 1964 1965 Logically Collective 1966 1967 Input Parameters: 1968 + ts - timestepping context 1969 - arktype - type of `TSARKIMEX` scheme 1970 1971 Options Database Key: 1972 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 1973 1974 Level: intermediate 1975 1976 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 1977 `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 1978 @*/ 1979 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 1980 { 1981 PetscFunctionBegin; 1982 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1983 PetscAssertPointer(arktype, 2); 1984 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 1985 PetscFunctionReturn(PETSC_SUCCESS); 1986 } 1987 1988 /*@C 1989 TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 1990 1991 Logically Collective 1992 1993 Input Parameter: 1994 . ts - timestepping context 1995 1996 Output Parameter: 1997 . arktype - type of `TSARKIMEX` scheme 1998 1999 Level: intermediate 2000 2001 .seealso: [](ch_ts), `TSARKIMEXc` 2002 @*/ 2003 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2004 { 2005 PetscFunctionBegin; 2006 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2007 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 2008 PetscFunctionReturn(PETSC_SUCCESS); 2009 } 2010 2011 /*@ 2012 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 2013 2014 Logically Collective 2015 2016 Input Parameters: 2017 + ts - timestepping context 2018 - flg - `PETSC_TRUE` for fully implicit 2019 2020 Level: intermediate 2021 2022 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 2023 @*/ 2024 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 2025 { 2026 PetscFunctionBegin; 2027 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2028 PetscValidLogicalCollectiveBool(ts, flg, 2); 2029 PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 2030 PetscFunctionReturn(PETSC_SUCCESS); 2031 } 2032 2033 /*@ 2034 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 2035 2036 Logically Collective 2037 2038 Input Parameter: 2039 . ts - timestepping context 2040 2041 Output Parameter: 2042 . flg - `PETSC_TRUE` for fully implicit 2043 2044 Level: intermediate 2045 2046 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 2047 @*/ 2048 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 2049 { 2050 PetscFunctionBegin; 2051 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2052 PetscAssertPointer(flg, 2); 2053 PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 2054 PetscFunctionReturn(PETSC_SUCCESS); 2055 } 2056 2057 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 2058 { 2059 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2060 2061 PetscFunctionBegin; 2062 *arktype = ark->tableau->name; 2063 PetscFunctionReturn(PETSC_SUCCESS); 2064 } 2065 2066 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 2067 { 2068 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2069 PetscBool match; 2070 ARKTableauLink link; 2071 2072 PetscFunctionBegin; 2073 if (ark->tableau) { 2074 PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 2075 if (match) PetscFunctionReturn(PETSC_SUCCESS); 2076 } 2077 for (link = ARKTableauList; link; link = link->next) { 2078 PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 2079 if (match) { 2080 if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 2081 ark->tableau = &link->tab; 2082 if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 2083 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 2084 PetscFunctionReturn(PETSC_SUCCESS); 2085 } 2086 } 2087 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 2088 } 2089 2090 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 2091 { 2092 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2093 2094 PetscFunctionBegin; 2095 ark->imex = (PetscBool)!flg; 2096 PetscFunctionReturn(PETSC_SUCCESS); 2097 } 2098 2099 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 2100 { 2101 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2102 2103 PetscFunctionBegin; 2104 *flg = (PetscBool)!ark->imex; 2105 PetscFunctionReturn(PETSC_SUCCESS); 2106 } 2107 2108 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 2109 { 2110 PetscFunctionBegin; 2111 PetscCall(TSReset_ARKIMEX(ts)); 2112 if (ts->dm) { 2113 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 2114 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2115 } 2116 PetscCall(PetscFree(ts->data)); 2117 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL)); 2118 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL)); 2119 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 2120 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 2121 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 2122 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 2123 PetscFunctionReturn(PETSC_SUCCESS); 2124 } 2125 2126 /* ------------------------------------------------------------ */ 2127 /*MC 2128 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 2129 2130 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2131 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2132 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2133 2134 Level: beginner 2135 2136 Notes: 2137 The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 2138 2139 If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2140 2141 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). 2142 2143 Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2144 2145 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2146 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2147 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 2148 M*/ 2149 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2150 { 2151 TS_ARKIMEX *ark; 2152 PetscBool dirk; 2153 2154 PetscFunctionBegin; 2155 PetscCall(TSARKIMEXInitializePackage()); 2156 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2157 2158 ts->ops->reset = TSReset_ARKIMEX; 2159 ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 2160 ts->ops->destroy = TSDestroy_ARKIMEX; 2161 ts->ops->view = TSView_ARKIMEX; 2162 ts->ops->load = TSLoad_ARKIMEX; 2163 ts->ops->setup = TSSetUp_ARKIMEX; 2164 ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 2165 ts->ops->step = TSStep_ARKIMEX; 2166 ts->ops->interpolate = TSInterpolate_ARKIMEX; 2167 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 2168 ts->ops->rollback = TSRollBack_ARKIMEX; 2169 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 2170 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 2171 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 2172 ts->ops->getstages = TSGetStages_ARKIMEX; 2173 ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 2174 2175 ts->usessnes = PETSC_TRUE; 2176 2177 PetscCall(PetscNew(&ark)); 2178 ts->data = (void *)ark; 2179 ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 2180 2181 ark->VecsDeltaLam = NULL; 2182 ark->VecsSensiTemp = NULL; 2183 ark->VecsSensiPTemp = NULL; 2184 2185 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2186 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2187 if (!dirk) { 2188 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 2189 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 2190 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2191 } 2192 PetscFunctionReturn(PETSC_SUCCESS); 2193 } 2194 2195 /* ------------------------------------------------------------ */ 2196 2197 static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2198 { 2199 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2200 2201 PetscFunctionBegin; 2202 PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2203 PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2204 PetscFunctionReturn(PETSC_SUCCESS); 2205 } 2206 2207 /*@C 2208 TSDIRKSetType - Set the type of `TSDIRK` scheme 2209 2210 Logically Collective 2211 2212 Input Parameters: 2213 + ts - timestepping context 2214 - dirktype - type of `TSDIRK` scheme 2215 2216 Options Database Key: 2217 . -ts_dirkimex_type - set `TSDIRK` scheme type 2218 2219 Level: intermediate 2220 2221 .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2222 @*/ 2223 PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2224 { 2225 PetscFunctionBegin; 2226 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2227 PetscAssertPointer(dirktype, 2); 2228 PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2229 PetscFunctionReturn(PETSC_SUCCESS); 2230 } 2231 2232 /*@C 2233 TSDIRKGetType - Get the type of `TSDIRK` scheme 2234 2235 Logically Collective 2236 2237 Input Parameter: 2238 . ts - timestepping context 2239 2240 Output Parameter: 2241 . dirktype - type of `TSDIRK` scheme 2242 2243 Level: intermediate 2244 2245 .seealso: [](ch_ts), `TSDIRKSetType()` 2246 @*/ 2247 PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2248 { 2249 PetscFunctionBegin; 2250 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2251 PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2252 PetscFunctionReturn(PETSC_SUCCESS); 2253 } 2254 2255 /*MC 2256 TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes 2257 2258 Level: beginner 2259 2260 Notes: 2261 The default is `TSDIRKS212`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type 2262 2263 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2264 M*/ 2265 PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2266 { 2267 PetscFunctionBegin; 2268 PetscCall(TSCreate_ARKIMEX(ts)); 2269 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2270 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2271 PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 2272 PetscFunctionReturn(PETSC_SUCCESS); 2273 } 2274