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