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