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