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