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