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 #include <../src/ts/impls/arkimex/arkimex.h> 15 #include <../src/ts/impls/arkimex/fsarkimex.h> 16 17 static ARKTableauLink ARKTableauList; 18 static TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3; 19 static TSDIRKType TSDIRKDefault = TSDIRKES213SAL; 20 static PetscBool TSARKIMEXRegisterAllCalled; 21 static PetscBool TSARKIMEXPackageInitialized; 22 static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec); 23 24 /*MC 25 TSARKIMEXARS122 - Second order ARK IMEX scheme, {cite}`ascher_1997` 26 27 This method has one explicit stage and one implicit stage. 28 29 Options Database Key: 30 . -ts_arkimex_type ars122 - set arkimex type to ars122 31 32 Level: advanced 33 34 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 35 M*/ 36 37 /*MC 38 TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part. 39 40 This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu. 41 42 Options Database Key: 43 . -ts_arkimex_type a2 - set arkimex type to a2 44 45 Level: advanced 46 47 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 48 M*/ 49 50 /*MC 51 TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part, {cite}`pareschi_2005` 52 53 This method has two implicit stages, and L-stable implicit scheme. 54 55 Options Database Key: 56 . -ts_arkimex_type l2 - set arkimex type to l2 57 58 Level: advanced 59 60 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 61 M*/ 62 63 /*MC 64 TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method. 65 66 This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used. 67 68 Options Database Key: 69 . -ts_arkimex_type 1bee - set arkimex type to 1bee 70 71 Level: advanced 72 73 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 74 M*/ 75 76 /*MC 77 TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part. 78 79 This method has one explicit stage and two implicit stages. The implicit part is the same as in TSARKIMEX2D and TSARKIMEX2E, but the explicit part has a larger stability region on the negative real axis. This method was provided by Emil Constantinescu. 80 81 Options Database Key: 82 . -ts_arkimex_type 2c - set arkimex type to 2c 83 84 Level: advanced 85 86 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 87 M*/ 88 89 /*MC 90 TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 91 92 This method has one explicit stage and two implicit stages. The stability function is independent of the explicit part in the infinity limit of the implicit component. This method was provided by Emil Constantinescu. 93 94 Options Database Key: 95 . -ts_arkimex_type 2d - set arkimex type to 2d 96 97 Level: advanced 98 99 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 100 M*/ 101 102 /*MC 103 TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 104 105 This method has one explicit stage and two implicit stages. It is an optimal method developed by Emil Constantinescu. 106 107 Options Database Key: 108 . -ts_arkimex_type 2e - set arkimex type to 2e 109 110 Level: advanced 111 112 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 113 M*/ 114 115 /*MC 116 TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme, {cite}`pareschi_2005` 117 118 This method has three implicit stages. 119 120 This method is referred to as SSP2-(3,3,2) in <https://arxiv.org/abs/1110.4375> 121 122 Options Database Key: 123 . -ts_arkimex_type prssp2 - set arkimex type to prssp2 124 125 Level: advanced 126 127 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 128 M*/ 129 130 /*MC 131 TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003` 132 133 This method has one explicit stage and three implicit stages. 134 135 Options Database Key: 136 . -ts_arkimex_type 3 - set arkimex type to 3 137 138 Level: advanced 139 140 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 141 M*/ 142 143 /*MC 144 TSARKIMEXARS443 - Third order ARK IMEX scheme, {cite}`ascher_1997` 145 146 This method has one explicit stage and four implicit stages. 147 148 Options Database Key: 149 . -ts_arkimex_type ars443 - set arkimex type to ars443 150 151 Level: advanced 152 153 Notes: 154 This method is referred to as ARS(4,4,3) in <https://arxiv.org/abs/1110.4375> 155 156 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 157 M*/ 158 159 /*MC 160 TSARKIMEXBPR3 - Third order ARK IMEX scheme. Referred to as ARK3 in <https://arxiv.org/abs/1110.4375> 161 162 This method has one explicit stage and four implicit stages. 163 164 Options Database Key: 165 . -ts_arkimex_type bpr3 - set arkimex type to bpr3 166 167 Level: advanced 168 169 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 170 M*/ 171 172 /*MC 173 TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`. 174 175 This method has one explicit stage and four implicit stages. 176 177 Options Database Key: 178 . -ts_arkimex_type 4 - set arkimex type to4 179 180 Level: advanced 181 182 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 183 M*/ 184 185 /*MC 186 TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part, {cite}`kennedy_2003`. 187 188 This method has one explicit stage and five implicit stages. 189 190 Options Database Key: 191 . -ts_arkimex_type 5 - set arkimex type to 5 192 193 Level: advanced 194 195 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 196 M*/ 197 198 /*MC 199 TSDIRKS212 - Second order DIRK scheme. 200 201 This method has two implicit stages with an embedded method of other 1. 202 See `TSDIRK` for additional details. 203 204 Options Database Key: 205 . -ts_dirk_type s212 - select this method. 206 207 Level: advanced 208 209 Note: 210 This is the default DIRK scheme in SUNDIALS. 211 212 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 213 M*/ 214 215 /*MC 216 TSDIRKES122SAL - First order DIRK scheme <https://arxiv.org/abs/1803.01613> 217 218 Uses backward Euler as advancing method and trapezoidal rule as embedded method. See `TSDIRK` for additional details. 219 220 Options Database Key: 221 . -ts_dirk_type es122sal - select this method. 222 223 Level: advanced 224 225 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 226 M*/ 227 228 /*MC 229 TSDIRKES213SAL - Second order DIRK scheme {cite}`kennedy2019diagonally`. Also known as TR-BDF2, see{cite}`hosea1996analysis` 230 231 See `TSDIRK` for additional details. 232 233 Options Database Key: 234 . -ts_dirk_type es213sal - select this method. 235 236 Level: advanced 237 238 Note: 239 This is the default DIRK scheme used in PETSc. 240 241 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 242 M*/ 243 244 /*MC 245 TSDIRKES324SAL - Third order DIRK scheme, {cite}`kennedy2019diagonally` 246 247 See `TSDIRK` for additional details. 248 249 Options Database Key: 250 . -ts_dirk_type es324sal - select this method. 251 252 Level: advanced 253 254 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 255 M*/ 256 257 /*MC 258 TSDIRKES325SAL - Third order DIRK scheme {cite}`kennedy2019diagonally`. 259 260 See `TSDIRK` for additional details. 261 262 Options Database Key: 263 . -ts_dirk_type es325sal - select this method. 264 265 Level: advanced 266 267 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 268 M*/ 269 270 /*MC 271 TSDIRK657A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 272 273 See `TSDIRK` for additional details. 274 275 Options Database Key: 276 . -ts_dirk_type 657a - select this method. 277 278 Level: advanced 279 280 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 281 M*/ 282 283 /*MC 284 TSDIRKES648SA - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 285 286 See `TSDIRK` for additional details. 287 288 Options Database Key: 289 . -ts_dirk_type es648sa - select this method. 290 291 Level: advanced 292 293 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 294 M*/ 295 296 /*MC 297 TSDIRK658A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 298 299 See `TSDIRK` for additional details. 300 301 Options Database Key: 302 . -ts_dirk_type 658a - select this method. 303 304 Level: advanced 305 306 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 307 M*/ 308 309 /*MC 310 TSDIRKS659A - Sixth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 311 312 See `TSDIRK` for additional details. 313 314 Options Database Key: 315 . -ts_dirk_type s659a - select this method. 316 317 Level: advanced 318 319 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 320 M*/ 321 322 /*MC 323 TSDIRK7510SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 324 325 See `TSDIRK` for additional details. 326 327 Options Database Key: 328 . -ts_dirk_type 7510sal - select this method. 329 330 Level: advanced 331 332 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 333 M*/ 334 335 /*MC 336 TSDIRKES7510SA - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 337 338 See `TSDIRK` for additional details. 339 340 Options Database Key: 341 . -ts_dirk_type es7510sa - select this method. 342 343 Level: advanced 344 345 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 346 M*/ 347 348 /*MC 349 TSDIRK759A - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 350 351 See `TSDIRK` for additional details. 352 353 Options Database Key: 354 . -ts_dirk_type 759a - select this method. 355 356 Level: advanced 357 358 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 359 M*/ 360 361 /*MC 362 TSDIRKS7511SAL - Seventh order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 363 364 See `TSDIRK` for additional details. 365 366 Options Database Key: 367 . -ts_dirk_type s7511sal - select this method. 368 369 Level: advanced 370 371 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 372 M*/ 373 374 /*MC 375 TSDIRK8614A - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 376 377 See `TSDIRK` for additional details. 378 379 Options Database Key: 380 . -ts_dirk_type 8614a - select this method. 381 382 Level: advanced 383 384 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 385 M*/ 386 387 /*MC 388 TSDIRK8616SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 389 390 See `TSDIRK` for additional details. 391 392 Options Database Key: 393 . -ts_dirk_type 8616sal - select this method. 394 395 Level: advanced 396 397 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 398 M*/ 399 400 /*MC 401 TSDIRKES8516SAL - Eighth order DIRK scheme <https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs> 402 403 See `TSDIRK` for additional details. 404 405 Options Database Key: 406 . -ts_dirk_type es8516sal - select this method. 407 408 Level: advanced 409 410 .seealso: [](ch_ts), `TSDIRK`, `TSDIRKType`, `TSDIRKSetType()` 411 M*/ 412 413 PetscErrorCode TSHasRHSFunction(TS ts, PetscBool *has) 414 { 415 TSRHSFunctionFn *func; 416 417 PetscFunctionBegin; 418 *has = PETSC_FALSE; 419 PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL)); 420 if (func) *has = PETSC_TRUE; 421 PetscFunctionReturn(PETSC_SUCCESS); 422 } 423 424 /*@C 425 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX` 426 427 Not Collective, but should be called by all processes which will need the schemes to be registered 428 429 Level: advanced 430 431 .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()` 432 @*/ 433 PetscErrorCode TSARKIMEXRegisterAll(void) 434 { 435 PetscFunctionBegin; 436 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 437 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 438 439 #define RC PetscRealConstant 440 #define s2 RC(1.414213562373095048802) /* PetscSqrtReal((PetscReal)2.0) */ 441 #define us2 RC(0.2928932188134524755992) /* 1.0-1.0/PetscSqrtReal((PetscReal)2.0); */ 442 443 /* Diagonally implicit methods */ 444 { 445 /* DIRK212, default of SUNDIALS */ 446 const PetscReal A[2][2] = { 447 {RC(1.0), RC(0.0)}, 448 {RC(-1.0), RC(1.0)} 449 }; 450 const PetscReal b[2] = {RC(0.5), RC(0.5)}; 451 const PetscReal bembed[2] = {RC(1.0), RC(0.0)}; 452 PetscCall(TSDIRKRegister(TSDIRKS212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 453 } 454 455 { 456 /* ESDIRK12 from https://arxiv.org/pdf/1803.01613.pdf */ 457 const PetscReal A[2][2] = { 458 {RC(0.0), RC(0.0)}, 459 {RC(0.0), RC(1.0)} 460 }; 461 const PetscReal b[2] = {RC(0.0), RC(1.0)}; 462 const PetscReal bembed[2] = {RC(0.5), RC(0.5)}; 463 PetscCall(TSDIRKRegister(TSDIRKES122SAL, 1, 2, &A[0][0], b, NULL, bembed, 1, b)); 464 } 465 466 { 467 /* ESDIRK213L[2]SA from KC16. 468 TR-BDF2 from Hosea and Shampine 469 ESDIRK23 in https://arxiv.org/pdf/1803.01613.pdf */ 470 const PetscReal A[3][3] = { 471 {RC(0.0), RC(0.0), RC(0.0)}, 472 {us2, us2, RC(0.0)}, 473 {s2 / RC(4.0), s2 / RC(4.0), us2 }, 474 }; 475 const PetscReal b[3] = {s2 / RC(4.0), s2 / RC(4.0), us2}; 476 const PetscReal bembed[3] = {(RC(1.0) - s2 / RC(4.0)) / RC(3.0), (RC(3.0) * s2 / RC(4.0) + RC(1.0)) / RC(3.0), us2 / RC(3.0)}; 477 PetscCall(TSDIRKRegister(TSDIRKES213SAL, 2, 3, &A[0][0], b, NULL, bembed, 1, b)); 478 } 479 480 { 481 #define g RC(0.43586652150845899941601945) 482 #define g2 PetscSqr(g) 483 #define g3 g *g2 484 #define g4 PetscSqr(g2) 485 #define g5 g *g4 486 #define c3 RC(1.0) 487 #define a32 c3 *(c3 - RC(2.0) * g) / (RC(4.0) * g) 488 #define b2 (-RC(2.0) + RC(3.0) * c3 + RC(6.0) * g * (RC(1.0) - c3)) / (RC(12.0) * g * (c3 - RC(2.0) * g)) 489 #define b3 (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) / (RC(3.0) * c3 * (c3 - RC(2.0) * g)) 490 #if 0 491 /* This is for c3 = 3/5 */ 492 #define bh2 \ 493 c3 * (-RC(1.0) + RC(6.0) * g - RC(23.0) * g3 + RC(12.0) * g4 - RC(6.0) * g5) / (RC(4.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) + (RC(3.0) - RC(27.0) * g + RC(68.0) * g2 - RC(55.0) * g3 + RC(21.0) * g4 - RC(6.0) * g5) / (RC(2.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) 494 #define bh3 -g * (-RC(2.0) + RC(21.0) * g - RC(68.0) * g2 + RC(79.0) * g3 - RC(33.0) * g4 + RC(12.0) * g5) / (RC(2.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) 495 #define bh4 -RC(3.0) * g2 * (-RC(1.0) + RC(4.0) * g - RC(2.0) * g2 + g3) / (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) 496 #else 497 /* This is for c3 = 1.0 */ 498 #define bh2 a32 499 #define bh3 g 500 #define bh4 RC(0.0) 501 #endif 502 /* ESDIRK3(2I)4L[2]SA from KC16 with c3 = 1.0 */ 503 /* Given by Kvaerno https://link.springer.com/article/10.1023/b:bitn.0000046811.70614.38 */ 504 const PetscReal A[4][4] = { 505 {RC(0.0), RC(0.0), RC(0.0), RC(0.0)}, 506 {g, g, RC(0.0), RC(0.0)}, 507 {c3 - a32 - g, a32, g, RC(0.0)}, 508 {RC(1.0) - b2 - b3 - g, b2, b3, g }, 509 }; 510 const PetscReal b[4] = {RC(1.0) - b2 - b3 - g, b2, b3, g}; 511 const PetscReal bembed[4] = {RC(1.0) - bh2 - bh3 - bh4, bh2, bh3, bh4}; 512 PetscCall(TSDIRKRegister(TSDIRKES324SAL, 3, 4, &A[0][0], b, NULL, bembed, 1, b)); 513 #undef g 514 #undef g2 515 #undef g3 516 #undef c3 517 #undef a32 518 #undef b2 519 #undef b3 520 #undef bh2 521 #undef bh3 522 #undef bh4 523 } 524 525 { 526 /* ESDIRK3(2I)5L[2]SA from KC16 */ 527 const PetscReal A[5][5] = { 528 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 529 {RC(9.0) / RC(40.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0), RC(0.0) }, 530 {RC(19.0) / RC(72.0), RC(14.0) / RC(45.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0) }, 531 {RC(3337.0) / RC(11520.0), RC(233.0) / RC(720.0), RC(207.0) / RC(1280.0), RC(9.0) / RC(40.0), RC(0.0) }, 532 {RC(7415.0) / RC(34776.0), RC(9920.0) / RC(30429.0), RC(4845.0) / RC(9016.0), -RC(5827.0) / RC(19320.0), RC(9.0) / RC(40.0)}, 533 }; 534 const PetscReal b[5] = {RC(7415.0) / RC(34776.0), RC(9920.0) / RC(30429.0), RC(4845.0) / RC(9016.0), -RC(5827.0) / RC(19320.0), RC(9.0) / RC(40.0)}; 535 const PetscReal bembed[5] = {RC(23705.0) / RC(104328.0), RC(29720.0) / RC(91287.0), RC(4225.0) / RC(9016.0), -RC(69304987.0) / RC(337732920.0), RC(42843.0) / RC(233080.0)}; 536 PetscCall(TSDIRKRegister(TSDIRKES325SAL, 3, 5, &A[0][0], b, NULL, bembed, 1, b)); 537 } 538 539 { 540 // DIRK(6,6)[1]A[(7,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 541 const PetscReal A[7][7] = { 542 {RC(0.303487844706747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 543 {RC(-0.279756492709814), RC(0.500032236020747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 544 {RC(0.280583215743895), RC(-0.438560061586751), RC(0.217250734515736), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 545 {RC(-0.0677678738539846), RC(0.984312781232293), RC(-0.266720192540149), RC(0.2476680834526), RC(0.0), RC(0.0), RC(0.0) }, 546 {RC(0.125671616147993), RC(-0.995401751002415), RC(0.761333109549059), RC(-0.210281837202208), RC(0.866743712636936), RC(0.0), RC(0.0) }, 547 {RC(-0.368056238801488), RC(-0.999928082701516), RC(0.534734253232519), RC(-0.174856916279082), RC(0.615007160285509), RC(0.696549912132029), RC(0.0) }, 548 {RC(-0.00570546839653984), RC(-0.113110431835656), RC(-0.000965563207671587), RC(-0.000130490084629567), RC(0.00111737736895673), RC(-0.279385587378871), RC(0.618455906845342)} 549 }; 550 const PetscReal b[7] = {RC(0.257561510484877), RC(0.234281287047716), RC(0.126658904241469), RC(0.252363215441784), RC(0.396701083526306), RC(-0.267566000742152), RC(0.0)}; 551 const PetscReal bembed[7] = {RC(0.257561510484945), RC(0.387312822934391), RC(0.126658904241468), RC(0.252363215441784), RC(0.396701083526306), RC(-0.267566000742225), RC(-0.153031535886669)}; 552 PetscCall(TSDIRKRegister(TSDIRK657A, 6, 7, &A[0][0], b, NULL, bembed, 1, b)); 553 } 554 { 555 // ESDIRK(8,6)[2]SA[(8,4)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 556 const PetscReal A[8][8] = { 557 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 558 {RC(0.333222149217725), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 559 {RC(0.0639743773182214), RC(-0.0830330224410214), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 560 {RC(-0.728522201369326), RC(-0.210414479522485), RC(0.532519916559342), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 561 {RC(-0.175135269272067), RC(0.666675582067552), RC(-0.304400907370867), RC(0.656797712445756), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0) }, 562 {RC(0.222695802705462), RC(-0.0948971794681061), RC(-0.0234336346686545), RC(-0.45385925012042), RC(0.0283910313826958), RC(0.333222149217725), RC(0.0), RC(0.0) }, 563 {RC(-0.132534078051299), RC(0.702597935004879), RC(-0.433316453128078), RC(0.893717488547587), RC(0.057381454791406), RC(-0.207798411552402), RC(0.333222149217725), RC(0.0) }, 564 {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)} 565 }; 566 const PetscReal b[8] = {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)}; 567 const PetscReal bembed[8] = {RC(0.0), RC(0.292331064554014), RC(0.409676102283681), RC(-0.002094718084982), RC(-0.282771520835975), RC(0.113862336644901), RC(0.181973572260693), RC(0.287023163177669)}; 568 PetscCall(TSDIRKRegister(TSDIRKES648SA, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 569 } 570 { 571 // DIRK(8,6)[1]SAL[(8,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 572 const PetscReal A[8][8] = { 573 {RC(0.477264457385826), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 574 {RC(-0.197052588415002), RC(0.476363428459584), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 575 {RC(-0.0347674430372966), RC(0.633051807335483), RC(0.193634310075028), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 576 {RC(0.0967797668578702), RC(-0.193533526466535), RC(-0.000207622945800473), RC(0.159572204849431), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 577 {RC(0.162527231819875), RC(-0.249672513547382), RC(-0.0459079972041795), RC(0.36579476400859), RC(0.255752838307699), RC(0.0), RC(0.0), RC(0.0) }, 578 {RC(-0.00707603197171262), RC(0.846299854860295), RC(0.344020016925018), RC(-0.0720926054548865), RC(-0.215492331980875), RC(0.104341097622161), RC(0.0), RC(0.0) }, 579 {RC(0.00176857935179744), RC(0.0779960013127515), RC(0.303333277564557), RC(0.213160806732836), RC(0.351769320319038), RC(-0.381545894386538), RC(0.433517909105558), RC(0.0) }, 580 {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)} 581 }; 582 const PetscReal b[8] = {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)}; 583 const PetscReal bembed[8] = {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.103483943222765), RC(-0.0103721771642262), RC(0.177302191576001)}; 584 PetscCall(TSDIRKRegister(TSDIRK658A, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 585 } 586 { 587 // SDIRK(9,6)[1]SAL[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 588 const PetscReal A[9][9] = { 589 {RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 590 {RC(-0.0903514856119419), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 591 {RC(0.172952039138937), RC(-0.35365501036282), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 592 {RC(0.511999875919193), RC(0.0289640332201925), RC(-0.0144030945657094), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 593 {RC(0.00465303495506782), RC(-0.075635818766597), RC(0.217273030786712), RC(-0.0206519428725472), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 594 {RC(0.896145501762472), RC(0.139267327700498), RC(-0.186920979752805), RC(0.0672971012371723), RC(-0.350891963442176), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0) }, 595 {RC(0.552959701885751), RC(-0.439360579793662), RC(0.333704002325091), RC(-0.0339426520778416), RC(-0.151947445912595), RC(0.0213825661026943), RC(0.218127781944908), RC(0.0), RC(0.0) }, 596 {RC(0.631360374036476), RC(0.724733619641466), RC(-0.432170625425258), RC(0.598611382182477), RC(-0.709087197034345), RC(-0.483986685696934), RC(0.378391562905131), RC(0.218127781944908), RC(0.0) }, 597 {RC(0.0), RC(-0.15504452530869), RC(0.194518478660789), RC(0.63515640279203), RC(0.81172278664173), RC(0.110736108691585), RC(-0.495304692414479), RC(-0.319912341007872), RC(0.218127781944908)} 598 }; 599 const PetscReal b[9] = {RC(0.0), RC(-0.15504452530869), RC(0.194518478660789), RC(0.63515640279203), RC(0.81172278664173), RC(0.110736108691585), RC(-0.495304692414479), RC(-0.319912341007872), RC(0.218127781944908)}; 600 const PetscReal bembed[9] = {RC(3.62671059311602e-16), RC(0.0736615558278942), RC(0.103527397262229), RC(1.00247481935499), RC(0.361377289250057), RC(-0.785425929961365), RC(-0.0170499047960784), RC(0.296321252214769), RC(-0.0348864791524953)}; 601 PetscCall(TSDIRKRegister(TSDIRKS659A, 6, 9, &A[0][0], b, NULL, bembed, 1, b)); 602 } 603 { 604 // DIRK(10,7)[1]SAL[(10,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 605 const PetscReal A[10][10] = { 606 {RC(0.233704632125264), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 607 {RC(-0.0739324813149407), RC(0.200056838146104), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 608 {RC(0.0943790344044812), RC(0.264056067701605), RC(0.133245202456465), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 609 {RC(0.269084810601201), RC(-0.503479002548384), RC(-0.00486736469695022), RC(0.251518716213569), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 610 {RC(0.145665801918423), RC(0.204983170463176), RC(0.407154634069484), RC(-0.0121039135200389), RC(0.190243622486334), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 611 {RC(0.985450198547345), RC(0.806942652811456), RC(-0.808130934167263), RC(-0.669035819439391), RC(0.0269384406756128), RC(0.462144080607327), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 612 {RC(0.163902957809563), RC(0.228315094960095), RC(0.0745971021260249), RC(0.000509793400156559), RC(0.0166533681378294), RC(-0.0229383879045797), RC(0.103505486637336), RC(0.0), RC(0.0), RC(0.0) }, 613 {RC(-0.162694156858437), RC(0.0453478837428434), RC(0.997443481211424), RC(0.200251514941093), RC(-0.000161755458839048), RC(-0.0848134335980281), RC(-0.36438666566666), RC(0.158604420136055), RC(0.0), RC(0.0) }, 614 {RC(0.200733156477425), RC(0.239686443444433), RC(0.303837014418929), RC(-0.0534390596279896), RC(0.0314067599640569), RC(-0.00764032790448536), RC(0.0609191260198661), RC(-0.0736319201590642), RC(0.204602530607021), RC(0.0) }, 615 {RC(0.0), RC(0.235563761744267), RC(0.658651488684319), RC(0.0308877804992098), RC(-0.906514945595336), RC(-0.0248488551739974), RC(-0.309967582365257), RC(0.191663316925525), RC(0.923933712199542), RC(0.200631323081727)} 616 }; 617 const PetscReal b[10] = {RC(0.0), RC(0.235563761744267), RC(0.658651488684319), RC(0.0308877804992098), RC(-0.906514945595336), RC(-0.0248488551739974), RC(-0.309967582365257), RC(0.191663316925525), RC(0.923933712199542), RC(0.200631323081727)}; 618 const PetscReal bembed[10] = 619 {RC(0.0), RC(0.222929376486581), RC(0.950668440138169), RC(0.0342694607044032), RC(0.362875840545746), RC(0.223572979288581), RC(-0.764361723526727), RC(0.563476909230026), RC(-0.690896961894185), RC(0.0974656790270323)}; 620 PetscCall(TSDIRKRegister(TSDIRK7510SAL, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 621 } 622 { 623 // ESDIRK(10,7)[2]SA[(10,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 624 const PetscReal A[10][10] = { 625 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 626 {RC(0.210055790203419), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 627 {RC(0.255781739921086), RC(0.239850916980976), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 628 {RC(0.286789624880437), RC(0.230494748834778), RC(0.263925149885491), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 629 {RC(-0.0219118128774335), RC(0.897684380345907), RC(-0.657954605498907), RC(0.124962304722633), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 630 {RC(-0.065614879584776), RC(-0.0565630711859497), RC(0.0254881105065311), RC(-0.00368981790650006), RC(-0.0115178258446329), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 631 {RC(0.399860851232098), RC(0.915588469718705), RC(-0.0758429094934412), RC(-0.263369154872759), RC(0.719687583564526), RC(-0.787410407015369), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0) }, 632 {RC(0.51693616104628), RC(1.00000540846973), RC(-0.0485110663289207), RC(-0.315208041581942), RC(0.749742806451587), RC(-0.990975090921248), RC(0.0159279583407308), RC(0.210055790203419), RC(0.0), RC(0.0) }, 633 {RC(-0.0303062129076945), RC(-0.297035174659034), RC(0.184724697462164), RC(-0.0351876079516183), RC(-0.00324668230690761), RC(0.216151004053531), RC(-0.126676252098317), RC(0.114040254365262), RC(0.210055790203419), RC(0.0) }, 634 {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)} 635 }; 636 const PetscReal b[10] = {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), 637 RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)}; 638 const PetscReal bembed[10] = {RC(-0.015494246543626), RC(0.167657963820093), RC(0.269858958144236), RC(-0.0443258997755156), RC(0.150049236875266), 639 RC(0.259452082755846), RC(0.244624573502521), RC(-0.215528446920284), RC(0.0487601760292619), RC(0.134945602112201)}; 640 PetscCall(TSDIRKRegister(TSDIRKES7510SA, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 641 } 642 { 643 // DIRK(9,7)[1]A[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 644 const PetscReal A[9][9] = { 645 {RC(0.179877789855839), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 646 {RC(-0.100405844885157), RC(0.214948590644819), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 647 {RC(0.112251360198995), RC(-0.206162139150298), RC(0.125159642941958), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 648 {RC(-0.0335164000768257), RC(0.999942349946143), RC(-0.491470853833294), RC(0.19820086325566), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 649 {RC(-0.0417345265478321), RC(0.187864510308215), RC(0.0533789224305102), RC(-0.00822060284862916), RC(0.127670843671646), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 650 {RC(-0.0278257925239257), RC(0.600979340683382), RC(-0.242632273241134), RC(-0.11318753652081), RC(0.164326917632931), RC(0.284116597781395), RC(0.0), RC(0.0), RC(0.0) }, 651 {RC(0.041465583858922), RC(0.429657872601836), RC(-0.381323410582524), RC(0.391934277498434), RC(-0.245918275501241), RC(-0.35960669741231), RC(0.184000022289158), RC(0.0), RC(0.0) }, 652 {RC(-0.105565651574538), RC(-0.0557833155018609), RC(0.358967568942643), RC(-0.13489263413921), RC(0.129553247260677), RC(0.0992493795371489), RC(-0.15716610563461), RC(0.17918862279814), RC(0.0) }, 653 {RC(0.00439696079965225), RC(0.960250486570491), RC(0.143558372286706), RC(0.0819015241056593), RC(0.999562318563625), RC(0.325203439314358), RC(-0.679013149331228), RC(-0.990589559837246), RC(0.0773648037639896)} 654 }; 655 656 const PetscReal b[9] = {RC(0.0), RC(0.179291520437966), RC(0.115310295273026), RC(-0.857943261453138), RC(0.654911318641998), RC(1.18713633508094), RC(-0.0949482361570542), RC(-0.37661430946407), RC(0.19285633764033)}; 657 const PetscReal bembed[9] = {RC(0.0), RC(0.1897135479408), RC(0.127461414808862), RC(-0.835810807663404), RC(0.665114177777166), RC(1.16481046518346), RC(-0.11661858889792), RC(-0.387303251022099), RC(0.192633041873135)}; 658 PetscCall(TSDIRKRegister(TSDIRK759A, 7, 9, &A[0][0], b, NULL, bembed, 1, b)); 659 } 660 { 661 // SDIRK(11,7)[1]SAL[(11,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 662 const PetscReal A[11][11] = { 663 {RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 664 {RC(-0.082947368165267), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 665 {RC(0.483452690540751), RC(0.0), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 666 {RC(0.771076453481321), RC(-0.22936926341842), RC(0.289733373208823), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 667 {RC(0.0329683054968892), RC(-0.162397421903366), RC(0.000951777538562805), RC(0.0), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 668 {RC(0.265888743485945), RC(0.606743151103931), RC(0.173443800537369), RC(-0.0433968261546912), RC(-0.385211017224481), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 669 {RC(0.220662294551146), RC(-0.0465078507657608), RC(-0.0333111995282464), RC(0.011801580836998), RC(0.169480801030105), RC(-0.0167974432139385), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 670 {RC(0.323099728365267), RC(0.0288371831672575), RC(-0.0543404318773196), RC(0.0137765831431662), RC(0.0516799019060702), RC(-0.0421359763835713), RC(0.181297932037826), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0) }, 671 {RC(-0.164226696476538), RC(0.187552004946792), RC(0.0628674420973025), RC(-0.0108886582703428), RC(-0.0117628641717889), RC(0.0432176880867965), RC(-0.0315206836275473), RC(-0.0846007021638797), RC(0.200252661187742), RC(0.0), RC(0.0) }, 672 {RC(0.651428598623771), RC(-0.10208078475356), RC(0.198305701801888), RC(-0.0117354096673789), RC(-0.0440385966743686), RC(-0.0358364455795087), RC(-0.0075408087654097), RC(0.160320941654639), RC(0.017940248694499), RC(0.200252661187742), RC(0.0) }, 673 {RC(0.0), RC(-0.266259448580236), RC(-0.615982357748271), RC(0.561474126687165), RC(0.266911112787025), RC(0.219775952207137), RC(0.387847665451514), RC(0.612483137773236), RC(0.330027015806089), RC(-0.6965298655714), RC(0.200252661187742)} 674 }; 675 const PetscReal b[11] = 676 {RC(0.0), RC(-0.266259448580236), RC(-0.615982357748271), RC(0.561474126687165), RC(0.266911112787025), RC(0.219775952207137), RC(0.387847665451514), RC(0.612483137773236), RC(0.330027015806089), RC(-0.6965298655714), RC(0.200252661187742)}; 677 const PetscReal bembed[11] = 678 {RC(0.0), RC(0.180185524442613), RC(-0.628869710835338), RC(0.186185675988647), RC(0.0484716652630425), RC(0.203927720607141), RC(0.44041662512573), RC(0.615710527731245), RC(0.0689648839032607), RC(-0.253599870605903), RC(0.138606958379488)}; 679 PetscCall(TSDIRKRegister(TSDIRKS7511SAL, 7, 11, &A[0][0], b, NULL, bembed, 1, b)); 680 } 681 { 682 // DIRK(13,8)[1]A[(14,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 683 const PetscReal A[14][14] = { 684 {RC(0.421050745442291), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 685 {RC(-0.0761079419591268), RC(0.264353986580857), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 686 {RC(0.0727106904170694), RC(-0.204265976977285), RC(0.181608196544136), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 687 {RC(0.55763054816611), RC(-0.409773579543499), RC(0.510926516886944), RC(0.259892204518476), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 688 {RC(0.0228083864844437), RC(-0.445569051836454), RC(-0.0915242778636248), RC(0.00450055909321655), RC(0.6397807199983), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 689 {RC(-0.135945849505152), RC(0.0946509646963754), RC(-0.236110197279175), RC(0.00318944206456517), RC(0.255453021028118), RC(0.174805219173446), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 690 {RC(-0.147960260670772), RC(-0.402188192230535), RC(-0.703014530043888), RC(0.00941974677418186), RC(0.885747111289207), RC(0.261314066449028), RC(0.16307697503668), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 691 {RC(0.165597241042244), RC(0.824182962188923), RC(-0.0280136160783609), RC(0.282372386631758), RC(-0.957721354131182), RC(0.489439550159977), RC(0.170094415598103), RC(0.0522519785718563), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 692 {RC(0.0335292011495618), RC(0.575750388029166), RC(0.223289855356637), RC(-0.00317458833242804), RC(-0.112890382135193), RC(-0.419809267954284), RC(0.0466136902102104), RC(-0.00115413813041085), RC(0.109685363692383), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 693 {RC(-0.0512616878252355), RC(0.699261265830807), RC(-0.117939611738769), RC(0.0021745241931243), RC(-0.00932826702640947), RC(-0.267575057469428), RC(0.126949139814065), RC(0.00330353204502163), RC(0.185949445053766), RC(0.0938215615963721), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 694 {RC(-0.106521517960343), RC(0.41835889096168), RC(0.353585905881916), RC(-0.0746474161579599), RC(-0.015450626460289), RC(-0.46224659192275), RC(-0.0576406327329181), RC(-0.00712066942504018), RC(0.377776558014452), RC(0.36890054338294), RC(0.0618488746331837), RC(0.0), RC(0.0), RC(0.0) }, 695 {RC(-0.163079104890997), RC(0.644561721693806), RC(0.636968661639572), RC(-0.122346720085377), RC(-0.333062564990312), RC(-0.3054226490478), RC(-0.357820712828352), RC(-0.0125510510334706), RC(0.371263681186311), RC(0.371979640363694), RC(0.0531090658708968), RC(0.0518279459132049), RC(0.0), RC(0.0) }, 696 {RC(0.579993784455521), RC(-0.188833728676494), RC(0.999975696843775), RC(0.0572810855901161), RC(-0.264374735003671), RC(0.165091739976854), RC(-0.546675809010452), RC(-0.0283821822291982), RC(-0.102639860418374), RC(-0.0343251040446405), RC(0.4762598462591), RC(-0.304153104931261), RC(0.0953911855943621), RC(0.0) }, 697 {RC(0.0848552694007844), RC(0.287193912340074), RC(0.543683503004232), RC(-0.081311059300692), RC(-0.0328661289388557), RC(-0.323456834372922), RC(-0.240378871658975), RC(-0.0189913019930369), RC(0.220663114082036), RC(0.253029984360864), RC(0.252011799370563), RC(-0.154882222605423), RC(0.0315202264687415), RC(0.0514095812104714)} 698 }; 699 const PetscReal b[14] = {RC(0.0), RC(0.516650324205117), RC(0.0773227217357826), RC(-0.12474204666975), RC(-0.0241052115180679), RC(-0.325821145180359), RC(0.0907237460123951), RC(0.0459271880596652), RC(0.221012259404702), RC(0.235510906761942), RC(0.491109674204385), RC(-0.323506525837343), RC(0.119918108821531), RC(0.0)}; 700 const PetscReal bembed[14] = {RC(2.32345691433618e-16), RC(0.499150900944401), RC(0.080991997189243), RC(-0.0359440417166322), RC(-0.0258910397441454), RC(-0.304540350278636), RC(0.0836627473632563), 701 RC(0.0417664613347638), RC(0.223636394275293), RC(0.231569156867596), RC(0.240526201277663), RC(-0.222933582911926), RC(-0.0111479879597561), RC(0.19915314335888)}; 702 PetscCall(TSDIRKRegister(TSDIRK8614A, 8, 14, &A[0][0], b, NULL, bembed, 1, b)); 703 } 704 { 705 // DIRK(15,8)[1]SAL[(16,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 706 const PetscReal A[16][16] = { 707 {RC(0.498904981271193), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 708 {RC(-0.303806037341816), RC(0.886299445992379), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 709 {RC(-0.581440223471476), RC(0.371003719460259), RC(0.43844717752802), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 710 {RC(0.531852638870051), RC(-0.339363014907108), RC(0.422373239795441), RC(0.223854203543397), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 711 {RC(0.118517891868867), RC(-0.0756235584174296), RC(-0.0864284870668712), RC(0.000536692838658312), RC(0.10101418329932), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 712 {RC(0.218733626116401), RC(-0.139568928299635), RC(0.30473612813488), RC(0.00354038623073564), RC(0.0932085751160559), RC(0.140161806097591), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 713 {RC(0.0692944686081835), RC(-0.0442152168939502), RC(-0.0903375348855603), RC(0.00259030241156141), RC(0.204514233679515), RC(-0.0245383758960002), RC(0.199289437094059), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 714 {RC(0.990640016505571), RC(-0.632104756315967), RC(0.856971425234221), RC(0.174494099232246), RC(-0.113715829680145), RC(-0.151494045307366), RC(-0.438268629569005), RC(0.120578398912139), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 715 {RC(-0.099415677713136), RC(0.211832014309207), RC(-0.245998265866888), RC(-0.182249672235861), RC(0.167897635713799), RC(0.212850335030069), RC(-0.391739299440123), RC(-0.0118718506876767), RC(0.526293701659093), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 716 {RC(0.383983914845461), RC(-0.245011361219604), RC(0.46717278554955), RC(-0.0361272447593202), RC(0.0742234660511333), RC(-0.0474816271948766), RC(-0.229859978525756), RC(0.0516283729206322), RC(0.0), RC(0.193823890777594), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 717 {RC(0.0967855003180134), RC(-0.0481037037916184), RC(0.191268138832434), RC(0.234977164564126), RC(0.0620265921753097), RC(0.403432826534738), RC(0.152403846687238), RC(-0.118420429237746), RC(0.0582141598685892), RC(-0.13924540906863), RC(0.106661313117545), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 718 {RC(0.133941307432055), RC(-0.0722076602896254), RC(0.217086297689275), RC(0.00495499602192887), RC(0.0306090174933995), RC(0.26483526755746), RC(0.204442440745605), RC(0.196883395136708), RC(0.056527012583996), RC(-0.150216381356784), RC(-0.217209415757333), RC(0.330353722743315), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 719 {RC(0.157014274561299), RC(-0.0883810256381874), RC(0.117193033885034), RC(-0.0362304243769466), RC(0.0169030211466111), RC(-0.169835753576141), RC(0.399749979234113), RC(0.31806704093008), RC(0.050340008347693), RC(0.120284837472214), RC(-0.235313193645423), RC(0.232488522208926), RC(0.117719679450729), RC(0.0), RC(0.0), RC(0.0) }, 720 {RC(0.00276453816875833), RC(-0.00366028255231782), RC(-0.331078914515559), RC(0.623377549031949), RC(0.167618142989491), RC(0.0748467945312516), RC(0.797629286699677), RC(-0.390714256799583), RC(-0.00808553925131555), RC(0.014840324980952), RC(-0.0856180410248133), RC(0.602943304937827), RC(-0.5771359338496), RC(0.112273026653282), RC(0.0), RC(0.0) }, 721 {RC(0.0), RC(0.0), RC(0.085283971980307), RC(0.51334393454179), RC(0.144355978013514), RC(0.255379109487853), RC(0.225075750790524), RC(-0.343241323394982), RC(0.0), RC(0.0798250392218852), RC(0.0528824734082655), RC(-0.0830350888900362), RC(0.022567388707279), RC(-0.0592631119040204), RC(0.106825878037621), RC(0.0) }, 722 {RC(0.173784481207652), RC(-0.110887906116241), RC(0.190052513365204), RC(-0.0688345422674029), RC(0.10326505079603), RC(0.267127097115219), RC(0.141703423176897), RC(0.0117966866651728), RC(-6.65650091812762e-15), RC(-0.0213725083662519), RC(-0.00931148598712566), RC(-0.10007679077114), RC(0.123471797451553), RC(0.00203684241073055), RC(-0.0294320891781173), RC(0.195746619921528)} 723 }; 724 const PetscReal b[16] = {RC(0.0), RC(0.0), RC(0.085283971980307), RC(0.51334393454179), RC(0.144355978013514), RC(0.255379109487853), RC(0.225075750790524), RC(-0.343241323394982), RC(0.0), RC(0.0798250392218852), RC(0.0528824734082655), RC(-0.0830350888900362), RC(0.022567388707279), RC(-0.0592631119040204), RC(0.106825878037621), RC(0.0)}; 725 const PetscReal bembed[16] = {RC(-1.31988512519898e-15), RC(7.53606601764004e-16), RC(0.0886789133915965), RC(0.0968726531622137), RC(0.143815375874267), RC(0.335214773313601), RC(0.221862366978063), RC(-0.147408947987273), 726 RC(4.16297599203445e-16), RC(0.000727276166520566), RC(-0.00284892677941246), RC(0.00512492274297611), RC(-0.000275595071215218), RC(0.0136014719350733), RC(0.0165190013607726), RC(0.228116714912817)}; 727 PetscCall(TSDIRKRegister(TSDIRK8616SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 728 } 729 { 730 // ESDIRK(16,8)[2]SAL[(16,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 731 const PetscReal A[16][16] = { 732 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 733 {RC(0.117318819358521), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 734 {RC(0.0557014605974616), RC(0.385525646638742), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 735 {RC(0.063493276428895), RC(0.373556126263681), RC(0.0082994166438953), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 736 {RC(0.0961351856230088), RC(0.335558324517178), RC(0.207077765910132), RC(-0.0581917140797146), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 737 {RC(0.0497669214238319), RC(0.384288616546039), RC(0.0821728117583936), RC(0.120337007107103), RC(0.202262782645888), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 738 {RC(0.00626710666809847), RC(0.496491452640725), RC(-0.111303249827358), RC(0.170478821683603), RC(0.166517073971103), RC(-0.0328669811542241), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 739 {RC(0.0463439767281591), RC(0.00306724391019652), RC(-0.00816305222386205), RC(-0.0353302599538294), RC(0.0139313601702569), RC(-0.00992014507967429), RC(0.0210087909090165), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 740 {RC(0.111574049232048), RC(0.467639166482209), RC(0.237773114804619), RC(0.0798895699267508), RC(0.109580615914593), RC(0.0307353103825936), RC(-0.0404391509541147), RC(-0.16942110744293), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 741 {RC(-0.0107072484863877), RC(-0.231376703354252), RC(0.017541113036611), RC(0.144871527682418), RC(-0.041855459769806), RC(0.0841832168332261), RC(-0.0850020937282192), RC(0.486170343825899), RC(-0.0526717116822739), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 742 {RC(-0.0142238262314935), RC(0.14752923682514), RC(0.238235830732566), RC(0.037950291904103), RC(0.252075123381518), RC(0.0474266904224567), RC(-0.00363139069342027), RC(0.274081442388563), RC(-0.0599166970745255), RC(-0.0527138812389185), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 743 {RC(-0.11837020183211), RC(-0.635712481821264), RC(0.239738832602538), RC(0.330058936651707), RC(-0.325784087988237), RC(-0.0506514314589253), RC(-0.281914404487009), RC(0.852596345144291), RC(0.651444614298805), RC(-0.103476387303591), RC(-0.354835880209975), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 744 {RC(-0.00458164025442349), RC(0.296219694015248), RC(0.322146049419995), RC(0.15917778285238), RC(0.284864871688843), RC(0.185509526463076), RC(-0.0784621067883274), RC(0.166312223692047), RC(-0.284152486083397), RC(-0.357125104338944), RC(0.078437074055306), RC(0.0884129667114481), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0) }, 745 {RC(-0.0545561913848106), RC(0.675785423442753), RC(0.423066443201941), RC(-0.000165300126841193), RC(0.104252994793763), RC(-0.105763019303021), RC(-0.15988308809318), RC(0.0515050001032011), RC(0.56013979290924), RC(-0.45781539708603), RC(-0.255870699752664), RC(0.026960254296416), RC(-0.0721245985053681), RC(0.117318819358521), RC(0.0), RC(0.0) }, 746 {RC(0.0649253995775223), RC(-0.0216056457922249), RC(-0.073738139377975), RC(0.0931033310077225), RC(-0.0194339577299149), RC(-0.0879623837313009), RC(0.057125517179467), RC(0.205120850488097), RC(0.132576503537441), RC(0.489416890627328), RC(-0.1106765720501), RC(-0.081038793996096), RC(0.0606031613503788), RC(-0.00241467937442272), RC(0.117318819358521), RC(0.0) }, 747 {RC(0.0459979286336779), RC(0.0780075394482806), RC(0.015021874148058), RC(0.195180277284195), RC(-0.00246643310153235), RC(0.0473977117068314), RC(-0.0682773558610363), RC(0.19568019123878), RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)} 748 }; 749 const PetscReal b[16] = {RC(0.0459979286336779), RC(0.0780075394482806), RC(0.015021874148058), RC(0.195180277284195), RC(-0.00246643310153235), RC(0.0473977117068314), RC(-0.0682773558610363), RC(0.19568019123878), 750 RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)}; 751 const PetscReal bembed[16] = {RC(0.0603373529853206), RC(0.175453809423998), RC(0.0537707777611352), RC(0.195309248607308), RC(0.0135893741970232), RC(-0.0221160259296707), RC(-0.00726526156430691), RC(0.102961059369124), 752 RC(0.000900215457460583), RC(0.0547959465692338), RC(-0.334995726863153), RC(0.0464409662093384), RC(0.301388101652194), RC(0.00524851570622031), RC(0.229538601845236), RC(0.124643044573514)}; 753 PetscCall(TSDIRKRegister(TSDIRKES8516SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 754 } 755 756 /* Additive methods */ 757 { 758 const PetscReal A[3][3] = { 759 {0.0, 0.0, 0.0}, 760 {0.0, 0.0, 0.0}, 761 {0.0, 0.5, 0.0} 762 }; 763 const PetscReal At[3][3] = { 764 {1.0, 0.0, 0.0}, 765 {0.0, 0.5, 0.0}, 766 {0.0, 0.5, 0.5} 767 }; 768 const PetscReal b[3] = {0.0, 0.5, 0.5}; 769 const PetscReal bembedt[3] = {1.0, 0.0, 0.0}; 770 PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 771 } 772 { 773 const PetscReal A[2][2] = { 774 {0.0, 0.0}, 775 {0.5, 0.0} 776 }; 777 const PetscReal At[2][2] = { 778 {0.0, 0.0}, 779 {0.0, 0.5} 780 }; 781 const PetscReal b[2] = {0.0, 1.0}; 782 const PetscReal bembedt[2] = {0.5, 0.5}; 783 /* binterpt[2][2] = {{1.0,-1.0},{0.0,1.0}}; second order dense output has poor stability properties and hence it is not currently in use */ 784 PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 785 } 786 { 787 const PetscReal A[2][2] = { 788 {0.0, 0.0}, 789 {1.0, 0.0} 790 }; 791 const PetscReal At[2][2] = { 792 {0.0, 0.0}, 793 {0.5, 0.5} 794 }; 795 const PetscReal b[2] = {0.5, 0.5}; 796 const PetscReal bembedt[2] = {0.0, 1.0}; 797 /* binterpt[2][2] = {{1.0,-0.5},{0.0,0.5}} second order dense output has poor stability properties and hence it is not currently in use */ 798 PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 799 } 800 { 801 const PetscReal A[2][2] = { 802 {0.0, 0.0}, 803 {1.0, 0.0} 804 }; 805 const PetscReal At[2][2] = { 806 {us2, 0.0}, 807 {1.0 - 2.0 * us2, us2} 808 }; 809 const PetscReal b[2] = {0.5, 0.5}; 810 const PetscReal bembedt[2] = {0.0, 1.0}; 811 const PetscReal binterpt[2][2] = { 812 {(us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))}, 813 {1 - (us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))} 814 }; 815 const PetscReal binterp[2][2] = { 816 {1.0, -0.5}, 817 {0.0, 0.5 } 818 }; 819 PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 820 } 821 { 822 const PetscReal A[3][3] = { 823 {0, 0, 0}, 824 {2 - s2, 0, 0}, 825 {0.5, 0.5, 0} 826 }; 827 const PetscReal At[3][3] = { 828 {0, 0, 0 }, 829 {1 - 1 / s2, 1 - 1 / s2, 0 }, 830 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 831 }; 832 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 833 const PetscReal binterpt[3][2] = { 834 {1.0 / s2, -1.0 / (2.0 * s2)}, 835 {1.0 / s2, -1.0 / (2.0 * s2)}, 836 {1.0 - s2, 1.0 / s2 } 837 }; 838 PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 839 } 840 { 841 const PetscReal A[3][3] = { 842 {0, 0, 0}, 843 {2 - s2, 0, 0}, 844 {0.75, 0.25, 0} 845 }; 846 const PetscReal At[3][3] = { 847 {0, 0, 0 }, 848 {1 - 1 / s2, 1 - 1 / s2, 0 }, 849 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 850 }; 851 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 852 const PetscReal binterpt[3][2] = { 853 {1.0 / s2, -1.0 / (2.0 * s2)}, 854 {1.0 / s2, -1.0 / (2.0 * s2)}, 855 {1.0 - s2, 1.0 / s2 } 856 }; 857 PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 858 } 859 { /* Optimal for linear implicit part */ 860 const PetscReal A[3][3] = { 861 {0, 0, 0}, 862 {2 - s2, 0, 0}, 863 {(3 - 2 * s2) / 6, (3 + 2 * s2) / 6, 0} 864 }; 865 const PetscReal At[3][3] = { 866 {0, 0, 0 }, 867 {1 - 1 / s2, 1 - 1 / s2, 0 }, 868 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 869 }; 870 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 871 const PetscReal binterpt[3][2] = { 872 {1.0 / s2, -1.0 / (2.0 * s2)}, 873 {1.0 / s2, -1.0 / (2.0 * s2)}, 874 {1.0 - s2, 1.0 / s2 } 875 }; 876 PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 877 } 878 { /* Optimal for linear implicit part */ 879 const PetscReal A[3][3] = { 880 {0, 0, 0}, 881 {0.5, 0, 0}, 882 {0.5, 0.5, 0} 883 }; 884 const PetscReal At[3][3] = { 885 {0.25, 0, 0 }, 886 {0, 0.25, 0 }, 887 {1. / 3, 1. / 3, 1. / 3} 888 }; 889 PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 890 } 891 { 892 const PetscReal A[4][4] = { 893 {0, 0, 0, 0}, 894 {1767732205903. / 2027836641118., 0, 0, 0}, 895 {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 896 {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 897 }; 898 const PetscReal At[4][4] = { 899 {0, 0, 0, 0 }, 900 {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0 }, 901 {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0 }, 902 {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.} 903 }; 904 const PetscReal bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}; 905 const PetscReal binterpt[4][2] = { 906 {4655552711362. / 22874653954995., -215264564351. / 13552729205753. }, 907 {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119. }, 908 {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, 909 {584795268549. / 6622622206610., 2508943948391. / 7218656332882. } 910 }; 911 PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 912 } 913 { 914 const PetscReal A[5][5] = { 915 {0, 0, 0, 0, 0}, 916 {1. / 2, 0, 0, 0, 0}, 917 {11. / 18, 1. / 18, 0, 0, 0}, 918 {5. / 6, -5. / 6, .5, 0, 0}, 919 {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 920 }; 921 const PetscReal At[5][5] = { 922 {0, 0, 0, 0, 0 }, 923 {0, 1. / 2, 0, 0, 0 }, 924 {0, 1. / 6, 1. / 2, 0, 0 }, 925 {0, -1. / 2, 1. / 2, 1. / 2, 0 }, 926 {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2} 927 }; 928 PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 929 } 930 { 931 const PetscReal A[5][5] = { 932 {0, 0, 0, 0, 0}, 933 {1, 0, 0, 0, 0}, 934 {4. / 9, 2. / 9, 0, 0, 0}, 935 {1. / 4, 0, 3. / 4, 0, 0}, 936 {1. / 4, 0, 3. / 5, 0, 0} 937 }; 938 const PetscReal At[5][5] = { 939 {0, 0, 0, 0, 0 }, 940 {.5, .5, 0, 0, 0 }, 941 {5. / 18, -1. / 9, .5, 0, 0 }, 942 {.5, 0, 0, .5, 0 }, 943 {.25, 0, .75, -.5, .5} 944 }; 945 PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 946 } 947 { 948 const PetscReal A[6][6] = { 949 {0, 0, 0, 0, 0, 0}, 950 {1. / 2, 0, 0, 0, 0, 0}, 951 {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 952 {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 953 {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 954 {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 955 }; 956 const PetscReal At[6][6] = { 957 {0, 0, 0, 0, 0, 0 }, 958 {1. / 4, 1. / 4, 0, 0, 0, 0 }, 959 {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0 }, 960 {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0 }, 961 {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0 }, 962 {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4} 963 }; 964 const PetscReal bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}; 965 const PetscReal binterpt[6][3] = { 966 {6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025. }, 967 {0, 0, 0 }, 968 {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035. }, 969 {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 970 {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469. }, 971 {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.} 972 }; 973 PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 974 } 975 { 976 const PetscReal A[8][8] = { 977 {0, 0, 0, 0, 0, 0, 0, 0}, 978 {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 979 {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 980 {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 981 {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 982 {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 983 {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 984 {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 985 }; 986 const PetscReal At[8][8] = { 987 {0, 0, 0, 0, 0, 0, 0, 0 }, 988 {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0 }, 989 {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0 }, 990 {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0 }, 991 {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0 }, 992 {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0 }, 993 {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0 }, 994 {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.} 995 }; 996 const PetscReal bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}; 997 const PetscReal binterpt[8][3] = { 998 {-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439. }, 999 {0, 0, 0 }, 1000 {0, 0, 0 }, 1001 {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, 1002 {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, 1003 {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, 1004 {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, 1005 {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460. } 1006 }; 1007 PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 1008 } 1009 #undef RC 1010 #undef us2 1011 #undef s2 1012 PetscFunctionReturn(PETSC_SUCCESS); 1013 } 1014 1015 /*@C 1016 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`. 1017 1018 Not Collective 1019 1020 Level: advanced 1021 1022 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 1023 @*/ 1024 PetscErrorCode TSARKIMEXRegisterDestroy(void) 1025 { 1026 ARKTableauLink link; 1027 1028 PetscFunctionBegin; 1029 while ((link = ARKTableauList)) { 1030 ARKTableau t = &link->tab; 1031 ARKTableauList = link->next; 1032 PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 1033 PetscCall(PetscFree2(t->bembedt, t->bembed)); 1034 PetscCall(PetscFree2(t->binterpt, t->binterp)); 1035 PetscCall(PetscFree(t->name)); 1036 PetscCall(PetscFree(link)); 1037 } 1038 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } 1041 1042 /*@C 1043 TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called 1044 from `TSInitializePackage()`. 1045 1046 Level: developer 1047 1048 .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()` 1049 @*/ 1050 PetscErrorCode TSARKIMEXInitializePackage(void) 1051 { 1052 PetscFunctionBegin; 1053 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 1054 TSARKIMEXPackageInitialized = PETSC_TRUE; 1055 PetscCall(TSARKIMEXRegisterAll()); 1056 PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 1057 PetscFunctionReturn(PETSC_SUCCESS); 1058 } 1059 1060 /*@C 1061 TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is 1062 called from `PetscFinalize()`. 1063 1064 Level: developer 1065 1066 .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()` 1067 @*/ 1068 PetscErrorCode TSARKIMEXFinalizePackage(void) 1069 { 1070 PetscFunctionBegin; 1071 TSARKIMEXPackageInitialized = PETSC_FALSE; 1072 PetscCall(TSARKIMEXRegisterDestroy()); 1073 PetscFunctionReturn(PETSC_SUCCESS); 1074 } 1075 1076 /*@C 1077 TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 1078 1079 Logically Collective 1080 1081 Input Parameters: 1082 + name - identifier for method 1083 . order - approximation order of method 1084 . s - number of stages, this is the dimension of the matrices below 1085 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 1086 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 1087 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 1088 . A - Non-stiff stage coefficients (dimension s*s, row-major) 1089 . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 1090 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 1091 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 1092 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 1093 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 1094 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 1095 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 1096 1097 Level: advanced 1098 1099 Note: 1100 Several `TSARKIMEX` methods are provided, this function is only needed to create new methods. 1101 1102 .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS` 1103 @*/ 1104 PetscErrorCode TSARKIMEXRegister(TSARKIMEXType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal A[], const PetscReal b[], const PetscReal c[], const PetscReal bembedt[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[], const PetscReal binterp[]) 1105 { 1106 ARKTableauLink link; 1107 ARKTableau t; 1108 PetscInt i, j; 1109 1110 PetscFunctionBegin; 1111 PetscCheck(s > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected number of stages s %" PetscInt_FMT " > 0", s); 1112 PetscCall(TSARKIMEXInitializePackage()); 1113 for (link = ARKTableauList; link; link = link->next) { 1114 PetscBool match; 1115 1116 PetscCall(PetscStrcmp(link->tab.name, name, &match)); 1117 PetscCheck(!match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Method with name \"%s\" already registered", name); 1118 } 1119 PetscCall(PetscNew(&link)); 1120 t = &link->tab; 1121 PetscCall(PetscStrallocpy(name, &t->name)); 1122 t->order = order; 1123 t->s = s; 1124 PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 1125 PetscCall(PetscArraycpy(t->At, At, s * s)); 1126 if (A) { 1127 PetscCall(PetscArraycpy(t->A, A, s * s)); 1128 t->additive = PETSC_TRUE; 1129 } 1130 1131 if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 1132 else 1133 for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 1134 1135 if (t->additive) { 1136 if (b) PetscCall(PetscArraycpy(t->b, b, s)); 1137 else 1138 for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 1139 } else PetscCall(PetscArrayzero(t->b, s)); 1140 1141 if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 1142 else 1143 for (i = 0; i < s; i++) 1144 for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 1145 1146 if (t->additive) { 1147 if (c) PetscCall(PetscArraycpy(t->c, c, s)); 1148 else 1149 for (i = 0; i < s; i++) 1150 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 1151 } else PetscCall(PetscArrayzero(t->c, s)); 1152 1153 t->stiffly_accurate = PETSC_TRUE; 1154 for (i = 0; i < s; i++) 1155 if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 1156 1157 t->explicit_first_stage = PETSC_TRUE; 1158 for (i = 0; i < s; i++) 1159 if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 1160 1161 /* def of FSAL can be made more precise */ 1162 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 1163 1164 if (bembedt) { 1165 PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 1166 PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 1167 PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 1168 } 1169 1170 t->pinterp = pinterp; 1171 PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 1172 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 1173 PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 1174 1175 link->next = ARKTableauList; 1176 ARKTableauList = link; 1177 PetscFunctionReturn(PETSC_SUCCESS); 1178 } 1179 1180 /*@C 1181 TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation 1182 1183 Logically Collective. 1184 1185 Input Parameters: 1186 + name - identifier for method 1187 . order - approximation order of method 1188 . s - number of stages, this is the dimension of the matrices below 1189 . At - Butcher table of stage coefficients (dimension `s`*`s`, row-major order) 1190 . bt - Butcher table for completing the step (dimension `s`; pass `NULL` to use the last row of `At`) 1191 . ct - Abscissa of each stage (dimension s, NULL to use row sums of At) 1192 . bembedt - Stiff part of completion table for embedded method (dimension s; `NULL` if not available) 1193 . pinterp - Order of the interpolation scheme, equal to the number of columns of `binterpt` and `binterp` 1194 - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 1195 1196 Level: advanced 1197 1198 Note: 1199 Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods. 1200 1201 .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS` 1202 @*/ 1203 PetscErrorCode TSDIRKRegister(TSDIRKType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal bembedt[], PetscInt pinterp, const PetscReal binterpt[]) 1204 { 1205 PetscFunctionBegin; 1206 PetscCall(TSARKIMEXRegister(name, order, s, At, bt, ct, NULL, NULL, NULL, bembedt, NULL, pinterp, binterpt, NULL)); 1207 PetscFunctionReturn(PETSC_SUCCESS); 1208 } 1209 1210 /* 1211 The step completion formula is 1212 1213 x1 = x0 - h bt^T YdotI + h b^T YdotRHS 1214 1215 This function can be called before or after ts->vec_sol has been updated. 1216 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 1217 We can write 1218 1219 x1e = x0 - h bet^T YdotI + h be^T YdotRHS 1220 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 1221 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 1222 1223 so we can evaluate the method with different order even after the step has been optimistically completed. 1224 */ 1225 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) 1226 { 1227 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1228 ARKTableau tab = ark->tableau; 1229 PetscScalar *w = ark->work; 1230 PetscReal h; 1231 PetscInt s = tab->s, j; 1232 PetscBool hasE; 1233 1234 PetscFunctionBegin; 1235 switch (ark->status) { 1236 case TS_STEP_INCOMPLETE: 1237 case TS_STEP_PENDING: 1238 h = ts->time_step; 1239 break; 1240 case TS_STEP_COMPLETE: 1241 h = ts->ptime - ts->ptime_prev; 1242 break; 1243 default: 1244 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1245 } 1246 if (order == tab->order) { 1247 if (ark->status == TS_STEP_INCOMPLETE) { 1248 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 1249 PetscCall(VecCopy(ark->Y[s - 1], X)); 1250 } else { /* Use the standard completion formula (bt,b) */ 1251 PetscCall(VecCopy(ts->vec_sol, X)); 1252 for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 1253 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1254 if (tab->additive && ark->imex) { /* Method is IMEX, complete the explicit formula */ 1255 PetscCall(TSHasRHSFunction(ts, &hasE)); 1256 if (hasE) { 1257 for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 1258 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1259 } 1260 } 1261 } 1262 } else PetscCall(VecCopy(ts->vec_sol, X)); 1263 if (done) *done = PETSC_TRUE; 1264 PetscFunctionReturn(PETSC_SUCCESS); 1265 } else if (order == tab->order - 1) { 1266 if (!tab->bembedt) goto unavailable; 1267 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 1268 PetscCall(VecCopy(ts->vec_sol, X)); 1269 for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 1270 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1271 if (tab->additive) { 1272 PetscCall(TSHasRHSFunction(ts, &hasE)); 1273 if (hasE) { 1274 for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 1275 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1276 } 1277 } 1278 } else { /* Rollback and re-complete using (bet-be,be-b) */ 1279 PetscCall(VecCopy(ts->vec_sol, X)); 1280 for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 1281 PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 1282 if (tab->additive) { 1283 PetscCall(TSHasRHSFunction(ts, &hasE)); 1284 if (hasE) { 1285 for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 1286 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1287 } 1288 } 1289 } 1290 if (done) *done = PETSC_TRUE; 1291 PetscFunctionReturn(PETSC_SUCCESS); 1292 } 1293 unavailable: 1294 PetscCheck(done, 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.", 1295 tab->name, tab->order, order); 1296 *done = PETSC_FALSE; 1297 PetscFunctionReturn(PETSC_SUCCESS); 1298 } 1299 1300 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 1301 { 1302 Vec Udot, Y1, Y2; 1303 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1304 PetscReal norm; 1305 1306 PetscFunctionBegin; 1307 PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 1308 PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 1309 PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 1310 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 1311 PetscCall(VecSetRandom(Udot, NULL)); 1312 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 1313 PetscCall(VecAXPY(Y2, -1.0, Y1)); 1314 PetscCall(VecAXPY(Y2, -1.0, Udot)); 1315 PetscCall(VecNorm(Y2, NORM_2, &norm)); 1316 if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 1317 *id = PETSC_TRUE; 1318 } else { 1319 *id = PETSC_FALSE; 1320 PetscCall(PetscInfo(ts, "IFunction(Udot = random) - IFunction(Udot = 0) is not near Udot, %g, suspect mass matrix implied in IFunction() is not the identity as required\n", (double)norm)); 1321 } 1322 PetscCall(VecDestroy(&Udot)); 1323 PetscCall(VecDestroy(&Y1)); 1324 PetscCall(VecDestroy(&Y2)); 1325 PetscFunctionReturn(PETSC_SUCCESS); 1326 } 1327 1328 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS, PetscReal, Vec, IS *); 1329 1330 static PetscErrorCode TSStep_ARKIMEX(TS ts) 1331 { 1332 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1333 ARKTableau tab = ark->tableau; 1334 const PetscInt s = tab->s; 1335 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 1336 PetscScalar *w = ark->work; 1337 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 1338 PetscBool extrapolate = ark->extrapolate; 1339 TSAdapt adapt; 1340 SNES snes; 1341 PetscInt i, j, its, lits; 1342 PetscInt rejections = 0; 1343 PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE; 1344 PetscReal next_time_step = ts->time_step; 1345 1346 PetscFunctionBegin; 1347 if (ark->extrapolate && !ark->Y_prev) { 1348 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1349 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1350 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1351 } 1352 1353 if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE)); 1354 if (!hasE) dirk = PETSC_TRUE; 1355 1356 if (!ts->steprollback) { 1357 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 1358 PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 1359 } 1360 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 1361 for (i = 0; i < s; i++) { 1362 PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 1363 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 1364 if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 1365 } 1366 } 1367 } 1368 1369 /* 1370 For fully implicit formulations we must solve the equations 1371 1372 F(t_n,x_n,xdot) = 0 1373 1374 for the explicit first stage. 1375 Here we call SNESSolve using PETSC_MAX_REAL as shift to flag it. 1376 Special handling is inside SNESTSFormFunction_ARKIMEX and SNESTSFormJacobian_ARKIMEX 1377 We compute Ydot0 if we restart the step or if we resized the problem after remeshing 1378 */ 1379 if (dirk && tab->explicit_first_stage && (ts->steprestart || ts->stepresize)) { 1380 ark->scoeff = PETSC_MAX_REAL; 1381 PetscCall(VecCopy(ts->vec_sol, Z)); 1382 if (!ark->alg_is) { 1383 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ts->ptime, Z, &ark->alg_is)); 1384 PetscCall(ISViewFromOptions(ark->alg_is, (PetscObject)ts, "-ts_arkimex_algebraic_is_view")); 1385 } 1386 PetscCall(TSGetSNES(ts, &snes)); 1387 PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, 1)); 1388 PetscCall(SNESSolve(snes, NULL, Ydot0)); 1389 if (ark->alg_is) PetscCall(VecISSet(Ydot0, ark->alg_is, 0.0)); 1390 PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, -1)); 1391 } 1392 1393 /* For IMEX we compute a step */ 1394 if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 1395 TS ts_start; 1396 if (PetscDefined(USE_DEBUG) && hasE) { 1397 PetscBool id = PETSC_FALSE; 1398 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1399 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "This scheme requires an identity mass matrix, however the TSIFunctionFn you provided does not utilize an identity mass matrix"); 1400 } 1401 PetscCall(TSClone(ts, &ts_start)); 1402 PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 1403 PetscCall(TSSetTime(ts_start, ts->ptime)); 1404 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 1405 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 1406 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 1407 PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 1408 PetscCall(TSSetType(ts_start, TSARKIMEX)); 1409 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 1410 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 1411 1412 PetscCall(TSRestartStep(ts_start)); 1413 PetscCall(TSSolve(ts_start, ts->vec_sol)); 1414 PetscCall(TSGetTime(ts_start, &ts->ptime)); 1415 PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 1416 1417 { /* Save the initial slope for the next step */ 1418 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 1419 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 1420 } 1421 ts->steps++; 1422 1423 /* Set the correct TS in SNES */ 1424 /* We'll try to bypass this by changing the method on the fly */ 1425 { 1426 PetscCall(TSGetSNES(ts, &snes)); 1427 PetscCall(TSSetSNES(ts, snes)); 1428 } 1429 PetscCall(TSDestroy(&ts_start)); 1430 } 1431 1432 ark->status = TS_STEP_INCOMPLETE; 1433 while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 1434 PetscReal t = ts->ptime; 1435 PetscReal h = ts->time_step; 1436 for (i = 0; i < s; i++) { 1437 ark->stage_time = t + h * ct[i]; 1438 PetscCall(TSPreStage(ts, ark->stage_time)); 1439 if (At[i * s + i] == 0) { /* This stage is explicit */ 1440 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"); 1441 PetscCall(VecCopy(ts->vec_sol, Y[i])); 1442 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1443 PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 1444 if (tab->additive && hasE) { 1445 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1446 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 1447 } 1448 PetscCall(TSGetSNES(ts, &snes)); 1449 PetscCall(SNESResetCounters(snes)); 1450 } else { 1451 ark->scoeff = 1. / At[i * s + i]; 1452 /* Ydot = shift*(Y-Z) */ 1453 PetscCall(VecCopy(ts->vec_sol, Z)); 1454 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1455 PetscCall(VecMAXPY(Z, i, w, YdotI)); 1456 if (tab->additive && hasE) { 1457 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1458 PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 1459 } 1460 if (extrapolate && !ts->steprestart) { 1461 /* Initial guess extrapolated from previous time step stage values */ 1462 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 1463 } else { 1464 /* Initial guess taken from last stage */ 1465 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 1466 } 1467 PetscCall(TSGetSNES(ts, &snes)); 1468 PetscCall(SNESSolve(snes, NULL, Y[i])); 1469 PetscCall(SNESGetIterationNumber(snes, &its)); 1470 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 1471 ts->snes_its += its; 1472 ts->ksp_its += lits; 1473 PetscCall(TSGetAdapt(ts, &adapt)); 1474 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 1475 if (!stageok) { 1476 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 1477 * use extrapolation to initialize the solves on the next attempt. */ 1478 extrapolate = PETSC_FALSE; 1479 goto reject_step; 1480 } 1481 } 1482 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { 1483 if (i == 0 && tab->explicit_first_stage) { 1484 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", 1485 ((PetscObject)ts)->type_name, ark->tableau->name); 1486 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 1487 } else { 1488 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1489 } 1490 } else { 1491 if (i == 0 && tab->explicit_first_stage) { 1492 PetscCall(VecZeroEntries(Ydot)); 1493 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 1494 PetscCall(VecScale(YdotI[i], -1.0)); 1495 } else { 1496 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1497 } 1498 if (hasE) { 1499 if (ark->imex) { 1500 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 1501 } else { 1502 PetscCall(VecZeroEntries(YdotRHS[i])); 1503 } 1504 } 1505 } 1506 PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 1507 } 1508 1509 ark->status = TS_STEP_INCOMPLETE; 1510 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 1511 ark->status = TS_STEP_PENDING; 1512 PetscCall(TSGetAdapt(ts, &adapt)); 1513 PetscCall(TSAdaptCandidatesClear(adapt)); 1514 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 1515 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1516 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1517 if (!accept) { /* Roll back the current step */ 1518 PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1519 ts->time_step = next_time_step; 1520 goto reject_step; 1521 } 1522 1523 ts->ptime += ts->time_step; 1524 ts->time_step = next_time_step; 1525 break; 1526 1527 reject_step: 1528 ts->reject++; 1529 accept = PETSC_FALSE; 1530 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1531 ts->reason = TS_DIVERGED_STEP_REJECTED; 1532 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1533 } 1534 } 1535 PetscFunctionReturn(PETSC_SUCCESS); 1536 } 1537 1538 /* 1539 This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 1540 Udot = H(t,U) + G(t,U) 1541 This corresponds to F(t,U,Udot) = Udot-H(t,U). 1542 1543 The complete adjoint equations are 1544 (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 1545 dGdU (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 1546 + dHdU (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 1547 lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 1548 mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 1549 + dGdP (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 1550 + dHdP (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 1551 mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 1552 where shift = 1/(at[i][i]*h) 1553 1554 If at[i][i] is 0, the first equation falls back to 1555 lambda_s[i] = h ( 1556 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 1557 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 1558 1559 */ 1560 static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 1561 { 1562 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1563 ARKTableau tab = ark->tableau; 1564 const PetscInt s = tab->s; 1565 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 1566 PetscScalar *w = ark->work; 1567 Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 1568 Mat Jex, Jim, Jimpre; 1569 PetscInt i, j, nadj; 1570 PetscReal t = ts->ptime, stage_time_ex; 1571 PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 1572 KSP ksp; 1573 1574 PetscFunctionBegin; 1575 ark->status = TS_STEP_INCOMPLETE; 1576 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1577 PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 1578 PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 1579 1580 for (i = s - 1; i >= 0; i--) { 1581 ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 1582 stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 1583 if (At[i * s + i] == 0) { // This stage is explicit 1584 ark->scoeff = 0.; 1585 } else { 1586 ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 1587 } 1588 PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 1589 PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 1590 if (ts->vecs_sensip) { 1591 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 1592 PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 1593 } 1594 /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 1595 for (nadj = 0; nadj < ts->numcost; nadj++) { 1596 /* build implicit part */ 1597 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1598 if (s - i - 1 > 0) { 1599 /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 1600 for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 1601 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1602 } 1603 /* Temp = Temp - bt[i] lambda_{n+1} */ 1604 PetscCall(VecAXPY(VecsSensiTemp[nadj], -bt[i], ts->vecs_sensi[nadj])); 1605 if (bt[i] || s - i - 1 > 0) { 1606 /* (shift I - dHdU) Temp */ 1607 PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 1608 /* cancel out shift Temp where shift=-scoeff/h */ 1609 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 1610 if (ts->vecs_sensip) { 1611 /* - dHdP Temp */ 1612 PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1613 /* mu_n += -h dHdP Temp */ 1614 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1615 } 1616 } else { 1617 PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 1618 } 1619 /* build explicit part */ 1620 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1621 if (s - i - 1 > 0) { 1622 /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 1623 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 1624 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1625 } 1626 /* Temp = Temp + b[i] lambda_{n+1} */ 1627 PetscCall(VecAXPY(VecsSensiTemp[nadj], b[i], ts->vecs_sensi[nadj])); 1628 if (b[i] || s - i - 1 > 0) { 1629 /* dGdU Temp */ 1630 PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1631 if (ts->vecs_sensip) { 1632 /* dGdP Temp */ 1633 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1634 /* mu_n += h dGdP Temp */ 1635 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1636 } 1637 } 1638 /* Build LHS for first-order adjoint */ 1639 if (At[i * s + i] == 0) { // This stage is explicit 1640 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 1641 } else { 1642 KSP ksp; 1643 KSPConvergedReason kspreason; 1644 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1645 PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 1646 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 1647 PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1648 PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 1649 if (kspreason < 0) { 1650 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 1651 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj)); 1652 } 1653 if (ts->vecs_sensip) { 1654 /* -dHdP lambda_s[i] */ 1655 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 1656 /* mu_n += h at[i][i] dHdP lambda_s[i] */ 1657 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 1658 } 1659 } 1660 } 1661 } 1662 for (j = 0; j < s; j++) w[j] = 1.0; 1663 for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 1664 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 1665 ark->status = TS_STEP_COMPLETE; 1666 PetscFunctionReturn(PETSC_SUCCESS); 1667 } 1668 1669 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1670 { 1671 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1672 ARKTableau tab = ark->tableau; 1673 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1674 PetscReal h; 1675 PetscReal tt, t; 1676 PetscScalar *bt = ark->work, *b = ark->work + s; 1677 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1678 1679 PetscFunctionBegin; 1680 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name); 1681 switch (ark->status) { 1682 case TS_STEP_INCOMPLETE: 1683 case TS_STEP_PENDING: 1684 h = ts->time_step; 1685 t = (itime - ts->ptime) / h; 1686 break; 1687 case TS_STEP_COMPLETE: 1688 h = ts->ptime - ts->ptime_prev; 1689 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1690 break; 1691 default: 1692 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1693 } 1694 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1695 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1696 for (i = 0; i < s; i++) { 1697 bt[i] += h * Bt[i * pinterp + j] * tt; 1698 b[i] += h * B[i * pinterp + j] * tt; 1699 } 1700 } 1701 PetscCall(VecCopy(ark->Y[0], X)); 1702 PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1703 if (tab->additive) { 1704 PetscBool hasE; 1705 PetscCall(TSHasRHSFunction(ts, &hasE)); 1706 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1707 } 1708 PetscFunctionReturn(PETSC_SUCCESS); 1709 } 1710 1711 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1712 { 1713 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1714 ARKTableau tab = ark->tableau; 1715 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1716 PetscReal h, h_prev, t, tt; 1717 PetscScalar *bt = ark->work, *b = ark->work + s; 1718 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1719 1720 PetscFunctionBegin; 1721 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 1722 h = ts->time_step; 1723 h_prev = ts->ptime - ts->ptime_prev; 1724 t = 1 + h / h_prev * c; 1725 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1726 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1727 for (i = 0; i < s; i++) { 1728 bt[i] += h * Bt[i * pinterp + j] * tt; 1729 b[i] += h * B[i * pinterp + j] * tt; 1730 } 1731 } 1732 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 1733 PetscCall(VecCopy(ark->Y_prev[0], X)); 1734 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1735 if (tab->additive) { 1736 PetscBool hasE; 1737 PetscCall(TSHasRHSFunction(ts, &hasE)); 1738 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1739 } 1740 PetscFunctionReturn(PETSC_SUCCESS); 1741 } 1742 1743 static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1744 { 1745 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1746 ARKTableau tab = ark->tableau; 1747 1748 PetscFunctionBegin; 1749 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1750 PetscCall(PetscFree(ark->work)); 1751 PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 1752 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 1753 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 1754 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 1755 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 1756 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 1757 PetscFunctionReturn(PETSC_SUCCESS); 1758 } 1759 1760 static PetscErrorCode TSReset_ARKIMEX(TS ts) 1761 { 1762 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1763 1764 PetscFunctionBegin; 1765 if (ark->fastslowsplit) { 1766 PetscTryMethod(ts, "TSReset_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 1767 } else { 1768 PetscCall(TSARKIMEXTableauReset(ts)); 1769 PetscCall(VecDestroy(&ark->Ydot)); 1770 PetscCall(VecDestroy(&ark->Ydot0)); 1771 PetscCall(VecDestroy(&ark->Z)); 1772 PetscCall(ISDestroy(&ark->alg_is)); 1773 } 1774 PetscFunctionReturn(PETSC_SUCCESS); 1775 } 1776 1777 static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 1778 { 1779 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1780 ARKTableau tab = ark->tableau; 1781 1782 PetscFunctionBegin; 1783 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 1784 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 1785 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 1786 PetscFunctionReturn(PETSC_SUCCESS); 1787 } 1788 1789 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1790 { 1791 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1792 1793 PetscFunctionBegin; 1794 if (Z) { 1795 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1796 else *Z = ax->Z; 1797 } 1798 if (Ydot) { 1799 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1800 else *Ydot = ax->Ydot; 1801 } 1802 PetscFunctionReturn(PETSC_SUCCESS); 1803 } 1804 1805 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1806 { 1807 PetscFunctionBegin; 1808 if (Z) { 1809 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1810 } 1811 if (Ydot) { 1812 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1813 } 1814 PetscFunctionReturn(PETSC_SUCCESS); 1815 } 1816 1817 /* 1818 DAEs need special handling for algebraic variables when restarting DIRK methods with explicit 1819 first stage. In particular, we need: 1820 - to zero the nonlinear function (in case the dual variables are not consistent in the first step) 1821 - to modify the matrix by calling MatZeroRows with identity on these variables. 1822 */ 1823 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is) 1824 { 1825 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1826 DM dm; 1827 Vec F, W, Xdot; 1828 const PetscScalar *w; 1829 PetscInt nz = 0, n, st; 1830 PetscInt *nzr; 1831 1832 PetscFunctionBegin; 1833 PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */ 1834 PetscCall(DMGetGlobalVector(dm, &Xdot)); 1835 PetscCall(DMGetGlobalVector(dm, &F)); 1836 PetscCall(DMGetGlobalVector(dm, &W)); 1837 PetscCall(VecSet(Xdot, 0.0)); 1838 PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex)); 1839 PetscCall(VecSetRandom(Xdot, NULL)); 1840 PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex)); 1841 PetscCall(VecAXPY(W, -1.0, F)); 1842 PetscCall(VecGetOwnershipRange(W, &st, NULL)); 1843 PetscCall(VecGetLocalSize(W, &n)); 1844 PetscCall(VecGetArrayRead(W, &w)); 1845 for (PetscInt i = 0; i < n; i++) 1846 if (w[i] == 0.0) nz++; 1847 PetscCall(PetscMalloc1(nz, &nzr)); 1848 nz = 0; 1849 for (PetscInt i = 0; i < n; i++) 1850 if (w[i] == 0.0) nzr[nz++] = i + st; 1851 PetscCall(VecRestoreArrayRead(W, &w)); 1852 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is)); 1853 PetscCall(DMRestoreGlobalVector(dm, &Xdot)); 1854 PetscCall(DMRestoreGlobalVector(dm, &F)); 1855 PetscCall(DMRestoreGlobalVector(dm, &W)); 1856 PetscFunctionReturn(PETSC_SUCCESS); 1857 } 1858 1859 /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure 1860 at the finest level, in the DM for coarser solves. */ 1861 static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is) 1862 { 1863 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1864 1865 PetscFunctionBegin; 1866 if (dm && dm != ts->dm) PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is)); 1867 else *alg_is = ax->alg_is; 1868 PetscFunctionReturn(PETSC_SUCCESS); 1869 } 1870 1871 /* This defines the nonlinear equation that is to be solved with SNES */ 1872 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1873 { 1874 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1875 DM dm, dmsave; 1876 Vec Z, Ydot; 1877 IS alg_is; 1878 1879 PetscFunctionBegin; 1880 PetscCall(SNESGetDM(snes, &dm)); 1881 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1882 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 1883 1884 dmsave = ts->dm; 1885 ts->dm = dm; 1886 1887 if (ark->scoeff == PETSC_MAX_REAL) { 1888 /* We are solving F(t_n,x_n,xdot) = 0 to start the method */ 1889 if (!alg_is) { 1890 PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 1891 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is)); 1892 PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is)); 1893 PetscCall(PetscObjectDereference((PetscObject)alg_is)); 1894 PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view")); 1895 } 1896 PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 1897 PetscCall(VecISSet(F, alg_is, 0.0)); 1898 } else { 1899 PetscReal shift = ark->scoeff / ts->time_step; 1900 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1901 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1902 } 1903 1904 ts->dm = dmsave; 1905 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1906 PetscFunctionReturn(PETSC_SUCCESS); 1907 } 1908 1909 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1910 { 1911 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1912 DM dm, dmsave; 1913 Vec Ydot, Z; 1914 PetscReal shift; 1915 IS alg_is; 1916 1917 PetscFunctionBegin; 1918 PetscCall(SNESGetDM(snes, &dm)); 1919 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1920 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1921 /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */ 1922 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 1923 1924 dmsave = ts->dm; 1925 ts->dm = dm; 1926 1927 if (ark->scoeff == PETSC_MAX_REAL) { 1928 PetscBool hasZeroRows; 1929 1930 /* We are solving F(t_n,x_n,xdot) = 0 to start the method 1931 We compute with a very large shift and then scale back the matrix */ 1932 shift = 1.0 / PETSC_MACHINE_EPSILON; 1933 PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1934 PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 1935 PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows)); 1936 if (hasZeroRows) { 1937 PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 1938 /* the default of AIJ is to not keep the pattern! We should probably change it someday */ 1939 PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 1940 PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL)); 1941 } 1942 PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view")); 1943 if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1944 } else { 1945 shift = ark->scoeff / ts->time_step; 1946 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1947 } 1948 ts->dm = dmsave; 1949 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1950 PetscFunctionReturn(PETSC_SUCCESS); 1951 } 1952 1953 static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 1954 { 1955 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1956 1957 PetscFunctionBegin; 1958 if (ns) *ns = ark->tableau->s; 1959 if (Y) *Y = ark->Y; 1960 PetscFunctionReturn(PETSC_SUCCESS); 1961 } 1962 1963 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1964 { 1965 PetscFunctionBegin; 1966 PetscFunctionReturn(PETSC_SUCCESS); 1967 } 1968 1969 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1970 { 1971 TS ts = (TS)ctx; 1972 Vec Z, Z_c; 1973 1974 PetscFunctionBegin; 1975 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 1976 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 1977 PetscCall(MatRestrict(restrct, Z, Z_c)); 1978 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 1979 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 1980 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 1981 PetscFunctionReturn(PETSC_SUCCESS); 1982 } 1983 1984 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1985 { 1986 PetscFunctionBegin; 1987 PetscFunctionReturn(PETSC_SUCCESS); 1988 } 1989 1990 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1991 { 1992 TS ts = (TS)ctx; 1993 Vec Z, Z_c; 1994 1995 PetscFunctionBegin; 1996 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 1997 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 1998 1999 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2000 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2001 2002 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 2003 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 2004 PetscFunctionReturn(PETSC_SUCCESS); 2005 } 2006 2007 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 2008 { 2009 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2010 ARKTableau tab = ark->tableau; 2011 2012 PetscFunctionBegin; 2013 PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 2014 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 2015 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 2016 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 2017 if (ark->extrapolate) { 2018 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 2019 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 2020 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 2021 } 2022 PetscFunctionReturn(PETSC_SUCCESS); 2023 } 2024 2025 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 2026 { 2027 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2028 DM dm; 2029 SNES snes; 2030 2031 PetscFunctionBegin; 2032 if (ark->fastslowsplit) { 2033 PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 2034 } else { 2035 PetscCall(TSARKIMEXTableauSetUp(ts)); 2036 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 2037 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 2038 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 2039 PetscCall(TSGetDM(ts, &dm)); 2040 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 2041 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2042 PetscCall(TSGetSNES(ts, &snes)); 2043 PetscCall(SNESSetDM(snes, dm)); 2044 } 2045 PetscFunctionReturn(PETSC_SUCCESS); 2046 } 2047 2048 static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 2049 { 2050 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2051 ARKTableau tab = ark->tableau; 2052 2053 PetscFunctionBegin; 2054 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 2055 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 2056 if (ts->vecs_sensip) PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); 2057 if (PetscDefined(USE_DEBUG)) { 2058 PetscBool id = PETSC_FALSE; 2059 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 2060 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Adjoint ARKIMEX requires an identity mass matrix, however the TSIFunctionFn you provided does not utilize an identity mass matrix"); 2061 } 2062 PetscFunctionReturn(PETSC_SUCCESS); 2063 } 2064 2065 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems PetscOptionsObject) 2066 { 2067 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2068 PetscBool dirk; 2069 2070 PetscFunctionBegin; 2071 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2072 PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 2073 { 2074 ARKTableauLink link; 2075 PetscInt count, choice; 2076 PetscBool flg; 2077 const char **namelist; 2078 for (link = ARKTableauList, count = 0; link; link = link->next) { 2079 if (!dirk && link->tab.additive) count++; 2080 if (dirk && !link->tab.additive) count++; 2081 } 2082 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 2083 for (link = ARKTableauList, count = 0; link; link = link->next) { 2084 if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 2085 if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 2086 } 2087 if (dirk) { 2088 PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2089 if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 2090 } else { 2091 PetscBool fastslowsplit; 2092 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2093 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 2094 flg = (PetscBool)!ark->imex; 2095 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 2096 ark->imex = (PetscBool)!flg; 2097 PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg)); 2098 if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit)); 2099 PetscCall(TSARKIMEXGetFastSlowSplit(ts, &fastslowsplit)); 2100 if (fastslowsplit) { 2101 SNES snes; 2102 2103 PetscCall(TSRHSSplitGetSNES(ts, &snes)); 2104 PetscCall(SNESSetFromOptions(snes)); 2105 } 2106 } 2107 PetscCall(PetscFree(namelist)); 2108 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)); 2109 } 2110 PetscOptionsHeadEnd(); 2111 PetscFunctionReturn(PETSC_SUCCESS); 2112 } 2113 2114 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 2115 { 2116 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2117 PetscBool isascii, dirk; 2118 2119 PetscFunctionBegin; 2120 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2121 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 2122 if (isascii) { 2123 PetscViewerFormat format; 2124 ARKTableau tab = ark->tableau; 2125 TSARKIMEXType arktype; 2126 char buf[2048]; 2127 PetscBool flg; 2128 2129 PetscCall(TSARKIMEXGetType(ts, &arktype)); 2130 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 2131 PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 2132 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 2133 PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 2134 PetscCall(PetscViewerGetFormat(viewer, &format)); 2135 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2136 PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 2137 for (PetscInt i = 0; i < tab->s; i++) { 2138 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 2139 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 2140 } 2141 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 2142 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 2143 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 2144 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 2145 } 2146 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 2147 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 2148 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 2149 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 2150 if (!dirk) { 2151 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 2152 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 2153 } 2154 } 2155 PetscFunctionReturn(PETSC_SUCCESS); 2156 } 2157 2158 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 2159 { 2160 SNES snes; 2161 TSAdapt adapt; 2162 2163 PetscFunctionBegin; 2164 PetscCall(TSGetAdapt(ts, &adapt)); 2165 PetscCall(TSAdaptLoad(adapt, viewer)); 2166 PetscCall(TSGetSNES(ts, &snes)); 2167 PetscCall(SNESLoad(snes, viewer)); 2168 /* function and Jacobian context for SNES when used with TS is always ts object */ 2169 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 2170 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 2171 PetscFunctionReturn(PETSC_SUCCESS); 2172 } 2173 2174 /*@ 2175 TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 2176 2177 Logically Collective 2178 2179 Input Parameters: 2180 + ts - timestepping context 2181 - arktype - type of `TSARKIMEX` scheme 2182 2183 Options Database Key: 2184 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 2185 2186 Level: intermediate 2187 2188 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, 2189 `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 2190 @*/ 2191 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 2192 { 2193 PetscFunctionBegin; 2194 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2195 PetscAssertPointer(arktype, 2); 2196 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 2197 PetscFunctionReturn(PETSC_SUCCESS); 2198 } 2199 2200 /*@ 2201 TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 2202 2203 Logically Collective 2204 2205 Input Parameter: 2206 . ts - timestepping context 2207 2208 Output Parameter: 2209 . arktype - type of `TSARKIMEX` scheme 2210 2211 Level: intermediate 2212 2213 .seealso: [](ch_ts), `TSARKIMEX` 2214 @*/ 2215 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2216 { 2217 PetscFunctionBegin; 2218 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2219 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 2220 PetscFunctionReturn(PETSC_SUCCESS); 2221 } 2222 2223 /*@ 2224 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 2225 2226 Logically Collective 2227 2228 Input Parameters: 2229 + ts - timestepping context 2230 - flg - `PETSC_TRUE` for fully implicit 2231 2232 Options Database Key: 2233 . -ts_arkimex_fully_implicit <true,false> - Solve both parts of the equation implicitly 2234 2235 Level: intermediate 2236 2237 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 2238 @*/ 2239 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 2240 { 2241 PetscFunctionBegin; 2242 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2243 PetscValidLogicalCollectiveBool(ts, flg, 2); 2244 PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 2245 PetscFunctionReturn(PETSC_SUCCESS); 2246 } 2247 2248 /*@ 2249 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 2250 2251 Logically Collective 2252 2253 Input Parameter: 2254 . ts - timestepping context 2255 2256 Output Parameter: 2257 . flg - `PETSC_TRUE` for fully implicit 2258 2259 Level: intermediate 2260 2261 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 2262 @*/ 2263 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 2264 { 2265 PetscFunctionBegin; 2266 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2267 PetscAssertPointer(flg, 2); 2268 PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 2269 PetscFunctionReturn(PETSC_SUCCESS); 2270 } 2271 2272 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 2273 { 2274 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2275 2276 PetscFunctionBegin; 2277 *arktype = ark->tableau->name; 2278 PetscFunctionReturn(PETSC_SUCCESS); 2279 } 2280 2281 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 2282 { 2283 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2284 PetscBool match; 2285 ARKTableauLink link; 2286 2287 PetscFunctionBegin; 2288 if (ark->tableau) { 2289 PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 2290 if (match) PetscFunctionReturn(PETSC_SUCCESS); 2291 } 2292 for (link = ARKTableauList; link; link = link->next) { 2293 PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 2294 if (match) { 2295 if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 2296 ark->tableau = &link->tab; 2297 if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 2298 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 2299 PetscFunctionReturn(PETSC_SUCCESS); 2300 } 2301 } 2302 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 2303 } 2304 2305 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 2306 { 2307 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2308 2309 PetscFunctionBegin; 2310 ark->imex = (PetscBool)!flg; 2311 PetscFunctionReturn(PETSC_SUCCESS); 2312 } 2313 2314 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 2315 { 2316 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2317 2318 PetscFunctionBegin; 2319 *flg = (PetscBool)!ark->imex; 2320 PetscFunctionReturn(PETSC_SUCCESS); 2321 } 2322 2323 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 2324 { 2325 PetscFunctionBegin; 2326 PetscCall(TSReset_ARKIMEX(ts)); 2327 if (ts->dm) { 2328 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 2329 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2330 } 2331 PetscCall(PetscFree(ts->data)); 2332 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL)); 2333 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL)); 2334 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 2335 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 2336 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 2337 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 2338 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", NULL)); 2339 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", NULL)); 2340 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", NULL)); 2341 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_ARKIMEX_FastSlowSplit_C", NULL)); 2342 PetscFunctionReturn(PETSC_SUCCESS); 2343 } 2344 2345 /*MC 2346 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 2347 2348 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2349 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2350 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2351 2352 Options Database Keys: 2353 + -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - Set `TSARKIMEX` scheme type 2354 . -ts_dirk_type <type> - Set `TSDIRK` scheme type 2355 . -ts_arkimex_fully_implicit <true,false> - Solve both parts of the equation implicitly 2356 . -ts_arkimex_fastslowsplit <true,false> - Enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise, 2357 see `TSRHSSplitSetIS()` 2358 - -ts_arkimex_initial_guess_extrapolate - Extrapolate the initial guess for the stage solution from stage values of the previous time step 2359 2360 Level: beginner 2361 2362 Notes: 2363 The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or `-ts_arkimex_type` 2364 2365 If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2366 2367 Methods with an explicit stage can only be used with ODE in which the stiff part $ G(t,X,\dot{X}) $ has the form $ \dot{X} + \hat{G}(t,X)$. 2368 2369 Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2370 2371 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2372 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2373 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 2374 M*/ 2375 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2376 { 2377 TS_ARKIMEX *ark; 2378 PetscBool dirk; 2379 2380 PetscFunctionBegin; 2381 PetscCall(TSARKIMEXInitializePackage()); 2382 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2383 2384 ts->ops->reset = TSReset_ARKIMEX; 2385 ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 2386 ts->ops->destroy = TSDestroy_ARKIMEX; 2387 ts->ops->view = TSView_ARKIMEX; 2388 ts->ops->load = TSLoad_ARKIMEX; 2389 ts->ops->setup = TSSetUp_ARKIMEX; 2390 ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 2391 ts->ops->step = TSStep_ARKIMEX; 2392 ts->ops->interpolate = TSInterpolate_ARKIMEX; 2393 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 2394 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 2395 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 2396 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 2397 ts->ops->getstages = TSGetStages_ARKIMEX; 2398 ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 2399 2400 ts->usessnes = PETSC_TRUE; 2401 2402 PetscCall(PetscNew(&ark)); 2403 ts->data = (void *)ark; 2404 ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 2405 2406 ark->VecsDeltaLam = NULL; 2407 ark->VecsSensiTemp = NULL; 2408 ark->VecsSensiPTemp = NULL; 2409 2410 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2411 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2412 if (!dirk) { 2413 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 2414 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 2415 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX)); 2416 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX)); 2417 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2418 } 2419 PetscFunctionReturn(PETSC_SUCCESS); 2420 } 2421 2422 static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2423 { 2424 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2425 2426 PetscFunctionBegin; 2427 PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2428 PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2429 PetscFunctionReturn(PETSC_SUCCESS); 2430 } 2431 2432 /*@ 2433 TSDIRKSetType - Set the type of `TSDIRK` scheme 2434 2435 Logically Collective 2436 2437 Input Parameters: 2438 + ts - timestepping context 2439 - dirktype - type of `TSDIRK` scheme 2440 2441 Options Database Key: 2442 . -ts_dirkimex_type - set `TSDIRK` scheme type 2443 2444 Level: intermediate 2445 2446 .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2447 @*/ 2448 PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2449 { 2450 PetscFunctionBegin; 2451 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2452 PetscAssertPointer(dirktype, 2); 2453 PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2454 PetscFunctionReturn(PETSC_SUCCESS); 2455 } 2456 2457 /*@ 2458 TSDIRKGetType - Get the type of `TSDIRK` scheme 2459 2460 Logically Collective 2461 2462 Input Parameter: 2463 . ts - timestepping context 2464 2465 Output Parameter: 2466 . dirktype - type of `TSDIRK` scheme 2467 2468 Level: intermediate 2469 2470 .seealso: [](ch_ts), `TSDIRKSetType()` 2471 @*/ 2472 PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2473 { 2474 PetscFunctionBegin; 2475 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2476 PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2477 PetscFunctionReturn(PETSC_SUCCESS); 2478 } 2479 2480 /*MC 2481 TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes. 2482 2483 Level: beginner 2484 2485 Notes: 2486 The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or `-ts_dirk_type`. 2487 The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with: 2488 + E - whether the method has an explicit first stage 2489 . S - whether the method is single diagonal 2490 . P - order of the advancing method 2491 . Q - order of the embedded method 2492 . S - number of stages 2493 . SA - whether the method is stiffly accurate 2494 . L - whether the method is L-stable 2495 - A - whether the method is A-stable 2496 2497 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2498 M*/ 2499 PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2500 { 2501 PetscFunctionBegin; 2502 PetscCall(TSCreate_ARKIMEX(ts)); 2503 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2504 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2505 PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 2506 PetscFunctionReturn(PETSC_SUCCESS); 2507 } 2508 2509 /*@ 2510 TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system 2511 2512 Logically Collective 2513 2514 Input Parameters: 2515 + ts - timestepping context 2516 - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise. 2517 2518 Options Database Key: 2519 . -ts_arkimex_fastslowsplit <true,false> - enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise 2520 2521 Level: intermediate 2522 2523 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()`, `TSRHSSplitSetIS()` 2524 @*/ 2525 PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow) 2526 { 2527 PetscFunctionBegin; 2528 PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow)); 2529 PetscFunctionReturn(PETSC_SUCCESS); 2530 } 2531 2532 /*@ 2533 TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system 2534 2535 Not Collective 2536 2537 Input Parameter: 2538 . ts - timestepping context 2539 2540 Output Parameter: 2541 . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise 2542 2543 Level: intermediate 2544 2545 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()` 2546 @*/ 2547 PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow) 2548 { 2549 PetscFunctionBegin; 2550 PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow)); 2551 PetscFunctionReturn(PETSC_SUCCESS); 2552 } 2553