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) { 1796 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1797 } else *Z = ax->Z; 1798 } 1799 if (Ydot) { 1800 if (dm && dm != ts->dm) { 1801 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1802 } else *Ydot = ax->Ydot; 1803 } 1804 PetscFunctionReturn(PETSC_SUCCESS); 1805 } 1806 1807 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1808 { 1809 PetscFunctionBegin; 1810 if (Z) { 1811 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1812 } 1813 if (Ydot) { 1814 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1815 } 1816 PetscFunctionReturn(PETSC_SUCCESS); 1817 } 1818 1819 /* 1820 DAEs need special handling for algebraic variables when restarting DIRK methods with explicit 1821 first stage. In particular, we need: 1822 - to zero the nonlinear function (in case the dual variables are not consistent in the first step) 1823 - to modify the matrix by calling MatZeroRows with identity on these variables. 1824 */ 1825 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is) 1826 { 1827 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1828 DM dm; 1829 Vec F, W, Xdot; 1830 const PetscScalar *w; 1831 PetscInt nz = 0, n, st; 1832 PetscInt *nzr; 1833 1834 PetscFunctionBegin; 1835 PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */ 1836 PetscCall(DMGetGlobalVector(dm, &Xdot)); 1837 PetscCall(DMGetGlobalVector(dm, &F)); 1838 PetscCall(DMGetGlobalVector(dm, &W)); 1839 PetscCall(VecSet(Xdot, 0.0)); 1840 PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex)); 1841 PetscCall(VecSetRandom(Xdot, NULL)); 1842 PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex)); 1843 PetscCall(VecAXPY(W, -1.0, F)); 1844 PetscCall(VecGetOwnershipRange(W, &st, NULL)); 1845 PetscCall(VecGetLocalSize(W, &n)); 1846 PetscCall(VecGetArrayRead(W, &w)); 1847 for (PetscInt i = 0; i < n; i++) 1848 if (w[i] == 0.0) nz++; 1849 PetscCall(PetscMalloc1(nz, &nzr)); 1850 nz = 0; 1851 for (PetscInt i = 0; i < n; i++) 1852 if (w[i] == 0.0) nzr[nz++] = i + st; 1853 PetscCall(VecRestoreArrayRead(W, &w)); 1854 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is)); 1855 PetscCall(DMRestoreGlobalVector(dm, &Xdot)); 1856 PetscCall(DMRestoreGlobalVector(dm, &F)); 1857 PetscCall(DMRestoreGlobalVector(dm, &W)); 1858 PetscFunctionReturn(PETSC_SUCCESS); 1859 } 1860 1861 /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure 1862 at the finest level, in the DM for coarser solves. */ 1863 static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is) 1864 { 1865 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1866 1867 PetscFunctionBegin; 1868 if (dm && dm != ts->dm) { 1869 PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is)); 1870 } else *alg_is = ax->alg_is; 1871 PetscFunctionReturn(PETSC_SUCCESS); 1872 } 1873 1874 /* This defines the nonlinear equation that is to be solved with SNES */ 1875 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1876 { 1877 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1878 DM dm, dmsave; 1879 Vec Z, Ydot; 1880 IS alg_is; 1881 1882 PetscFunctionBegin; 1883 PetscCall(SNESGetDM(snes, &dm)); 1884 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1885 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 1886 1887 dmsave = ts->dm; 1888 ts->dm = dm; 1889 1890 if (ark->scoeff == PETSC_MAX_REAL) { 1891 /* We are solving F(t_n,x_n,xdot) = 0 to start the method */ 1892 if (!alg_is) { 1893 PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 1894 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is)); 1895 PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is)); 1896 PetscCall(PetscObjectDereference((PetscObject)alg_is)); 1897 PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view")); 1898 } 1899 PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 1900 PetscCall(VecISSet(F, alg_is, 0.0)); 1901 } else { 1902 PetscReal shift = ark->scoeff / ts->time_step; 1903 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1904 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1905 } 1906 1907 ts->dm = dmsave; 1908 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1909 PetscFunctionReturn(PETSC_SUCCESS); 1910 } 1911 1912 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1913 { 1914 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1915 DM dm, dmsave; 1916 Vec Ydot, Z; 1917 PetscReal shift; 1918 IS alg_is; 1919 1920 PetscFunctionBegin; 1921 PetscCall(SNESGetDM(snes, &dm)); 1922 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1923 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1924 /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */ 1925 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 1926 1927 dmsave = ts->dm; 1928 ts->dm = dm; 1929 1930 if (ark->scoeff == PETSC_MAX_REAL) { 1931 PetscBool hasZeroRows; 1932 1933 /* We are solving F(t_n,x_n,xdot) = 0 to start the method 1934 We compute with a very large shift and then scale back the matrix */ 1935 shift = 1.0 / PETSC_MACHINE_EPSILON; 1936 PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1937 PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 1938 PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows)); 1939 if (hasZeroRows) { 1940 PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 1941 /* the default of AIJ is to not keep the pattern! We should probably change it someday */ 1942 PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 1943 PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL)); 1944 } 1945 PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view")); 1946 if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1947 } else { 1948 shift = ark->scoeff / ts->time_step; 1949 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1950 } 1951 ts->dm = dmsave; 1952 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1953 PetscFunctionReturn(PETSC_SUCCESS); 1954 } 1955 1956 static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 1957 { 1958 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1959 1960 PetscFunctionBegin; 1961 if (ns) *ns = ark->tableau->s; 1962 if (Y) *Y = ark->Y; 1963 PetscFunctionReturn(PETSC_SUCCESS); 1964 } 1965 1966 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1967 { 1968 PetscFunctionBegin; 1969 PetscFunctionReturn(PETSC_SUCCESS); 1970 } 1971 1972 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1973 { 1974 TS ts = (TS)ctx; 1975 Vec Z, Z_c; 1976 1977 PetscFunctionBegin; 1978 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 1979 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 1980 PetscCall(MatRestrict(restrct, Z, Z_c)); 1981 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 1982 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 1983 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 1984 PetscFunctionReturn(PETSC_SUCCESS); 1985 } 1986 1987 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1988 { 1989 PetscFunctionBegin; 1990 PetscFunctionReturn(PETSC_SUCCESS); 1991 } 1992 1993 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1994 { 1995 TS ts = (TS)ctx; 1996 Vec Z, Z_c; 1997 1998 PetscFunctionBegin; 1999 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 2000 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 2001 2002 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2003 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2004 2005 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 2006 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 2007 PetscFunctionReturn(PETSC_SUCCESS); 2008 } 2009 2010 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 2011 { 2012 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2013 ARKTableau tab = ark->tableau; 2014 2015 PetscFunctionBegin; 2016 PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 2017 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 2018 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 2019 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 2020 if (ark->extrapolate) { 2021 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 2022 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 2023 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 2024 } 2025 PetscFunctionReturn(PETSC_SUCCESS); 2026 } 2027 2028 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 2029 { 2030 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2031 DM dm; 2032 SNES snes; 2033 2034 PetscFunctionBegin; 2035 if (ark->fastslowsplit) { 2036 PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 2037 } else { 2038 PetscCall(TSARKIMEXTableauSetUp(ts)); 2039 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 2040 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 2041 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 2042 PetscCall(TSGetDM(ts, &dm)); 2043 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 2044 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2045 PetscCall(TSGetSNES(ts, &snes)); 2046 PetscCall(SNESSetDM(snes, dm)); 2047 } 2048 PetscFunctionReturn(PETSC_SUCCESS); 2049 } 2050 2051 static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 2052 { 2053 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2054 ARKTableau tab = ark->tableau; 2055 2056 PetscFunctionBegin; 2057 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 2058 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 2059 if (ts->vecs_sensip) PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); 2060 if (PetscDefined(USE_DEBUG)) { 2061 PetscBool id = PETSC_FALSE; 2062 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 2063 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"); 2064 } 2065 PetscFunctionReturn(PETSC_SUCCESS); 2066 } 2067 2068 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems PetscOptionsObject) 2069 { 2070 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2071 PetscBool dirk; 2072 2073 PetscFunctionBegin; 2074 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2075 PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 2076 { 2077 ARKTableauLink link; 2078 PetscInt count, choice; 2079 PetscBool flg; 2080 const char **namelist; 2081 for (link = ARKTableauList, count = 0; link; link = link->next) { 2082 if (!dirk && link->tab.additive) count++; 2083 if (dirk && !link->tab.additive) count++; 2084 } 2085 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 2086 for (link = ARKTableauList, count = 0; link; link = link->next) { 2087 if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 2088 if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 2089 } 2090 if (dirk) { 2091 PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2092 if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 2093 } else { 2094 PetscBool fastslowsplit; 2095 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2096 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 2097 flg = (PetscBool)!ark->imex; 2098 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 2099 ark->imex = (PetscBool)!flg; 2100 PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg)); 2101 if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit)); 2102 PetscCall(TSARKIMEXGetFastSlowSplit(ts, &fastslowsplit)); 2103 if (fastslowsplit) { 2104 SNES snes; 2105 2106 PetscCall(TSRHSSplitGetSNES(ts, &snes)); 2107 PetscCall(SNESSetFromOptions(snes)); 2108 } 2109 } 2110 PetscCall(PetscFree(namelist)); 2111 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)); 2112 } 2113 PetscOptionsHeadEnd(); 2114 PetscFunctionReturn(PETSC_SUCCESS); 2115 } 2116 2117 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 2118 { 2119 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2120 PetscBool isascii, dirk; 2121 2122 PetscFunctionBegin; 2123 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2124 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 2125 if (isascii) { 2126 PetscViewerFormat format; 2127 ARKTableau tab = ark->tableau; 2128 TSARKIMEXType arktype; 2129 char buf[2048]; 2130 PetscBool flg; 2131 2132 PetscCall(TSARKIMEXGetType(ts, &arktype)); 2133 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 2134 PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 2135 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 2136 PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 2137 PetscCall(PetscViewerGetFormat(viewer, &format)); 2138 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2139 PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 2140 for (PetscInt i = 0; i < tab->s; i++) { 2141 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 2142 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 2143 } 2144 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 2145 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 2146 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 2147 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 2148 } 2149 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 2150 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 2151 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 2152 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 2153 if (!dirk) { 2154 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 2155 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 2156 } 2157 } 2158 PetscFunctionReturn(PETSC_SUCCESS); 2159 } 2160 2161 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 2162 { 2163 SNES snes; 2164 TSAdapt adapt; 2165 2166 PetscFunctionBegin; 2167 PetscCall(TSGetAdapt(ts, &adapt)); 2168 PetscCall(TSAdaptLoad(adapt, viewer)); 2169 PetscCall(TSGetSNES(ts, &snes)); 2170 PetscCall(SNESLoad(snes, viewer)); 2171 /* function and Jacobian context for SNES when used with TS is always ts object */ 2172 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 2173 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 2174 PetscFunctionReturn(PETSC_SUCCESS); 2175 } 2176 2177 /*@ 2178 TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 2179 2180 Logically Collective 2181 2182 Input Parameters: 2183 + ts - timestepping context 2184 - arktype - type of `TSARKIMEX` scheme 2185 2186 Options Database Key: 2187 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 2188 2189 Level: intermediate 2190 2191 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 2192 `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 2193 @*/ 2194 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 2195 { 2196 PetscFunctionBegin; 2197 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2198 PetscAssertPointer(arktype, 2); 2199 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 2200 PetscFunctionReturn(PETSC_SUCCESS); 2201 } 2202 2203 /*@ 2204 TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 2205 2206 Logically Collective 2207 2208 Input Parameter: 2209 . ts - timestepping context 2210 2211 Output Parameter: 2212 . arktype - type of `TSARKIMEX` scheme 2213 2214 Level: intermediate 2215 2216 .seealso: [](ch_ts), `TSARKIMEX` 2217 @*/ 2218 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2219 { 2220 PetscFunctionBegin; 2221 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2222 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 2223 PetscFunctionReturn(PETSC_SUCCESS); 2224 } 2225 2226 /*@ 2227 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 2228 2229 Logically Collective 2230 2231 Input Parameters: 2232 + ts - timestepping context 2233 - flg - `PETSC_TRUE` for fully implicit 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 /* ------------------------------------------------------------ */ 2346 /*MC 2347 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 2348 2349 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2350 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2351 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2352 2353 Level: beginner 2354 2355 Notes: 2356 The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 2357 2358 If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2359 2360 Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 2361 2362 Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2363 2364 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2365 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2366 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 2367 M*/ 2368 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2369 { 2370 TS_ARKIMEX *ark; 2371 PetscBool dirk; 2372 2373 PetscFunctionBegin; 2374 PetscCall(TSARKIMEXInitializePackage()); 2375 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2376 2377 ts->ops->reset = TSReset_ARKIMEX; 2378 ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 2379 ts->ops->destroy = TSDestroy_ARKIMEX; 2380 ts->ops->view = TSView_ARKIMEX; 2381 ts->ops->load = TSLoad_ARKIMEX; 2382 ts->ops->setup = TSSetUp_ARKIMEX; 2383 ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 2384 ts->ops->step = TSStep_ARKIMEX; 2385 ts->ops->interpolate = TSInterpolate_ARKIMEX; 2386 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 2387 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 2388 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 2389 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 2390 ts->ops->getstages = TSGetStages_ARKIMEX; 2391 ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 2392 2393 ts->usessnes = PETSC_TRUE; 2394 2395 PetscCall(PetscNew(&ark)); 2396 ts->data = (void *)ark; 2397 ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 2398 2399 ark->VecsDeltaLam = NULL; 2400 ark->VecsSensiTemp = NULL; 2401 ark->VecsSensiPTemp = NULL; 2402 2403 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2404 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2405 if (!dirk) { 2406 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 2407 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 2408 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX)); 2409 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX)); 2410 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2411 } 2412 PetscFunctionReturn(PETSC_SUCCESS); 2413 } 2414 2415 /* ------------------------------------------------------------ */ 2416 2417 static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2418 { 2419 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2420 2421 PetscFunctionBegin; 2422 PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2423 PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2424 PetscFunctionReturn(PETSC_SUCCESS); 2425 } 2426 2427 /*@ 2428 TSDIRKSetType - Set the type of `TSDIRK` scheme 2429 2430 Logically Collective 2431 2432 Input Parameters: 2433 + ts - timestepping context 2434 - dirktype - type of `TSDIRK` scheme 2435 2436 Options Database Key: 2437 . -ts_dirkimex_type - set `TSDIRK` scheme type 2438 2439 Level: intermediate 2440 2441 .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2442 @*/ 2443 PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2444 { 2445 PetscFunctionBegin; 2446 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2447 PetscAssertPointer(dirktype, 2); 2448 PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2449 PetscFunctionReturn(PETSC_SUCCESS); 2450 } 2451 2452 /*@ 2453 TSDIRKGetType - Get the type of `TSDIRK` scheme 2454 2455 Logically Collective 2456 2457 Input Parameter: 2458 . ts - timestepping context 2459 2460 Output Parameter: 2461 . dirktype - type of `TSDIRK` scheme 2462 2463 Level: intermediate 2464 2465 .seealso: [](ch_ts), `TSDIRKSetType()` 2466 @*/ 2467 PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2468 { 2469 PetscFunctionBegin; 2470 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2471 PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2472 PetscFunctionReturn(PETSC_SUCCESS); 2473 } 2474 2475 /*MC 2476 TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes. 2477 2478 Level: beginner 2479 2480 Notes: 2481 The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type. 2482 The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with: 2483 + E - whether the method has an explicit first stage 2484 . S - whether the method is single diagonal 2485 . P - order of the advancing method 2486 . Q - order of the embedded method 2487 . S - number of stages 2488 . SA - whether the method is stiffly accurate 2489 . L - whether the method is L-stable 2490 - A - whether the method is A-stable 2491 2492 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2493 M*/ 2494 PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2495 { 2496 PetscFunctionBegin; 2497 PetscCall(TSCreate_ARKIMEX(ts)); 2498 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2499 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2500 PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 2501 PetscFunctionReturn(PETSC_SUCCESS); 2502 } 2503 2504 /*@ 2505 TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system 2506 2507 Logically Collective 2508 2509 Input Parameters: 2510 + ts - timestepping context 2511 - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise. 2512 2513 Options Database Key: 2514 . -ts_arkimex_fastslowsplit - <true,false> 2515 2516 Level: intermediate 2517 2518 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()` 2519 @*/ 2520 PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow) 2521 { 2522 PetscFunctionBegin; 2523 PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow)); 2524 PetscFunctionReturn(PETSC_SUCCESS); 2525 } 2526 2527 /*@ 2528 TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system 2529 2530 Not Collective 2531 2532 Input Parameter: 2533 . ts - timestepping context 2534 2535 Output Parameter: 2536 . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise 2537 2538 Level: intermediate 2539 2540 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()` 2541 @*/ 2542 PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow) 2543 { 2544 PetscFunctionBegin; 2545 PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow)); 2546 PetscFunctionReturn(PETSC_SUCCESS); 2547 } 2548