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