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 if (done) *done = PETSC_FALSE; 1295 else 1296 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "ARKIMEX '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.", tab->name, 1297 tab->order, order); 1298 PetscFunctionReturn(PETSC_SUCCESS); 1299 } 1300 1301 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 1302 { 1303 Vec Udot, Y1, Y2; 1304 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1305 PetscReal norm; 1306 1307 PetscFunctionBegin; 1308 PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 1309 PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 1310 PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 1311 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 1312 PetscCall(VecSetRandom(Udot, NULL)); 1313 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 1314 PetscCall(VecAXPY(Y2, -1.0, Y1)); 1315 PetscCall(VecAXPY(Y2, -1.0, Udot)); 1316 PetscCall(VecNorm(Y2, NORM_2, &norm)); 1317 if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 1318 *id = PETSC_TRUE; 1319 } else { 1320 *id = PETSC_FALSE; 1321 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)); 1322 } 1323 PetscCall(VecDestroy(&Udot)); 1324 PetscCall(VecDestroy(&Y1)); 1325 PetscCall(VecDestroy(&Y2)); 1326 PetscFunctionReturn(PETSC_SUCCESS); 1327 } 1328 1329 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS, PetscReal, Vec, IS *); 1330 1331 static PetscErrorCode TSStep_ARKIMEX(TS ts) 1332 { 1333 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1334 ARKTableau tab = ark->tableau; 1335 const PetscInt s = tab->s; 1336 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 1337 PetscScalar *w = ark->work; 1338 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 1339 PetscBool extrapolate = ark->extrapolate; 1340 TSAdapt adapt; 1341 SNES snes; 1342 PetscInt i, j, its, lits; 1343 PetscInt rejections = 0; 1344 PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE; 1345 PetscReal next_time_step = ts->time_step; 1346 1347 PetscFunctionBegin; 1348 if (ark->extrapolate && !ark->Y_prev) { 1349 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1350 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1351 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1352 } 1353 1354 if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE)); 1355 if (!hasE) dirk = PETSC_TRUE; 1356 1357 if (!ts->steprollback) { 1358 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 1359 PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 1360 } 1361 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 1362 for (i = 0; i < s; i++) { 1363 PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 1364 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 1365 if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 1366 } 1367 } 1368 } 1369 1370 /* 1371 For fully implicit formulations we must solve the equations 1372 1373 F(t_n,x_n,xdot) = 0 1374 1375 for the explicit first stage. 1376 Here we call SNESSolve using PETSC_MAX_REAL as shift to flag it. 1377 Special handling is inside SNESTSFormFunction_ARKIMEX and SNESTSFormJacobian_ARKIMEX 1378 We compute Ydot0 if we restart the step or if we resized the problem after remeshing 1379 */ 1380 if (dirk && tab->explicit_first_stage && (ts->steprestart || ts->stepresize)) { 1381 ark->scoeff = PETSC_MAX_REAL; 1382 PetscCall(VecCopy(ts->vec_sol, Z)); 1383 if (!ark->alg_is) { 1384 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ts->ptime, Z, &ark->alg_is)); 1385 PetscCall(ISViewFromOptions(ark->alg_is, (PetscObject)ts, "-ts_arkimex_algebraic_is_view")); 1386 } 1387 PetscCall(TSGetSNES(ts, &snes)); 1388 PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, 1)); 1389 PetscCall(SNESSolve(snes, NULL, Ydot0)); 1390 if (ark->alg_is) PetscCall(VecISSet(Ydot0, ark->alg_is, 0.0)); 1391 PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes, (PetscObject)snes, -1)); 1392 } 1393 1394 /* For IMEX we compute a step */ 1395 if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 1396 TS ts_start; 1397 if (PetscDefined(USE_DEBUG) && hasE) { 1398 PetscBool id = PETSC_FALSE; 1399 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1400 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"); 1401 } 1402 PetscCall(TSClone(ts, &ts_start)); 1403 PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 1404 PetscCall(TSSetTime(ts_start, ts->ptime)); 1405 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 1406 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 1407 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 1408 PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 1409 PetscCall(TSSetType(ts_start, TSARKIMEX)); 1410 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 1411 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 1412 1413 PetscCall(TSRestartStep(ts_start)); 1414 PetscCall(TSSolve(ts_start, ts->vec_sol)); 1415 PetscCall(TSGetTime(ts_start, &ts->ptime)); 1416 PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 1417 1418 { /* Save the initial slope for the next step */ 1419 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 1420 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 1421 } 1422 ts->steps++; 1423 1424 /* Set the correct TS in SNES */ 1425 /* We'll try to bypass this by changing the method on the fly */ 1426 { 1427 PetscCall(TSGetSNES(ts, &snes)); 1428 PetscCall(TSSetSNES(ts, snes)); 1429 } 1430 PetscCall(TSDestroy(&ts_start)); 1431 } 1432 1433 ark->status = TS_STEP_INCOMPLETE; 1434 while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 1435 PetscReal t = ts->ptime; 1436 PetscReal h = ts->time_step; 1437 for (i = 0; i < s; i++) { 1438 ark->stage_time = t + h * ct[i]; 1439 PetscCall(TSPreStage(ts, ark->stage_time)); 1440 if (At[i * s + i] == 0) { /* This stage is explicit */ 1441 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"); 1442 PetscCall(VecCopy(ts->vec_sol, Y[i])); 1443 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1444 PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 1445 if (tab->additive && hasE) { 1446 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1447 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 1448 } 1449 PetscCall(TSGetSNES(ts, &snes)); 1450 PetscCall(SNESResetCounters(snes)); 1451 } else { 1452 ark->scoeff = 1. / At[i * s + i]; 1453 /* Ydot = shift*(Y-Z) */ 1454 PetscCall(VecCopy(ts->vec_sol, Z)); 1455 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1456 PetscCall(VecMAXPY(Z, i, w, YdotI)); 1457 if (tab->additive && hasE) { 1458 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1459 PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 1460 } 1461 if (extrapolate && !ts->steprestart) { 1462 /* Initial guess extrapolated from previous time step stage values */ 1463 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 1464 } else { 1465 /* Initial guess taken from last stage */ 1466 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 1467 } 1468 PetscCall(TSGetSNES(ts, &snes)); 1469 PetscCall(SNESSolve(snes, NULL, Y[i])); 1470 PetscCall(SNESGetIterationNumber(snes, &its)); 1471 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 1472 ts->snes_its += its; 1473 ts->ksp_its += lits; 1474 PetscCall(TSGetAdapt(ts, &adapt)); 1475 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 1476 if (!stageok) { 1477 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 1478 * use extrapolation to initialize the solves on the next attempt. */ 1479 extrapolate = PETSC_FALSE; 1480 goto reject_step; 1481 } 1482 } 1483 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { 1484 if (i == 0 && tab->explicit_first_stage) { 1485 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", 1486 ((PetscObject)ts)->type_name, ark->tableau->name); 1487 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 1488 } else { 1489 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1490 } 1491 } else { 1492 if (i == 0 && tab->explicit_first_stage) { 1493 PetscCall(VecZeroEntries(Ydot)); 1494 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 1495 PetscCall(VecScale(YdotI[i], -1.0)); 1496 } else { 1497 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1498 } 1499 if (hasE) { 1500 if (ark->imex) { 1501 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 1502 } else { 1503 PetscCall(VecZeroEntries(YdotRHS[i])); 1504 } 1505 } 1506 } 1507 PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 1508 } 1509 1510 ark->status = TS_STEP_INCOMPLETE; 1511 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 1512 ark->status = TS_STEP_PENDING; 1513 PetscCall(TSGetAdapt(ts, &adapt)); 1514 PetscCall(TSAdaptCandidatesClear(adapt)); 1515 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 1516 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1517 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1518 if (!accept) { /* Roll back the current step */ 1519 PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1520 ts->time_step = next_time_step; 1521 goto reject_step; 1522 } 1523 1524 ts->ptime += ts->time_step; 1525 ts->time_step = next_time_step; 1526 break; 1527 1528 reject_step: 1529 ts->reject++; 1530 accept = PETSC_FALSE; 1531 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1532 ts->reason = TS_DIVERGED_STEP_REJECTED; 1533 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1534 } 1535 } 1536 PetscFunctionReturn(PETSC_SUCCESS); 1537 } 1538 1539 /* 1540 This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 1541 Udot = H(t,U) + G(t,U) 1542 This corresponds to F(t,U,Udot) = Udot-H(t,U). 1543 1544 The complete adjoint equations are 1545 (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 1546 dGdU (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 1547 + dHdU (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 1548 lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 1549 mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 1550 + dGdP (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 1551 + dHdP (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 1552 mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 1553 where shift = 1/(at[i][i]*h) 1554 1555 If at[i][i] is 0, the first equation falls back to 1556 lambda_s[i] = h ( 1557 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 1558 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 1559 1560 */ 1561 static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 1562 { 1563 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1564 ARKTableau tab = ark->tableau; 1565 const PetscInt s = tab->s; 1566 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 1567 PetscScalar *w = ark->work; 1568 Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 1569 Mat Jex, Jim, Jimpre; 1570 PetscInt i, j, nadj; 1571 PetscReal t = ts->ptime, stage_time_ex; 1572 PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 1573 KSP ksp; 1574 1575 PetscFunctionBegin; 1576 ark->status = TS_STEP_INCOMPLETE; 1577 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1578 PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 1579 PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 1580 1581 for (i = s - 1; i >= 0; i--) { 1582 ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 1583 stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 1584 if (At[i * s + i] == 0) { // This stage is explicit 1585 ark->scoeff = 0.; 1586 } else { 1587 ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 1588 } 1589 PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 1590 PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 1591 if (ts->vecs_sensip) { 1592 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 1593 PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 1594 } 1595 /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 1596 for (nadj = 0; nadj < ts->numcost; nadj++) { 1597 /* build implicit part */ 1598 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1599 if (s - i - 1 > 0) { 1600 /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 1601 for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 1602 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1603 } 1604 /* Temp = Temp - bt[i] lambda_{n+1} */ 1605 PetscCall(VecAXPY(VecsSensiTemp[nadj], -bt[i], ts->vecs_sensi[nadj])); 1606 if (bt[i] || s - i - 1 > 0) { 1607 /* (shift I - dHdU) Temp */ 1608 PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 1609 /* cancel out shift Temp where shift=-scoeff/h */ 1610 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 1611 if (ts->vecs_sensip) { 1612 /* - dHdP Temp */ 1613 PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1614 /* mu_n += -h dHdP Temp */ 1615 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1616 } 1617 } else { 1618 PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 1619 } 1620 /* build explicit part */ 1621 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1622 if (s - i - 1 > 0) { 1623 /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 1624 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 1625 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1626 } 1627 /* Temp = Temp + b[i] lambda_{n+1} */ 1628 PetscCall(VecAXPY(VecsSensiTemp[nadj], b[i], ts->vecs_sensi[nadj])); 1629 if (b[i] || s - i - 1 > 0) { 1630 /* dGdU Temp */ 1631 PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1632 if (ts->vecs_sensip) { 1633 /* dGdP Temp */ 1634 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1635 /* mu_n += h dGdP Temp */ 1636 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1637 } 1638 } 1639 /* Build LHS for first-order adjoint */ 1640 if (At[i * s + i] == 0) { // This stage is explicit 1641 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 1642 } else { 1643 KSP ksp; 1644 KSPConvergedReason kspreason; 1645 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1646 PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 1647 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 1648 PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1649 PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 1650 if (kspreason < 0) { 1651 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 1652 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj)); 1653 } 1654 if (ts->vecs_sensip) { 1655 /* -dHdP lambda_s[i] */ 1656 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 1657 /* mu_n += h at[i][i] dHdP lambda_s[i] */ 1658 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 1659 } 1660 } 1661 } 1662 } 1663 for (j = 0; j < s; j++) w[j] = 1.0; 1664 for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 1665 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 1666 ark->status = TS_STEP_COMPLETE; 1667 PetscFunctionReturn(PETSC_SUCCESS); 1668 } 1669 1670 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1671 { 1672 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1673 ARKTableau tab = ark->tableau; 1674 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1675 PetscReal h; 1676 PetscReal tt, t; 1677 PetscScalar *bt = ark->work, *b = ark->work + s; 1678 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1679 1680 PetscFunctionBegin; 1681 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name); 1682 switch (ark->status) { 1683 case TS_STEP_INCOMPLETE: 1684 case TS_STEP_PENDING: 1685 h = ts->time_step; 1686 t = (itime - ts->ptime) / h; 1687 break; 1688 case TS_STEP_COMPLETE: 1689 h = ts->ptime - ts->ptime_prev; 1690 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1691 break; 1692 default: 1693 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1694 } 1695 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1696 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1697 for (i = 0; i < s; i++) { 1698 bt[i] += h * Bt[i * pinterp + j] * tt; 1699 b[i] += h * B[i * pinterp + j] * tt; 1700 } 1701 } 1702 PetscCall(VecCopy(ark->Y[0], X)); 1703 PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1704 if (tab->additive) { 1705 PetscBool hasE; 1706 PetscCall(TSHasRHSFunction(ts, &hasE)); 1707 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1708 } 1709 PetscFunctionReturn(PETSC_SUCCESS); 1710 } 1711 1712 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1713 { 1714 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1715 ARKTableau tab = ark->tableau; 1716 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1717 PetscReal h, h_prev, t, tt; 1718 PetscScalar *bt = ark->work, *b = ark->work + s; 1719 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1720 1721 PetscFunctionBegin; 1722 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 1723 h = ts->time_step; 1724 h_prev = ts->ptime - ts->ptime_prev; 1725 t = 1 + h / h_prev * c; 1726 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1727 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1728 for (i = 0; i < s; i++) { 1729 bt[i] += h * Bt[i * pinterp + j] * tt; 1730 b[i] += h * B[i * pinterp + j] * tt; 1731 } 1732 } 1733 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 1734 PetscCall(VecCopy(ark->Y_prev[0], X)); 1735 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1736 if (tab->additive) { 1737 PetscBool hasE; 1738 PetscCall(TSHasRHSFunction(ts, &hasE)); 1739 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1740 } 1741 PetscFunctionReturn(PETSC_SUCCESS); 1742 } 1743 1744 static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1745 { 1746 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1747 ARKTableau tab = ark->tableau; 1748 1749 PetscFunctionBegin; 1750 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1751 PetscCall(PetscFree(ark->work)); 1752 PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 1753 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 1754 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 1755 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 1756 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 1757 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 1758 PetscFunctionReturn(PETSC_SUCCESS); 1759 } 1760 1761 static PetscErrorCode TSReset_ARKIMEX(TS ts) 1762 { 1763 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1764 1765 PetscFunctionBegin; 1766 if (ark->fastslowsplit) { 1767 PetscTryMethod(ts, "TSReset_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 1768 } else { 1769 PetscCall(TSARKIMEXTableauReset(ts)); 1770 PetscCall(VecDestroy(&ark->Ydot)); 1771 PetscCall(VecDestroy(&ark->Ydot0)); 1772 PetscCall(VecDestroy(&ark->Z)); 1773 PetscCall(ISDestroy(&ark->alg_is)); 1774 } 1775 PetscFunctionReturn(PETSC_SUCCESS); 1776 } 1777 1778 static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 1779 { 1780 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1781 ARKTableau tab = ark->tableau; 1782 1783 PetscFunctionBegin; 1784 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 1785 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 1786 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 1787 PetscFunctionReturn(PETSC_SUCCESS); 1788 } 1789 1790 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1791 { 1792 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1793 1794 PetscFunctionBegin; 1795 if (Z) { 1796 if (dm && dm != ts->dm) { 1797 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1798 } else *Z = ax->Z; 1799 } 1800 if (Ydot) { 1801 if (dm && dm != ts->dm) { 1802 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1803 } else *Ydot = ax->Ydot; 1804 } 1805 PetscFunctionReturn(PETSC_SUCCESS); 1806 } 1807 1808 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1809 { 1810 PetscFunctionBegin; 1811 if (Z) { 1812 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1813 } 1814 if (Ydot) { 1815 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1816 } 1817 PetscFunctionReturn(PETSC_SUCCESS); 1818 } 1819 1820 /* 1821 DAEs need special handling for algebraic variables when restarting DIRK methods with explicit 1822 first stage. In particular, we need: 1823 - to zero the nonlinear function (in case the dual variables are not consistent in the first step) 1824 - to modify the preconditioning matrix by calling MatZeroRows with identity on these variables. 1825 */ 1826 static PetscErrorCode TSARKIMEXComputeAlgebraicIS(TS ts, PetscReal time, Vec X, IS *alg_is) 1827 { 1828 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1829 DM dm; 1830 Vec F, W, Xdot; 1831 const PetscScalar *w; 1832 PetscInt nz = 0, n, st; 1833 PetscInt *nzr; 1834 1835 PetscFunctionBegin; 1836 PetscCall(TSGetDM(ts, &dm)); /* may be already from SNES */ 1837 PetscCall(DMGetGlobalVector(dm, &Xdot)); 1838 PetscCall(DMGetGlobalVector(dm, &F)); 1839 PetscCall(DMGetGlobalVector(dm, &W)); 1840 PetscCall(VecSet(Xdot, 0.0)); 1841 PetscCall(TSComputeIFunction(ts, time, X, Xdot, F, ark->imex)); 1842 PetscCall(VecSetRandom(Xdot, NULL)); 1843 PetscCall(TSComputeIFunction(ts, time, X, Xdot, W, ark->imex)); 1844 PetscCall(VecAXPY(W, -1.0, F)); 1845 PetscCall(VecGetOwnershipRange(W, &st, NULL)); 1846 PetscCall(VecGetLocalSize(W, &n)); 1847 PetscCall(VecGetArrayRead(W, &w)); 1848 for (PetscInt i = 0; i < n; i++) 1849 if (w[i] == 0.0) nz++; 1850 PetscCall(PetscMalloc1(nz, &nzr)); 1851 nz = 0; 1852 for (PetscInt i = 0; i < n; i++) 1853 if (w[i] == 0.0) nzr[nz++] = i + st; 1854 PetscCall(VecRestoreArrayRead(W, &w)); 1855 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nz, nzr, PETSC_OWN_POINTER, alg_is)); 1856 PetscCall(DMRestoreGlobalVector(dm, &Xdot)); 1857 PetscCall(DMRestoreGlobalVector(dm, &F)); 1858 PetscCall(DMRestoreGlobalVector(dm, &W)); 1859 PetscFunctionReturn(PETSC_SUCCESS); 1860 } 1861 1862 /* As for the method specific Z and Ydot, we store the algebraic IS in the ARKIMEX data structure 1863 at the finest level, in the DM for coarser solves. */ 1864 static PetscErrorCode TSARKIMEXGetAlgebraicIS(TS ts, DM dm, IS *alg_is) 1865 { 1866 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1867 1868 PetscFunctionBegin; 1869 if (dm && dm != ts->dm) { 1870 PetscCall(PetscObjectQuery((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject *)alg_is)); 1871 } else *alg_is = ax->alg_is; 1872 PetscFunctionReturn(PETSC_SUCCESS); 1873 } 1874 1875 /* This defines the nonlinear equation that is to be solved with SNES */ 1876 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1877 { 1878 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1879 DM dm, dmsave; 1880 Vec Z, Ydot; 1881 IS alg_is; 1882 1883 PetscFunctionBegin; 1884 PetscCall(SNESGetDM(snes, &dm)); 1885 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1886 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 1887 1888 dmsave = ts->dm; 1889 ts->dm = dm; 1890 1891 if (ark->scoeff == PETSC_MAX_REAL) { 1892 /* We are solving F(t_n,x_n,xdot) = 0 to start the method */ 1893 if (!alg_is) { 1894 PetscCheck(dmsave != ts->dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 1895 PetscCall(TSARKIMEXComputeAlgebraicIS(ts, ark->stage_time, Z, &alg_is)); 1896 PetscCall(PetscObjectCompose((PetscObject)dm, "TSARKIMEX_ALG_IS", (PetscObject)alg_is)); 1897 PetscCall(PetscObjectDereference((PetscObject)alg_is)); 1898 PetscCall(ISViewFromOptions(alg_is, (PetscObject)snes, "-ts_arkimex_algebraic_is_view")); 1899 } 1900 PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 1901 PetscCall(VecISSet(F, alg_is, 0.0)); 1902 } else { 1903 PetscReal shift = ark->scoeff / ts->time_step; 1904 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1905 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1906 } 1907 1908 ts->dm = dmsave; 1909 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1910 PetscFunctionReturn(PETSC_SUCCESS); 1911 } 1912 1913 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1914 { 1915 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1916 DM dm, dmsave; 1917 Vec Ydot, Z; 1918 PetscReal shift; 1919 IS alg_is; 1920 1921 PetscFunctionBegin; 1922 PetscCall(SNESGetDM(snes, &dm)); 1923 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1924 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1925 /* alg_is has been computed in SNESTSFormFunction_ARKIMEX */ 1926 if (ark->scoeff == PETSC_MAX_REAL) PetscCall(TSARKIMEXGetAlgebraicIS(ts, dm, &alg_is)); 1927 1928 dmsave = ts->dm; 1929 ts->dm = dm; 1930 1931 if (ark->scoeff == PETSC_MAX_REAL) { 1932 PetscBool hasZeroRows; 1933 1934 /* We are solving F(t_n,x_n,xdot) = 0 to start the method 1935 We compute with a very large shift and then scale back the matrix */ 1936 shift = 1.0 / PETSC_MACHINE_EPSILON; 1937 PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1938 PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 1939 PetscCall(MatHasOperation(B, MATOP_ZERO_ROWS, &hasZeroRows)); 1940 if (hasZeroRows) { 1941 PetscCheck(alg_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing algebraic IS"); 1942 /* the default of AIJ is to not keep the pattern! We should probably change it someday */ 1943 PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 1944 PetscCall(MatZeroRowsIS(B, alg_is, 1.0, NULL, NULL)); 1945 } 1946 PetscCall(MatViewFromOptions(B, (PetscObject)snes, "-ts_arkimex_alg_mat_view")); 1947 if (A != B) PetscCall(MatScale(A, PETSC_MACHINE_EPSILON)); 1948 } else { 1949 shift = ark->scoeff / ts->time_step; 1950 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1951 } 1952 ts->dm = dmsave; 1953 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1954 PetscFunctionReturn(PETSC_SUCCESS); 1955 } 1956 1957 static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 1958 { 1959 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1960 1961 PetscFunctionBegin; 1962 if (ns) *ns = ark->tableau->s; 1963 if (Y) *Y = ark->Y; 1964 PetscFunctionReturn(PETSC_SUCCESS); 1965 } 1966 1967 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1968 { 1969 PetscFunctionBegin; 1970 PetscFunctionReturn(PETSC_SUCCESS); 1971 } 1972 1973 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1974 { 1975 TS ts = (TS)ctx; 1976 Vec Z, Z_c; 1977 1978 PetscFunctionBegin; 1979 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 1980 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 1981 PetscCall(MatRestrict(restrct, Z, Z_c)); 1982 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 1983 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 1984 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 1985 PetscFunctionReturn(PETSC_SUCCESS); 1986 } 1987 1988 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1989 { 1990 PetscFunctionBegin; 1991 PetscFunctionReturn(PETSC_SUCCESS); 1992 } 1993 1994 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1995 { 1996 TS ts = (TS)ctx; 1997 Vec Z, Z_c; 1998 1999 PetscFunctionBegin; 2000 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 2001 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 2002 2003 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2004 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 2005 2006 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 2007 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 2008 PetscFunctionReturn(PETSC_SUCCESS); 2009 } 2010 2011 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 2012 { 2013 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2014 ARKTableau tab = ark->tableau; 2015 2016 PetscFunctionBegin; 2017 PetscCall(PetscMalloc1(2 * tab->s, &ark->work)); 2018 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 2019 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 2020 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 2021 if (ark->extrapolate) { 2022 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 2023 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 2024 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 2025 } 2026 PetscFunctionReturn(PETSC_SUCCESS); 2027 } 2028 2029 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 2030 { 2031 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2032 DM dm; 2033 SNES snes; 2034 2035 PetscFunctionBegin; 2036 if (ark->fastslowsplit) { 2037 PetscTryMethod(ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", (TS), (ts)); 2038 } else { 2039 PetscCall(TSARKIMEXTableauSetUp(ts)); 2040 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 2041 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 2042 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 2043 PetscCall(TSGetDM(ts, &dm)); 2044 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 2045 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2046 PetscCall(TSGetSNES(ts, &snes)); 2047 PetscCall(SNESSetDM(snes, dm)); 2048 } 2049 PetscFunctionReturn(PETSC_SUCCESS); 2050 } 2051 2052 static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 2053 { 2054 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2055 ARKTableau tab = ark->tableau; 2056 2057 PetscFunctionBegin; 2058 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 2059 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 2060 if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); } 2061 if (PetscDefined(USE_DEBUG)) { 2062 PetscBool id = PETSC_FALSE; 2063 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 2064 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"); 2065 } 2066 PetscFunctionReturn(PETSC_SUCCESS); 2067 } 2068 2069 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems PetscOptionsObject) 2070 { 2071 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2072 PetscBool dirk; 2073 2074 PetscFunctionBegin; 2075 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2076 PetscOptionsHeadBegin(PetscOptionsObject, dirk ? "DIRK ODE solver options" : "ARKIMEX ODE solver options"); 2077 { 2078 ARKTableauLink link; 2079 PetscInt count, choice; 2080 PetscBool flg; 2081 const char **namelist; 2082 for (link = ARKTableauList, count = 0; link; link = link->next) { 2083 if (!dirk && link->tab.additive) count++; 2084 if (dirk && !link->tab.additive) count++; 2085 } 2086 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 2087 for (link = ARKTableauList, count = 0; link; link = link->next) { 2088 if (!dirk && link->tab.additive) namelist[count++] = link->tab.name; 2089 if (dirk && !link->tab.additive) namelist[count++] = link->tab.name; 2090 } 2091 if (dirk) { 2092 PetscCall(PetscOptionsEList("-ts_dirk_type", "Family of DIRK method", "TSDIRKSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2093 if (flg) PetscCall(TSDIRKSetType(ts, namelist[choice])); 2094 } else { 2095 PetscBool fastslowsplit; 2096 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 2097 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 2098 flg = (PetscBool)!ark->imex; 2099 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 2100 ark->imex = (PetscBool)!flg; 2101 PetscCall(PetscOptionsBool("-ts_arkimex_fastslowsplit", "Use ARK IMEX for fast-slow systems", "TSARKIMEXSetFastSlowSplit", ark->fastslowsplit, &fastslowsplit, &flg)); 2102 if (flg) PetscCall(TSARKIMEXSetFastSlowSplit(ts, fastslowsplit)); 2103 PetscCall(TSARKIMEXGetFastSlowSplit(ts, &fastslowsplit)); 2104 if (fastslowsplit) { 2105 SNES snes; 2106 2107 PetscCall(TSRHSSplitGetSNES(ts, &snes)); 2108 PetscCall(SNESSetFromOptions(snes)); 2109 } 2110 } 2111 PetscCall(PetscFree(namelist)); 2112 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)); 2113 } 2114 PetscOptionsHeadEnd(); 2115 PetscFunctionReturn(PETSC_SUCCESS); 2116 } 2117 2118 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 2119 { 2120 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2121 PetscBool iascii, dirk; 2122 2123 PetscFunctionBegin; 2124 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2125 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 2126 if (iascii) { 2127 PetscViewerFormat format; 2128 ARKTableau tab = ark->tableau; 2129 TSARKIMEXType arktype; 2130 char buf[2048]; 2131 PetscBool flg; 2132 2133 PetscCall(TSARKIMEXGetType(ts, &arktype)); 2134 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 2135 PetscCall(PetscViewerASCIIPrintf(viewer, " %s %s\n", dirk ? "DIRK" : "ARK IMEX", arktype)); 2136 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 2137 PetscCall(PetscViewerASCIIPrintf(viewer, " %sabscissa ct = %s\n", dirk ? "" : "Stiff ", buf)); 2138 PetscCall(PetscViewerGetFormat(viewer, &format)); 2139 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 2140 PetscCall(PetscViewerASCIIPrintf(viewer, " %sAt =\n", dirk ? "" : "Stiff ")); 2141 for (PetscInt i = 0; i < tab->s; i++) { 2142 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->At + i * tab->s)); 2143 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", buf)); 2144 } 2145 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bt)); 2146 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbt = %s\n", dirk ? "" : "Stiff ", buf)); 2147 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->bembedt)); 2148 PetscCall(PetscViewerASCIIPrintf(viewer, " %sbet = %s\n", dirk ? "" : "Stiff ", buf)); 2149 } 2150 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 2151 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 2152 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 2153 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 2154 if (!dirk) { 2155 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 2156 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 2157 } 2158 } 2159 PetscFunctionReturn(PETSC_SUCCESS); 2160 } 2161 2162 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 2163 { 2164 SNES snes; 2165 TSAdapt adapt; 2166 2167 PetscFunctionBegin; 2168 PetscCall(TSGetAdapt(ts, &adapt)); 2169 PetscCall(TSAdaptLoad(adapt, viewer)); 2170 PetscCall(TSGetSNES(ts, &snes)); 2171 PetscCall(SNESLoad(snes, viewer)); 2172 /* function and Jacobian context for SNES when used with TS is always ts object */ 2173 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 2174 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 2175 PetscFunctionReturn(PETSC_SUCCESS); 2176 } 2177 2178 /*@ 2179 TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 2180 2181 Logically Collective 2182 2183 Input Parameters: 2184 + ts - timestepping context 2185 - arktype - type of `TSARKIMEX` scheme 2186 2187 Options Database Key: 2188 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 2189 2190 Level: intermediate 2191 2192 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 2193 `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 2194 @*/ 2195 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 2196 { 2197 PetscFunctionBegin; 2198 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2199 PetscAssertPointer(arktype, 2); 2200 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 2201 PetscFunctionReturn(PETSC_SUCCESS); 2202 } 2203 2204 /*@ 2205 TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 2206 2207 Logically Collective 2208 2209 Input Parameter: 2210 . ts - timestepping context 2211 2212 Output Parameter: 2213 . arktype - type of `TSARKIMEX` scheme 2214 2215 Level: intermediate 2216 2217 .seealso: [](ch_ts), `TSARKIMEX` 2218 @*/ 2219 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 2220 { 2221 PetscFunctionBegin; 2222 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2223 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 2224 PetscFunctionReturn(PETSC_SUCCESS); 2225 } 2226 2227 /*@ 2228 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 2229 2230 Logically Collective 2231 2232 Input Parameters: 2233 + ts - timestepping context 2234 - flg - `PETSC_TRUE` for fully implicit 2235 2236 Level: intermediate 2237 2238 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 2239 @*/ 2240 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 2241 { 2242 PetscFunctionBegin; 2243 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2244 PetscValidLogicalCollectiveBool(ts, flg, 2); 2245 PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 2246 PetscFunctionReturn(PETSC_SUCCESS); 2247 } 2248 2249 /*@ 2250 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 2251 2252 Logically Collective 2253 2254 Input Parameter: 2255 . ts - timestepping context 2256 2257 Output Parameter: 2258 . flg - `PETSC_TRUE` for fully implicit 2259 2260 Level: intermediate 2261 2262 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 2263 @*/ 2264 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 2265 { 2266 PetscFunctionBegin; 2267 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2268 PetscAssertPointer(flg, 2); 2269 PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 2270 PetscFunctionReturn(PETSC_SUCCESS); 2271 } 2272 2273 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 2274 { 2275 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2276 2277 PetscFunctionBegin; 2278 *arktype = ark->tableau->name; 2279 PetscFunctionReturn(PETSC_SUCCESS); 2280 } 2281 2282 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 2283 { 2284 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2285 PetscBool match; 2286 ARKTableauLink link; 2287 2288 PetscFunctionBegin; 2289 if (ark->tableau) { 2290 PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 2291 if (match) PetscFunctionReturn(PETSC_SUCCESS); 2292 } 2293 for (link = ARKTableauList; link; link = link->next) { 2294 PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 2295 if (match) { 2296 if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 2297 ark->tableau = &link->tab; 2298 if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 2299 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 2300 PetscFunctionReturn(PETSC_SUCCESS); 2301 } 2302 } 2303 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 2304 } 2305 2306 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 2307 { 2308 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2309 2310 PetscFunctionBegin; 2311 ark->imex = (PetscBool)!flg; 2312 PetscFunctionReturn(PETSC_SUCCESS); 2313 } 2314 2315 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 2316 { 2317 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2318 2319 PetscFunctionBegin; 2320 *flg = (PetscBool)!ark->imex; 2321 PetscFunctionReturn(PETSC_SUCCESS); 2322 } 2323 2324 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 2325 { 2326 PetscFunctionBegin; 2327 PetscCall(TSReset_ARKIMEX(ts)); 2328 if (ts->dm) { 2329 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 2330 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 2331 } 2332 PetscCall(PetscFree(ts->data)); 2333 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", NULL)); 2334 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", NULL)); 2335 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 2336 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 2337 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 2338 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 2339 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", NULL)); 2340 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", NULL)); 2341 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_ARKIMEX_FastSlowSplit_C", NULL)); 2342 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_ARKIMEX_FastSlowSplit_C", NULL)); 2343 PetscFunctionReturn(PETSC_SUCCESS); 2344 } 2345 2346 /* ------------------------------------------------------------ */ 2347 /*MC 2348 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 2349 2350 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 2351 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 2352 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 2353 2354 Level: beginner 2355 2356 Notes: 2357 The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 2358 2359 If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 2360 2361 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). 2362 2363 Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 2364 2365 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 2366 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 2367 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 2368 M*/ 2369 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 2370 { 2371 TS_ARKIMEX *ark; 2372 PetscBool dirk; 2373 2374 PetscFunctionBegin; 2375 PetscCall(TSARKIMEXInitializePackage()); 2376 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2377 2378 ts->ops->reset = TSReset_ARKIMEX; 2379 ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 2380 ts->ops->destroy = TSDestroy_ARKIMEX; 2381 ts->ops->view = TSView_ARKIMEX; 2382 ts->ops->load = TSLoad_ARKIMEX; 2383 ts->ops->setup = TSSetUp_ARKIMEX; 2384 ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 2385 ts->ops->step = TSStep_ARKIMEX; 2386 ts->ops->interpolate = TSInterpolate_ARKIMEX; 2387 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 2388 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 2389 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 2390 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 2391 ts->ops->getstages = TSGetStages_ARKIMEX; 2392 ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 2393 2394 ts->usessnes = PETSC_TRUE; 2395 2396 PetscCall(PetscNew(&ark)); 2397 ts->data = (void *)ark; 2398 ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 2399 2400 ark->VecsDeltaLam = NULL; 2401 ark->VecsSensiTemp = NULL; 2402 ark->VecsSensiPTemp = NULL; 2403 2404 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2405 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2406 if (!dirk) { 2407 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 2408 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 2409 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFastSlowSplit_C", TSARKIMEXSetFastSlowSplit_ARKIMEX)); 2410 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFastSlowSplit_C", TSARKIMEXGetFastSlowSplit_ARKIMEX)); 2411 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2412 } 2413 PetscFunctionReturn(PETSC_SUCCESS); 2414 } 2415 2416 /* ------------------------------------------------------------ */ 2417 2418 static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2419 { 2420 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2421 2422 PetscFunctionBegin; 2423 PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2424 PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2425 PetscFunctionReturn(PETSC_SUCCESS); 2426 } 2427 2428 /*@ 2429 TSDIRKSetType - Set the type of `TSDIRK` scheme 2430 2431 Logically Collective 2432 2433 Input Parameters: 2434 + ts - timestepping context 2435 - dirktype - type of `TSDIRK` scheme 2436 2437 Options Database Key: 2438 . -ts_dirkimex_type - set `TSDIRK` scheme type 2439 2440 Level: intermediate 2441 2442 .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2443 @*/ 2444 PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2445 { 2446 PetscFunctionBegin; 2447 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2448 PetscAssertPointer(dirktype, 2); 2449 PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2450 PetscFunctionReturn(PETSC_SUCCESS); 2451 } 2452 2453 /*@ 2454 TSDIRKGetType - Get the type of `TSDIRK` scheme 2455 2456 Logically Collective 2457 2458 Input Parameter: 2459 . ts - timestepping context 2460 2461 Output Parameter: 2462 . dirktype - type of `TSDIRK` scheme 2463 2464 Level: intermediate 2465 2466 .seealso: [](ch_ts), `TSDIRKSetType()` 2467 @*/ 2468 PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2469 { 2470 PetscFunctionBegin; 2471 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2472 PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2473 PetscFunctionReturn(PETSC_SUCCESS); 2474 } 2475 2476 /*MC 2477 TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes. 2478 2479 Level: beginner 2480 2481 Notes: 2482 The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type. 2483 The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with: 2484 + E - whether the method has an explicit first stage 2485 . S - whether the method is single diagonal 2486 . P - order of the advancing method 2487 . Q - order of the embedded method 2488 . S - number of stages 2489 . SA - whether the method is stiffly accurate 2490 . L - whether the method is L-stable 2491 - A - whether the method is A-stable 2492 2493 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2494 M*/ 2495 PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2496 { 2497 PetscFunctionBegin; 2498 PetscCall(TSCreate_ARKIMEX(ts)); 2499 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2500 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2501 PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 2502 PetscFunctionReturn(PETSC_SUCCESS); 2503 } 2504 2505 /*@ 2506 TSARKIMEXSetFastSlowSplit - Use `TSARKIMEX` for solving a fast-slow system 2507 2508 Logically Collective 2509 2510 Input Parameters: 2511 + ts - timestepping context 2512 - fastslow - `PETSC_TRUE` enables the `TSARKIMEX` solver for a fast-slow system where the RHS is split component-wise. 2513 2514 Options Database Key: 2515 . -ts_arkimex_fastslowsplit - <true,false> 2516 2517 Level: intermediate 2518 2519 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetFastSlowSplit()` 2520 @*/ 2521 PetscErrorCode TSARKIMEXSetFastSlowSplit(TS ts, PetscBool fastslow) 2522 { 2523 PetscFunctionBegin; 2524 PetscTryMethod(ts, "TSARKIMEXSetFastSlowSplit_C", (TS, PetscBool), (ts, fastslow)); 2525 PetscFunctionReturn(PETSC_SUCCESS); 2526 } 2527 2528 /*@ 2529 TSARKIMEXGetFastSlowSplit - Gets whether to use `TSARKIMEX` for a fast-slow system 2530 2531 Not Collective 2532 2533 Input Parameter: 2534 . ts - timestepping context 2535 2536 Output Parameter: 2537 . fastslow - `PETSC_TRUE` if `TSARKIMEX` will be used for solving a fast-slow system, `PETSC_FALSE` otherwise 2538 2539 Level: intermediate 2540 2541 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXSetFastSlowSplit()` 2542 @*/ 2543 PetscErrorCode TSARKIMEXGetFastSlowSplit(TS ts, PetscBool *fastslow) 2544 { 2545 PetscFunctionBegin; 2546 PetscUseMethod(ts, "TSARKIMEXGetFastSlowSplit_C", (TS, PetscBool *), (ts, fastslow)); 2547 PetscFunctionReturn(PETSC_SUCCESS); 2548 } 2549