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 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 TSRHSFunction func; 460 461 PetscFunctionBegin; 462 *has = PETSC_FALSE; 463 PetscCall(DMTSGetRHSFunction(ts->dm, &func, NULL)); 464 if (func) *has = PETSC_TRUE; 465 PetscFunctionReturn(PETSC_SUCCESS); 466 } 467 468 /*@C 469 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX` 470 471 Not Collective, but should be called by all processes which will need the schemes to be registered 472 473 Level: advanced 474 475 .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()` 476 @*/ 477 PetscErrorCode TSARKIMEXRegisterAll(void) 478 { 479 PetscFunctionBegin; 480 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 481 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 482 483 #define RC PetscRealConstant 484 #define s2 RC(1.414213562373095048802) /* PetscSqrtReal((PetscReal)2.0) */ 485 #define us2 RC(0.2928932188134524755992) /* 1.0-1.0/PetscSqrtReal((PetscReal)2.0); */ 486 487 /* Diagonally implicit methods */ 488 { 489 /* DIRK212, default of SUNDIALS */ 490 const PetscReal A[2][2] = { 491 {RC(1.0), RC(0.0)}, 492 {RC(-1.0), RC(1.0)} 493 }; 494 const PetscReal b[2] = {RC(0.5), RC(0.5)}; 495 const PetscReal bembed[2] = {RC(1.0), RC(0.0)}; 496 PetscCall(TSDIRKRegister(TSDIRKS212, 2, 2, &A[0][0], b, NULL, bembed, 1, b)); 497 } 498 499 { 500 /* ESDIRK12 from https://arxiv.org/pdf/1803.01613.pdf */ 501 const PetscReal A[2][2] = { 502 {RC(0.0), RC(0.0)}, 503 {RC(0.0), RC(1.0)} 504 }; 505 const PetscReal b[2] = {RC(0.0), RC(1.0)}; 506 const PetscReal bembed[2] = {RC(0.5), RC(0.5)}; 507 PetscCall(TSDIRKRegister(TSDIRKES122SAL, 1, 2, &A[0][0], b, NULL, bembed, 1, b)); 508 } 509 510 { 511 /* ESDIRK213L[2]SA from KC16. 512 TR-BDF2 from Hosea and Shampine 513 ESDIRK23 in https://arxiv.org/pdf/1803.01613.pdf */ 514 const PetscReal A[3][3] = { 515 {RC(0.0), RC(0.0), RC(0.0)}, 516 {us2, us2, RC(0.0)}, 517 {s2 / RC(4.0), s2 / RC(4.0), us2 }, 518 }; 519 const PetscReal b[3] = {s2 / RC(4.0), s2 / RC(4.0), us2}; 520 const PetscReal bembed[3] = {(RC(1.0) - s2 / RC(4.0)) / RC(3.0), (RC(3.0) * s2 / RC(4.0) + RC(1.0)) / RC(3.0), us2 / RC(3.0)}; 521 PetscCall(TSDIRKRegister(TSDIRKES213SAL, 2, 3, &A[0][0], b, NULL, bembed, 1, b)); 522 } 523 524 { 525 #define g RC(0.43586652150845899941601945) 526 #define g2 PetscSqr(g) 527 #define g3 g *g2 528 #define g4 PetscSqr(g2) 529 #define g5 g *g4 530 #define c3 RC(1.0) 531 #define a32 c3 *(c3 - RC(2.0) * g) / (RC(4.0) * g) 532 #define b2 (-RC(2.0) + RC(3.0) * c3 + RC(6.0) * g * (RC(1.0) - c3)) / (RC(12.0) * g * (c3 - RC(2.0) * g)) 533 #define b3 (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) / (RC(3.0) * c3 * (c3 - RC(2.0) * g)) 534 #if 0 535 /* This is for c3 = 3/5 */ 536 #define bh2 \ 537 c3 * (-RC(1.0) + RC(6.0) * g - RC(23.0) * g3 + RC(12.0) * g4 - RC(6.0) * g5) / (RC(4.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) + (RC(3.0) - RC(27.0) * g + RC(68.0) * g2 - RC(55.0) * g3 + RC(21.0) * g4 - RC(6.0) * g5) / (RC(2.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) 538 #define bh3 -g * (-RC(2.0) + RC(21.0) * g - RC(68.0) * g2 + RC(79.0) * g3 - RC(33.0) * g4 + RC(12.0) * g5) / (RC(2.0) * (RC(2.0) * g - c3) * (RC(1.0) - RC(6.0) * g + RC(6.0) * g2)) 539 #define bh4 -RC(3.0) * g2 * (-RC(1.0) + RC(4.0) * g - RC(2.0) * g2 + g3) / (RC(1.0) - RC(6.0) * g + RC(6.0) * g2) 540 #else 541 /* This is for c3 = 1.0 */ 542 #define bh2 a32 543 #define bh3 g 544 #define bh4 RC(0.0) 545 #endif 546 /* ESDIRK3(2I)4L[2]SA from KC16 with c3 = 1.0 */ 547 /* Given by Kvaerno https://link.springer.com/article/10.1023/b:bitn.0000046811.70614.38 */ 548 const PetscReal A[4][4] = { 549 {RC(0.0), RC(0.0), RC(0.0), RC(0.0)}, 550 {g, g, RC(0.0), RC(0.0)}, 551 {c3 - a32 - g, a32, g, RC(0.0)}, 552 {RC(1.0) - b2 - b3 - g, b2, b3, g }, 553 }; 554 const PetscReal b[4] = {RC(1.0) - b2 - b3 - g, b2, b3, g}; 555 const PetscReal bembed[4] = {RC(1.0) - bh2 - bh3 - bh4, bh2, bh3, bh4}; 556 PetscCall(TSDIRKRegister(TSDIRKES324SAL, 3, 4, &A[0][0], b, NULL, bembed, 1, b)); 557 #undef g 558 #undef g2 559 #undef g3 560 #undef c3 561 #undef a32 562 #undef b2 563 #undef b3 564 #undef bh2 565 #undef bh3 566 #undef bh4 567 } 568 569 { 570 /* ESDIRK3(2I)5L[2]SA from KC16 */ 571 const PetscReal A[5][5] = { 572 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 573 {RC(9.0) / RC(40.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0), RC(0.0) }, 574 {RC(19.0) / RC(72.0), RC(14.0) / RC(45.0), RC(9.0) / RC(40.0), RC(0.0), RC(0.0) }, 575 {RC(3337.0) / RC(11520.0), RC(233.0) / RC(720.0), RC(207.0) / RC(1280.0), RC(9.0) / RC(40.0), RC(0.0) }, 576 {RC(7415.0) / RC(34776.0), RC(9920.0) / RC(30429.0), RC(4845.0) / RC(9016.0), -RC(5827.0) / RC(19320.0), RC(9.0) / RC(40.0)}, 577 }; 578 const PetscReal b[5] = {RC(7415.0) / RC(34776.0), RC(9920.0) / RC(30429.0), RC(4845.0) / RC(9016.0), -RC(5827.0) / RC(19320.0), RC(9.0) / RC(40.0)}; 579 const PetscReal bembed[5] = {RC(23705.0) / RC(104328.0), RC(29720.0) / RC(91287.0), RC(4225.0) / RC(9016.0), -RC(69304987.0) / RC(337732920.0), RC(42843.0) / RC(233080.0)}; 580 PetscCall(TSDIRKRegister(TSDIRKES325SAL, 3, 5, &A[0][0], b, NULL, bembed, 1, b)); 581 } 582 583 { 584 // DIRK(6,6)[1]A[(7,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 585 const PetscReal A[7][7] = { 586 {RC(0.303487844706747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 587 {RC(-0.279756492709814), RC(0.500032236020747), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 588 {RC(0.280583215743895), RC(-0.438560061586751), RC(0.217250734515736), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 589 {RC(-0.0677678738539846), RC(0.984312781232293), RC(-0.266720192540149), RC(0.2476680834526), RC(0.0), RC(0.0), RC(0.0) }, 590 {RC(0.125671616147993), RC(-0.995401751002415), RC(0.761333109549059), RC(-0.210281837202208), RC(0.866743712636936), RC(0.0), RC(0.0) }, 591 {RC(-0.368056238801488), RC(-0.999928082701516), RC(0.534734253232519), RC(-0.174856916279082), RC(0.615007160285509), RC(0.696549912132029), RC(0.0) }, 592 {RC(-0.00570546839653984), RC(-0.113110431835656), RC(-0.000965563207671587), RC(-0.000130490084629567), RC(0.00111737736895673), RC(-0.279385587378871), RC(0.618455906845342)} 593 }; 594 const PetscReal b[7] = {RC(0.257561510484877), RC(0.234281287047716), RC(0.126658904241469), RC(0.252363215441784), RC(0.396701083526306), RC(-0.267566000742152), RC(0.0)}; 595 const PetscReal bembed[7] = {RC(0.257561510484945), RC(0.387312822934391), RC(0.126658904241468), RC(0.252363215441784), RC(0.396701083526306), RC(-0.267566000742225), RC(-0.153031535886669)}; 596 PetscCall(TSDIRKRegister(TSDIRK657A, 6, 7, &A[0][0], b, NULL, bembed, 1, b)); 597 } 598 { 599 // ESDIRK(8,6)[2]SA[(8,4)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 600 const PetscReal A[8][8] = { 601 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 602 {RC(0.333222149217725), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 603 {RC(0.0639743773182214), RC(-0.0830330224410214), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 604 {RC(-0.728522201369326), RC(-0.210414479522485), RC(0.532519916559342), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 605 {RC(-0.175135269272067), RC(0.666675582067552), RC(-0.304400907370867), RC(0.656797712445756), RC(0.333222149217725), RC(0.0), RC(0.0), RC(0.0) }, 606 {RC(0.222695802705462), RC(-0.0948971794681061), RC(-0.0234336346686545), RC(-0.45385925012042), RC(0.0283910313826958), RC(0.333222149217725), RC(0.0), RC(0.0) }, 607 {RC(-0.132534078051299), RC(0.702597935004879), RC(-0.433316453128078), RC(0.893717488547587), RC(0.057381454791406), RC(-0.207798411552402), RC(0.333222149217725), RC(0.0) }, 608 {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)} 609 }; 610 const PetscReal b[8] = {RC(0.0802253121418085), RC(0.281196044671022), RC(0.406758926172157), RC(-0.01945708512416), RC(-0.41785600088526), RC(0.0545342658870322), RC(0.281376387919675), RC(0.333222149217725)}; 611 const PetscReal bembed[8] = {RC(0.0), RC(0.292331064554014), RC(0.409676102283681), RC(-0.002094718084982), RC(-0.282771520835975), RC(0.113862336644901), RC(0.181973572260693), RC(0.287023163177669)}; 612 PetscCall(TSDIRKRegister(TSDIRKES648SA, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 613 } 614 { 615 // DIRK(8,6)[1]SAL[(8,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 616 const PetscReal A[8][8] = { 617 {RC(0.477264457385826), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 618 {RC(-0.197052588415002), RC(0.476363428459584), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 619 {RC(-0.0347674430372966), RC(0.633051807335483), RC(0.193634310075028), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 620 {RC(0.0967797668578702), RC(-0.193533526466535), RC(-0.000207622945800473), RC(0.159572204849431), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 621 {RC(0.162527231819875), RC(-0.249672513547382), RC(-0.0459079972041795), RC(0.36579476400859), RC(0.255752838307699), RC(0.0), RC(0.0), RC(0.0) }, 622 {RC(-0.00707603197171262), RC(0.846299854860295), RC(0.344020016925018), RC(-0.0720926054548865), RC(-0.215492331980875), RC(0.104341097622161), RC(0.0), RC(0.0) }, 623 {RC(0.00176857935179744), RC(0.0779960013127515), RC(0.303333277564557), RC(0.213160806732836), RC(0.351769320319038), RC(-0.381545894386538), RC(0.433517909105558), RC(0.0) }, 624 {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)} 625 }; 626 const PetscReal b[8] = {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.120953626732831), RC(-0.0802678473399899), RC(0.264667545261832)}; 627 const PetscReal bembed[8] = {RC(0.0), RC(0.22732353410559), RC(0.308415837980118), RC(0.157263419573007), RC(0.243551137152275), RC(-0.103483943222765), RC(-0.0103721771642262), RC(0.177302191576001)}; 628 PetscCall(TSDIRKRegister(TSDIRK658A, 6, 8, &A[0][0], b, NULL, bembed, 1, b)); 629 } 630 { 631 // SDIRK(9,6)[1]SAL[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 632 const PetscReal A[9][9] = { 633 {RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 634 {RC(-0.0903514856119419), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 635 {RC(0.172952039138937), RC(-0.35365501036282), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 636 {RC(0.511999875919193), RC(0.0289640332201925), RC(-0.0144030945657094), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 637 {RC(0.00465303495506782), RC(-0.075635818766597), RC(0.217273030786712), RC(-0.0206519428725472), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 638 {RC(0.896145501762472), RC(0.139267327700498), RC(-0.186920979752805), RC(0.0672971012371723), RC(-0.350891963442176), RC(0.218127781944908), RC(0.0), RC(0.0), RC(0.0) }, 639 {RC(0.552959701885751), RC(-0.439360579793662), RC(0.333704002325091), RC(-0.0339426520778416), RC(-0.151947445912595), RC(0.0213825661026943), RC(0.218127781944908), RC(0.0), RC(0.0) }, 640 {RC(0.631360374036476), RC(0.724733619641466), RC(-0.432170625425258), RC(0.598611382182477), RC(-0.709087197034345), RC(-0.483986685696934), RC(0.378391562905131), RC(0.218127781944908), RC(0.0) }, 641 {RC(0.0), RC(-0.15504452530869), RC(0.194518478660789), RC(0.63515640279203), RC(0.81172278664173), RC(0.110736108691585), RC(-0.495304692414479), RC(-0.319912341007872), RC(0.218127781944908)} 642 }; 643 const PetscReal b[9] = {RC(0.0), RC(-0.15504452530869), RC(0.194518478660789), RC(0.63515640279203), RC(0.81172278664173), RC(0.110736108691585), RC(-0.495304692414479), RC(-0.319912341007872), RC(0.218127781944908)}; 644 const PetscReal bembed[9] = {RC(3.62671059311602e-16), RC(0.0736615558278942), RC(0.103527397262229), RC(1.00247481935499), RC(0.361377289250057), RC(-0.785425929961365), RC(-0.0170499047960784), RC(0.296321252214769), RC(-0.0348864791524953)}; 645 PetscCall(TSDIRKRegister(TSDIRKS659A, 6, 9, &A[0][0], b, NULL, bembed, 1, b)); 646 } 647 { 648 // DIRK(10,7)[1]SAL[(10,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 649 const PetscReal A[10][10] = { 650 {RC(0.233704632125264), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 651 {RC(-0.0739324813149407), RC(0.200056838146104), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 652 {RC(0.0943790344044812), RC(0.264056067701605), RC(0.133245202456465), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 653 {RC(0.269084810601201), RC(-0.503479002548384), RC(-0.00486736469695022), RC(0.251518716213569), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 654 {RC(0.145665801918423), RC(0.204983170463176), RC(0.407154634069484), RC(-0.0121039135200389), RC(0.190243622486334), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 655 {RC(0.985450198547345), RC(0.806942652811456), RC(-0.808130934167263), RC(-0.669035819439391), RC(0.0269384406756128), RC(0.462144080607327), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 656 {RC(0.163902957809563), RC(0.228315094960095), RC(0.0745971021260249), RC(0.000509793400156559), RC(0.0166533681378294), RC(-0.0229383879045797), RC(0.103505486637336), RC(0.0), RC(0.0), RC(0.0) }, 657 {RC(-0.162694156858437), RC(0.0453478837428434), RC(0.997443481211424), RC(0.200251514941093), RC(-0.000161755458839048), RC(-0.0848134335980281), RC(-0.36438666566666), RC(0.158604420136055), RC(0.0), RC(0.0) }, 658 {RC(0.200733156477425), RC(0.239686443444433), RC(0.303837014418929), RC(-0.0534390596279896), RC(0.0314067599640569), RC(-0.00764032790448536), RC(0.0609191260198661), RC(-0.0736319201590642), RC(0.204602530607021), RC(0.0) }, 659 {RC(0.0), RC(0.235563761744267), RC(0.658651488684319), RC(0.0308877804992098), RC(-0.906514945595336), RC(-0.0248488551739974), RC(-0.309967582365257), RC(0.191663316925525), RC(0.923933712199542), RC(0.200631323081727)} 660 }; 661 const PetscReal b[10] = {RC(0.0), RC(0.235563761744267), RC(0.658651488684319), RC(0.0308877804992098), RC(-0.906514945595336), RC(-0.0248488551739974), RC(-0.309967582365257), RC(0.191663316925525), RC(0.923933712199542), RC(0.200631323081727)}; 662 const PetscReal bembed[10] = 663 {RC(0.0), RC(0.222929376486581), RC(0.950668440138169), RC(0.0342694607044032), RC(0.362875840545746), RC(0.223572979288581), RC(-0.764361723526727), RC(0.563476909230026), RC(-0.690896961894185), RC(0.0974656790270323)}; 664 PetscCall(TSDIRKRegister(TSDIRK7510SAL, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 665 } 666 { 667 // ESDIRK(10,7)[2]SA[(10,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 668 const PetscReal A[10][10] = { 669 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 670 {RC(0.210055790203419), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 671 {RC(0.255781739921086), RC(0.239850916980976), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 672 {RC(0.286789624880437), RC(0.230494748834778), RC(0.263925149885491), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 673 {RC(-0.0219118128774335), RC(0.897684380345907), RC(-0.657954605498907), RC(0.124962304722633), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 674 {RC(-0.065614879584776), RC(-0.0565630711859497), RC(0.0254881105065311), RC(-0.00368981790650006), RC(-0.0115178258446329), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 675 {RC(0.399860851232098), RC(0.915588469718705), RC(-0.0758429094934412), RC(-0.263369154872759), RC(0.719687583564526), RC(-0.787410407015369), RC(0.210055790203419), RC(0.0), RC(0.0), RC(0.0) }, 676 {RC(0.51693616104628), RC(1.00000540846973), RC(-0.0485110663289207), RC(-0.315208041581942), RC(0.749742806451587), RC(-0.990975090921248), RC(0.0159279583407308), RC(0.210055790203419), RC(0.0), RC(0.0) }, 677 {RC(-0.0303062129076945), RC(-0.297035174659034), RC(0.184724697462164), RC(-0.0351876079516183), RC(-0.00324668230690761), RC(0.216151004053531), RC(-0.126676252098317), RC(0.114040254365262), RC(0.210055790203419), RC(0.0) }, 678 {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)} 679 }; 680 const PetscReal b[10] = {RC(0.0705997961586714), RC(-0.0281516061956374), RC(0.314600470734633), RC(-0.0907057557963371), RC(0.168078953957742), 681 RC(-0.00655694984590575), RC(0.0505384497804303), RC(-0.0569572058725042), RC(0.368498056875488), RC(0.210055790203419)}; 682 const PetscReal bembed[10] = {RC(-0.015494246543626), RC(0.167657963820093), RC(0.269858958144236), RC(-0.0443258997755156), RC(0.150049236875266), 683 RC(0.259452082755846), RC(0.244624573502521), RC(-0.215528446920284), RC(0.0487601760292619), RC(0.134945602112201)}; 684 PetscCall(TSDIRKRegister(TSDIRKES7510SA, 7, 10, &A[0][0], b, NULL, bembed, 1, b)); 685 } 686 { 687 // DIRK(9,7)[1]A[(9,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 688 const PetscReal A[9][9] = { 689 {RC(0.179877789855839), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 690 {RC(-0.100405844885157), RC(0.214948590644819), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 691 {RC(0.112251360198995), RC(-0.206162139150298), RC(0.125159642941958), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 692 {RC(-0.0335164000768257), RC(0.999942349946143), RC(-0.491470853833294), RC(0.19820086325566), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 693 {RC(-0.0417345265478321), RC(0.187864510308215), RC(0.0533789224305102), RC(-0.00822060284862916), RC(0.127670843671646), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 694 {RC(-0.0278257925239257), RC(0.600979340683382), RC(-0.242632273241134), RC(-0.11318753652081), RC(0.164326917632931), RC(0.284116597781395), RC(0.0), RC(0.0), RC(0.0) }, 695 {RC(0.041465583858922), RC(0.429657872601836), RC(-0.381323410582524), RC(0.391934277498434), RC(-0.245918275501241), RC(-0.35960669741231), RC(0.184000022289158), RC(0.0), RC(0.0) }, 696 {RC(-0.105565651574538), RC(-0.0557833155018609), RC(0.358967568942643), RC(-0.13489263413921), RC(0.129553247260677), RC(0.0992493795371489), RC(-0.15716610563461), RC(0.17918862279814), RC(0.0) }, 697 {RC(0.00439696079965225), RC(0.960250486570491), RC(0.143558372286706), RC(0.0819015241056593), RC(0.999562318563625), RC(0.325203439314358), RC(-0.679013149331228), RC(-0.990589559837246), RC(0.0773648037639896)} 698 }; 699 700 const PetscReal b[9] = {RC(0.0), RC(0.179291520437966), RC(0.115310295273026), RC(-0.857943261453138), RC(0.654911318641998), RC(1.18713633508094), RC(-0.0949482361570542), RC(-0.37661430946407), RC(0.19285633764033)}; 701 const PetscReal bembed[9] = {RC(0.0), RC(0.1897135479408), RC(0.127461414808862), RC(-0.835810807663404), RC(0.665114177777166), RC(1.16481046518346), RC(-0.11661858889792), RC(-0.387303251022099), RC(0.192633041873135)}; 702 PetscCall(TSDIRKRegister(TSDIRK759A, 7, 9, &A[0][0], b, NULL, bembed, 1, b)); 703 } 704 { 705 // SDIRK(11,7)[1]SAL[(11,5)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 706 const PetscReal A[11][11] = { 707 {RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 708 {RC(-0.082947368165267), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 709 {RC(0.483452690540751), RC(0.0), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 710 {RC(0.771076453481321), RC(-0.22936926341842), RC(0.289733373208823), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 711 {RC(0.0329683054968892), RC(-0.162397421903366), RC(0.000951777538562805), RC(0.0), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 712 {RC(0.265888743485945), RC(0.606743151103931), RC(0.173443800537369), RC(-0.0433968261546912), RC(-0.385211017224481), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 713 {RC(0.220662294551146), RC(-0.0465078507657608), RC(-0.0333111995282464), RC(0.011801580836998), RC(0.169480801030105), RC(-0.0167974432139385), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 714 {RC(0.323099728365267), RC(0.0288371831672575), RC(-0.0543404318773196), RC(0.0137765831431662), RC(0.0516799019060702), RC(-0.0421359763835713), RC(0.181297932037826), RC(0.200252661187742), RC(0.0), RC(0.0), RC(0.0) }, 715 {RC(-0.164226696476538), RC(0.187552004946792), RC(0.0628674420973025), RC(-0.0108886582703428), RC(-0.0117628641717889), RC(0.0432176880867965), RC(-0.0315206836275473), RC(-0.0846007021638797), RC(0.200252661187742), RC(0.0), RC(0.0) }, 716 {RC(0.651428598623771), RC(-0.10208078475356), RC(0.198305701801888), RC(-0.0117354096673789), RC(-0.0440385966743686), RC(-0.0358364455795087), RC(-0.0075408087654097), RC(0.160320941654639), RC(0.017940248694499), RC(0.200252661187742), RC(0.0) }, 717 {RC(0.0), RC(-0.266259448580236), RC(-0.615982357748271), RC(0.561474126687165), RC(0.266911112787025), RC(0.219775952207137), RC(0.387847665451514), RC(0.612483137773236), RC(0.330027015806089), RC(-0.6965298655714), RC(0.200252661187742)} 718 }; 719 const PetscReal b[11] = 720 {RC(0.0), RC(-0.266259448580236), RC(-0.615982357748271), RC(0.561474126687165), RC(0.266911112787025), RC(0.219775952207137), RC(0.387847665451514), RC(0.612483137773236), RC(0.330027015806089), RC(-0.6965298655714), RC(0.200252661187742)}; 721 const PetscReal bembed[11] = 722 {RC(0.0), RC(0.180185524442613), RC(-0.628869710835338), RC(0.186185675988647), RC(0.0484716652630425), RC(0.203927720607141), RC(0.44041662512573), RC(0.615710527731245), RC(0.0689648839032607), RC(-0.253599870605903), RC(0.138606958379488)}; 723 PetscCall(TSDIRKRegister(TSDIRKS7511SAL, 7, 11, &A[0][0], b, NULL, bembed, 1, b)); 724 } 725 { 726 // DIRK(13,8)[1]A[(14,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 727 const PetscReal A[14][14] = { 728 {RC(0.421050745442291), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 729 {RC(-0.0761079419591268), RC(0.264353986580857), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 730 {RC(0.0727106904170694), RC(-0.204265976977285), RC(0.181608196544136), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 731 {RC(0.55763054816611), RC(-0.409773579543499), RC(0.510926516886944), RC(0.259892204518476), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 732 {RC(0.0228083864844437), RC(-0.445569051836454), RC(-0.0915242778636248), RC(0.00450055909321655), RC(0.6397807199983), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 733 {RC(-0.135945849505152), RC(0.0946509646963754), RC(-0.236110197279175), RC(0.00318944206456517), RC(0.255453021028118), RC(0.174805219173446), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 734 {RC(-0.147960260670772), RC(-0.402188192230535), RC(-0.703014530043888), RC(0.00941974677418186), RC(0.885747111289207), RC(0.261314066449028), RC(0.16307697503668), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 735 {RC(0.165597241042244), RC(0.824182962188923), RC(-0.0280136160783609), RC(0.282372386631758), RC(-0.957721354131182), RC(0.489439550159977), RC(0.170094415598103), RC(0.0522519785718563), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 736 {RC(0.0335292011495618), RC(0.575750388029166), RC(0.223289855356637), RC(-0.00317458833242804), RC(-0.112890382135193), RC(-0.419809267954284), RC(0.0466136902102104), RC(-0.00115413813041085), RC(0.109685363692383), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 737 {RC(-0.0512616878252355), RC(0.699261265830807), RC(-0.117939611738769), RC(0.0021745241931243), RC(-0.00932826702640947), RC(-0.267575057469428), RC(0.126949139814065), RC(0.00330353204502163), RC(0.185949445053766), RC(0.0938215615963721), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 738 {RC(-0.106521517960343), RC(0.41835889096168), RC(0.353585905881916), RC(-0.0746474161579599), RC(-0.015450626460289), RC(-0.46224659192275), RC(-0.0576406327329181), RC(-0.00712066942504018), RC(0.377776558014452), RC(0.36890054338294), RC(0.0618488746331837), RC(0.0), RC(0.0), RC(0.0) }, 739 {RC(-0.163079104890997), RC(0.644561721693806), RC(0.636968661639572), RC(-0.122346720085377), RC(-0.333062564990312), RC(-0.3054226490478), RC(-0.357820712828352), RC(-0.0125510510334706), RC(0.371263681186311), RC(0.371979640363694), RC(0.0531090658708968), RC(0.0518279459132049), RC(0.0), RC(0.0) }, 740 {RC(0.579993784455521), RC(-0.188833728676494), RC(0.999975696843775), RC(0.0572810855901161), RC(-0.264374735003671), RC(0.165091739976854), RC(-0.546675809010452), RC(-0.0283821822291982), RC(-0.102639860418374), RC(-0.0343251040446405), RC(0.4762598462591), RC(-0.304153104931261), RC(0.0953911855943621), RC(0.0) }, 741 {RC(0.0848552694007844), RC(0.287193912340074), RC(0.543683503004232), RC(-0.081311059300692), RC(-0.0328661289388557), RC(-0.323456834372922), RC(-0.240378871658975), RC(-0.0189913019930369), RC(0.220663114082036), RC(0.253029984360864), RC(0.252011799370563), RC(-0.154882222605423), RC(0.0315202264687415), RC(0.0514095812104714)} 742 }; 743 const PetscReal b[14] = {RC(0.0), RC(0.516650324205117), RC(0.0773227217357826), RC(-0.12474204666975), RC(-0.0241052115180679), RC(-0.325821145180359), RC(0.0907237460123951), RC(0.0459271880596652), RC(0.221012259404702), RC(0.235510906761942), RC(0.491109674204385), RC(-0.323506525837343), RC(0.119918108821531), RC(0.0)}; 744 const PetscReal bembed[14] = {RC(2.32345691433618e-16), RC(0.499150900944401), RC(0.080991997189243), RC(-0.0359440417166322), RC(-0.0258910397441454), RC(-0.304540350278636), RC(0.0836627473632563), 745 RC(0.0417664613347638), RC(0.223636394275293), RC(0.231569156867596), RC(0.240526201277663), RC(-0.222933582911926), RC(-0.0111479879597561), RC(0.19915314335888)}; 746 PetscCall(TSDIRKRegister(TSDIRK8614A, 8, 14, &A[0][0], b, NULL, bembed, 1, b)); 747 } 748 { 749 // DIRK(15,8)[1]SAL[(16,6)A] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 750 const PetscReal A[16][16] = { 751 {RC(0.498904981271193), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 752 {RC(-0.303806037341816), RC(0.886299445992379), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 753 {RC(-0.581440223471476), RC(0.371003719460259), RC(0.43844717752802), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 754 {RC(0.531852638870051), RC(-0.339363014907108), RC(0.422373239795441), RC(0.223854203543397), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 755 {RC(0.118517891868867), RC(-0.0756235584174296), RC(-0.0864284870668712), RC(0.000536692838658312), RC(0.10101418329932), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 756 {RC(0.218733626116401), RC(-0.139568928299635), RC(0.30473612813488), RC(0.00354038623073564), RC(0.0932085751160559), RC(0.140161806097591), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 757 {RC(0.0692944686081835), RC(-0.0442152168939502), RC(-0.0903375348855603), RC(0.00259030241156141), RC(0.204514233679515), RC(-0.0245383758960002), RC(0.199289437094059), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 758 {RC(0.990640016505571), RC(-0.632104756315967), RC(0.856971425234221), RC(0.174494099232246), RC(-0.113715829680145), RC(-0.151494045307366), RC(-0.438268629569005), RC(0.120578398912139), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 759 {RC(-0.099415677713136), RC(0.211832014309207), RC(-0.245998265866888), RC(-0.182249672235861), RC(0.167897635713799), RC(0.212850335030069), RC(-0.391739299440123), RC(-0.0118718506876767), RC(0.526293701659093), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 760 {RC(0.383983914845461), RC(-0.245011361219604), RC(0.46717278554955), RC(-0.0361272447593202), RC(0.0742234660511333), RC(-0.0474816271948766), RC(-0.229859978525756), RC(0.0516283729206322), RC(0.0), RC(0.193823890777594), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 761 {RC(0.0967855003180134), RC(-0.0481037037916184), RC(0.191268138832434), RC(0.234977164564126), RC(0.0620265921753097), RC(0.403432826534738), RC(0.152403846687238), RC(-0.118420429237746), RC(0.0582141598685892), RC(-0.13924540906863), RC(0.106661313117545), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 762 {RC(0.133941307432055), RC(-0.0722076602896254), RC(0.217086297689275), RC(0.00495499602192887), RC(0.0306090174933995), RC(0.26483526755746), RC(0.204442440745605), RC(0.196883395136708), RC(0.056527012583996), RC(-0.150216381356784), RC(-0.217209415757333), RC(0.330353722743315), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 763 {RC(0.157014274561299), RC(-0.0883810256381874), RC(0.117193033885034), RC(-0.0362304243769466), RC(0.0169030211466111), RC(-0.169835753576141), RC(0.399749979234113), RC(0.31806704093008), RC(0.050340008347693), RC(0.120284837472214), RC(-0.235313193645423), RC(0.232488522208926), RC(0.117719679450729), RC(0.0), RC(0.0), RC(0.0) }, 764 {RC(0.00276453816875833), RC(-0.00366028255231782), RC(-0.331078914515559), RC(0.623377549031949), RC(0.167618142989491), RC(0.0748467945312516), RC(0.797629286699677), RC(-0.390714256799583), RC(-0.00808553925131555), RC(0.014840324980952), RC(-0.0856180410248133), RC(0.602943304937827), RC(-0.5771359338496), RC(0.112273026653282), RC(0.0), RC(0.0) }, 765 {RC(0.0), RC(0.0), RC(0.085283971980307), RC(0.51334393454179), RC(0.144355978013514), RC(0.255379109487853), RC(0.225075750790524), RC(-0.343241323394982), RC(0.0), RC(0.0798250392218852), RC(0.0528824734082655), RC(-0.0830350888900362), RC(0.022567388707279), RC(-0.0592631119040204), RC(0.106825878037621), RC(0.0) }, 766 {RC(0.173784481207652), RC(-0.110887906116241), RC(0.190052513365204), RC(-0.0688345422674029), RC(0.10326505079603), RC(0.267127097115219), RC(0.141703423176897), RC(0.0117966866651728), RC(-6.65650091812762e-15), RC(-0.0213725083662519), RC(-0.00931148598712566), RC(-0.10007679077114), RC(0.123471797451553), RC(0.00203684241073055), RC(-0.0294320891781173), RC(0.195746619921528)} 767 }; 768 const PetscReal b[16] = {RC(0.0), RC(0.0), RC(0.085283971980307), RC(0.51334393454179), RC(0.144355978013514), RC(0.255379109487853), RC(0.225075750790524), RC(-0.343241323394982), RC(0.0), RC(0.0798250392218852), RC(0.0528824734082655), RC(-0.0830350888900362), RC(0.022567388707279), RC(-0.0592631119040204), RC(0.106825878037621), RC(0.0)}; 769 const PetscReal bembed[16] = {RC(-1.31988512519898e-15), RC(7.53606601764004e-16), RC(0.0886789133915965), RC(0.0968726531622137), RC(0.143815375874267), RC(0.335214773313601), RC(0.221862366978063), RC(-0.147408947987273), 770 RC(4.16297599203445e-16), RC(0.000727276166520566), RC(-0.00284892677941246), RC(0.00512492274297611), RC(-0.000275595071215218), RC(0.0136014719350733), RC(0.0165190013607726), RC(0.228116714912817)}; 771 PetscCall(TSDIRKRegister(TSDIRK8616SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 772 } 773 { 774 // ESDIRK(16,8)[2]SAL[(16,5)] from https://github.com/yousefalamri55/High_Order_DIRK_Methods_Coeffs 775 const PetscReal A[16][16] = { 776 {RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 777 {RC(0.117318819358521), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 778 {RC(0.0557014605974616), RC(0.385525646638742), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 779 {RC(0.063493276428895), RC(0.373556126263681), RC(0.0082994166438953), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 780 {RC(0.0961351856230088), RC(0.335558324517178), RC(0.207077765910132), RC(-0.0581917140797146), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 781 {RC(0.0497669214238319), RC(0.384288616546039), RC(0.0821728117583936), RC(0.120337007107103), RC(0.202262782645888), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 782 {RC(0.00626710666809847), RC(0.496491452640725), RC(-0.111303249827358), RC(0.170478821683603), RC(0.166517073971103), RC(-0.0328669811542241), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 783 {RC(0.0463439767281591), RC(0.00306724391019652), RC(-0.00816305222386205), RC(-0.0353302599538294), RC(0.0139313601702569), RC(-0.00992014507967429), RC(0.0210087909090165), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 784 {RC(0.111574049232048), RC(0.467639166482209), RC(0.237773114804619), RC(0.0798895699267508), RC(0.109580615914593), RC(0.0307353103825936), RC(-0.0404391509541147), RC(-0.16942110744293), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 785 {RC(-0.0107072484863877), RC(-0.231376703354252), RC(0.017541113036611), RC(0.144871527682418), RC(-0.041855459769806), RC(0.0841832168332261), RC(-0.0850020937282192), RC(0.486170343825899), RC(-0.0526717116822739), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 786 {RC(-0.0142238262314935), RC(0.14752923682514), RC(0.238235830732566), RC(0.037950291904103), RC(0.252075123381518), RC(0.0474266904224567), RC(-0.00363139069342027), RC(0.274081442388563), RC(-0.0599166970745255), RC(-0.0527138812389185), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 787 {RC(-0.11837020183211), RC(-0.635712481821264), RC(0.239738832602538), RC(0.330058936651707), RC(-0.325784087988237), RC(-0.0506514314589253), RC(-0.281914404487009), RC(0.852596345144291), RC(0.651444614298805), RC(-0.103476387303591), RC(-0.354835880209975), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0), RC(0.0) }, 788 {RC(-0.00458164025442349), RC(0.296219694015248), RC(0.322146049419995), RC(0.15917778285238), RC(0.284864871688843), RC(0.185509526463076), RC(-0.0784621067883274), RC(0.166312223692047), RC(-0.284152486083397), RC(-0.357125104338944), RC(0.078437074055306), RC(0.0884129667114481), RC(0.117318819358521), RC(0.0), RC(0.0), RC(0.0) }, 789 {RC(-0.0545561913848106), RC(0.675785423442753), RC(0.423066443201941), RC(-0.000165300126841193), RC(0.104252994793763), RC(-0.105763019303021), RC(-0.15988308809318), RC(0.0515050001032011), RC(0.56013979290924), RC(-0.45781539708603), RC(-0.255870699752664), RC(0.026960254296416), RC(-0.0721245985053681), RC(0.117318819358521), RC(0.0), RC(0.0) }, 790 {RC(0.0649253995775223), RC(-0.0216056457922249), RC(-0.073738139377975), RC(0.0931033310077225), RC(-0.0194339577299149), RC(-0.0879623837313009), RC(0.057125517179467), RC(0.205120850488097), RC(0.132576503537441), RC(0.489416890627328), RC(-0.1106765720501), RC(-0.081038793996096), RC(0.0606031613503788), RC(-0.00241467937442272), RC(0.117318819358521), RC(0.0) }, 791 {RC(0.0459979286336779), RC(0.0780075394482806), RC(0.015021874148058), RC(0.195180277284195), RC(-0.00246643310153235), RC(0.0473977117068314), RC(-0.0682773558610363), RC(0.19568019123878), RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)} 792 }; 793 const PetscReal b[16] = {RC(0.0459979286336779), RC(0.0780075394482806), RC(0.015021874148058), RC(0.195180277284195), RC(-0.00246643310153235), RC(0.0473977117068314), RC(-0.0682773558610363), RC(0.19568019123878), 794 RC(-0.0876765449323747), RC(0.177874852409192), RC(-0.337519251582222), RC(-0.0123255553640736), RC(0.311573291192553), RC(0.0458604327754991), RC(0.278352222645651), RC(0.117318819358521)}; 795 const PetscReal bembed[16] = {RC(0.0603373529853206), RC(0.175453809423998), RC(0.0537707777611352), RC(0.195309248607308), RC(0.0135893741970232), RC(-0.0221160259296707), RC(-0.00726526156430691), RC(0.102961059369124), 796 RC(0.000900215457460583), RC(0.0547959465692338), RC(-0.334995726863153), RC(0.0464409662093384), RC(0.301388101652194), RC(0.00524851570622031), RC(0.229538601845236), RC(0.124643044573514)}; 797 PetscCall(TSDIRKRegister(TSDIRKES8516SAL, 8, 16, &A[0][0], b, NULL, bembed, 1, b)); 798 } 799 800 /* Additive methods */ 801 { 802 const PetscReal A[3][3] = { 803 {0.0, 0.0, 0.0}, 804 {0.0, 0.0, 0.0}, 805 {0.0, 0.5, 0.0} 806 }; 807 const PetscReal At[3][3] = { 808 {1.0, 0.0, 0.0}, 809 {0.0, 0.5, 0.0}, 810 {0.0, 0.5, 0.5} 811 }; 812 const PetscReal b[3] = {0.0, 0.5, 0.5}; 813 const PetscReal bembedt[3] = {1.0, 0.0, 0.0}; 814 PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 815 } 816 { 817 const PetscReal A[2][2] = { 818 {0.0, 0.0}, 819 {0.5, 0.0} 820 }; 821 const PetscReal At[2][2] = { 822 {0.0, 0.0}, 823 {0.0, 0.5} 824 }; 825 const PetscReal b[2] = {0.0, 1.0}; 826 const PetscReal bembedt[2] = {0.5, 0.5}; 827 /* binterpt[2][2] = {{1.0,-1.0},{0.0,1.0}}; second order dense output has poor stability properties and hence it is not currently in use */ 828 PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 829 } 830 { 831 const PetscReal A[2][2] = { 832 {0.0, 0.0}, 833 {1.0, 0.0} 834 }; 835 const PetscReal At[2][2] = { 836 {0.0, 0.0}, 837 {0.5, 0.5} 838 }; 839 const PetscReal b[2] = {0.5, 0.5}; 840 const PetscReal bembedt[2] = {0.0, 1.0}; 841 /* binterpt[2][2] = {{1.0,-0.5},{0.0,0.5}} second order dense output has poor stability properties and hence it is not currently in use */ 842 PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 843 } 844 { 845 const PetscReal A[2][2] = { 846 {0.0, 0.0}, 847 {1.0, 0.0} 848 }; 849 const PetscReal At[2][2] = { 850 {us2, 0.0}, 851 {1.0 - 2.0 * us2, us2} 852 }; 853 const PetscReal b[2] = {0.5, 0.5}; 854 const PetscReal bembedt[2] = {0.0, 1.0}; 855 const PetscReal binterpt[2][2] = { 856 {(us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))}, 857 {1 - (us2 - 1.0) / (2.0 * us2 - 1.0), -1 / (2.0 * (1.0 - 2.0 * us2))} 858 }; 859 const PetscReal binterp[2][2] = { 860 {1.0, -0.5}, 861 {0.0, 0.5 } 862 }; 863 PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 864 } 865 { 866 const PetscReal A[3][3] = { 867 {0, 0, 0}, 868 {2 - s2, 0, 0}, 869 {0.5, 0.5, 0} 870 }; 871 const PetscReal At[3][3] = { 872 {0, 0, 0 }, 873 {1 - 1 / s2, 1 - 1 / s2, 0 }, 874 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 875 }; 876 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 877 const PetscReal binterpt[3][2] = { 878 {1.0 / s2, -1.0 / (2.0 * s2)}, 879 {1.0 / s2, -1.0 / (2.0 * s2)}, 880 {1.0 - s2, 1.0 / s2 } 881 }; 882 PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 883 } 884 { 885 const PetscReal A[3][3] = { 886 {0, 0, 0}, 887 {2 - s2, 0, 0}, 888 {0.75, 0.25, 0} 889 }; 890 const PetscReal At[3][3] = { 891 {0, 0, 0 }, 892 {1 - 1 / s2, 1 - 1 / s2, 0 }, 893 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 894 }; 895 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 896 const PetscReal binterpt[3][2] = { 897 {1.0 / s2, -1.0 / (2.0 * s2)}, 898 {1.0 / s2, -1.0 / (2.0 * s2)}, 899 {1.0 - s2, 1.0 / s2 } 900 }; 901 PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 902 } 903 { /* Optimal for linear implicit part */ 904 const PetscReal A[3][3] = { 905 {0, 0, 0}, 906 {2 - s2, 0, 0}, 907 {(3 - 2 * s2) / 6, (3 + 2 * s2) / 6, 0} 908 }; 909 const PetscReal At[3][3] = { 910 {0, 0, 0 }, 911 {1 - 1 / s2, 1 - 1 / s2, 0 }, 912 {1 / (2 * s2), 1 / (2 * s2), 1 - 1 / s2} 913 }; 914 const PetscReal bembedt[3] = {(4. - s2) / 8., (4. - s2) / 8., 1 / (2. * s2)}; 915 const PetscReal binterpt[3][2] = { 916 {1.0 / s2, -1.0 / (2.0 * s2)}, 917 {1.0 / s2, -1.0 / (2.0 * s2)}, 918 {1.0 - s2, 1.0 / s2 } 919 }; 920 PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 921 } 922 { /* Optimal for linear implicit part */ 923 const PetscReal A[3][3] = { 924 {0, 0, 0}, 925 {0.5, 0, 0}, 926 {0.5, 0.5, 0} 927 }; 928 const PetscReal At[3][3] = { 929 {0.25, 0, 0 }, 930 {0, 0.25, 0 }, 931 {1. / 3, 1. / 3, 1. / 3} 932 }; 933 PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 934 } 935 { 936 const PetscReal A[4][4] = { 937 {0, 0, 0, 0}, 938 {1767732205903. / 2027836641118., 0, 0, 0}, 939 {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 940 {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 941 }; 942 const PetscReal At[4][4] = { 943 {0, 0, 0, 0 }, 944 {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0 }, 945 {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0 }, 946 {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.} 947 }; 948 const PetscReal bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}; 949 const PetscReal binterpt[4][2] = { 950 {4655552711362. / 22874653954995., -215264564351. / 13552729205753. }, 951 {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119. }, 952 {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, 953 {584795268549. / 6622622206610., 2508943948391. / 7218656332882. } 954 }; 955 PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 956 } 957 { 958 const PetscReal A[5][5] = { 959 {0, 0, 0, 0, 0}, 960 {1. / 2, 0, 0, 0, 0}, 961 {11. / 18, 1. / 18, 0, 0, 0}, 962 {5. / 6, -5. / 6, .5, 0, 0}, 963 {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 964 }; 965 const PetscReal At[5][5] = { 966 {0, 0, 0, 0, 0 }, 967 {0, 1. / 2, 0, 0, 0 }, 968 {0, 1. / 6, 1. / 2, 0, 0 }, 969 {0, -1. / 2, 1. / 2, 1. / 2, 0 }, 970 {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2} 971 }; 972 PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 973 } 974 { 975 const PetscReal A[5][5] = { 976 {0, 0, 0, 0, 0}, 977 {1, 0, 0, 0, 0}, 978 {4. / 9, 2. / 9, 0, 0, 0}, 979 {1. / 4, 0, 3. / 4, 0, 0}, 980 {1. / 4, 0, 3. / 5, 0, 0} 981 }; 982 const PetscReal At[5][5] = { 983 {0, 0, 0, 0, 0 }, 984 {.5, .5, 0, 0, 0 }, 985 {5. / 18, -1. / 9, .5, 0, 0 }, 986 {.5, 0, 0, .5, 0 }, 987 {.25, 0, .75, -.5, .5} 988 }; 989 PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 990 } 991 { 992 const PetscReal A[6][6] = { 993 {0, 0, 0, 0, 0, 0}, 994 {1. / 2, 0, 0, 0, 0, 0}, 995 {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 996 {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 997 {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 998 {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 999 }; 1000 const PetscReal At[6][6] = { 1001 {0, 0, 0, 0, 0, 0 }, 1002 {1. / 4, 1. / 4, 0, 0, 0, 0 }, 1003 {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0 }, 1004 {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0 }, 1005 {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0 }, 1006 {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4} 1007 }; 1008 const PetscReal bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}; 1009 const PetscReal binterpt[6][3] = { 1010 {6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025. }, 1011 {0, 0, 0 }, 1012 {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035. }, 1013 {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 1014 {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469. }, 1015 {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.} 1016 }; 1017 PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 1018 } 1019 { 1020 const PetscReal A[8][8] = { 1021 {0, 0, 0, 0, 0, 0, 0, 0}, 1022 {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 1023 {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 1024 {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 1025 {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 1026 {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 1027 {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 1028 {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 1029 }; 1030 const PetscReal At[8][8] = { 1031 {0, 0, 0, 0, 0, 0, 0, 0 }, 1032 {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0 }, 1033 {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0 }, 1034 {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0 }, 1035 {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0 }, 1036 {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0 }, 1037 {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0 }, 1038 {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.} 1039 }; 1040 const PetscReal bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}; 1041 const PetscReal binterpt[8][3] = { 1042 {-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439. }, 1043 {0, 0, 0 }, 1044 {0, 0, 0 }, 1045 {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, 1046 {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, 1047 {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, 1048 {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, 1049 {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460. } 1050 }; 1051 PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 1052 } 1053 #undef RC 1054 #undef us2 1055 #undef s2 1056 PetscFunctionReturn(PETSC_SUCCESS); 1057 } 1058 1059 /*@C 1060 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`. 1061 1062 Not Collective 1063 1064 Level: advanced 1065 1066 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 1067 @*/ 1068 PetscErrorCode TSARKIMEXRegisterDestroy(void) 1069 { 1070 ARKTableauLink link; 1071 1072 PetscFunctionBegin; 1073 while ((link = ARKTableauList)) { 1074 ARKTableau t = &link->tab; 1075 ARKTableauList = link->next; 1076 PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 1077 PetscCall(PetscFree2(t->bembedt, t->bembed)); 1078 PetscCall(PetscFree2(t->binterpt, t->binterp)); 1079 PetscCall(PetscFree(t->name)); 1080 PetscCall(PetscFree(link)); 1081 } 1082 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 1083 PetscFunctionReturn(PETSC_SUCCESS); 1084 } 1085 1086 /*@C 1087 TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called 1088 from `TSInitializePackage()`. 1089 1090 Level: developer 1091 1092 .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()` 1093 @*/ 1094 PetscErrorCode TSARKIMEXInitializePackage(void) 1095 { 1096 PetscFunctionBegin; 1097 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 1098 TSARKIMEXPackageInitialized = PETSC_TRUE; 1099 PetscCall(TSARKIMEXRegisterAll()); 1100 PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 1101 PetscFunctionReturn(PETSC_SUCCESS); 1102 } 1103 1104 /*@C 1105 TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is 1106 called from `PetscFinalize()`. 1107 1108 Level: developer 1109 1110 .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()` 1111 @*/ 1112 PetscErrorCode TSARKIMEXFinalizePackage(void) 1113 { 1114 PetscFunctionBegin; 1115 TSARKIMEXPackageInitialized = PETSC_FALSE; 1116 PetscCall(TSARKIMEXRegisterDestroy()); 1117 PetscFunctionReturn(PETSC_SUCCESS); 1118 } 1119 1120 /*@C 1121 TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 1122 1123 Logically Collective. 1124 1125 Input Parameters: 1126 + name - identifier for method 1127 . order - approximation order of method 1128 . s - number of stages, this is the dimension of the matrices below 1129 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 1130 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 1131 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 1132 . A - Non-stiff stage coefficients (dimension s*s, row-major) 1133 . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 1134 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 1135 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 1136 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 1137 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 1138 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 1139 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 1140 1141 Level: advanced 1142 1143 Note: 1144 Several `TSARKIMEX` methods are provided, this function is only needed to create new methods. 1145 1146 .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS` 1147 @*/ 1148 PetscErrorCode TSARKIMEXRegister(TSARKIMEXType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal A[], const PetscReal b[], const PetscReal c[], const PetscReal bembedt[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[], const PetscReal binterp[]) 1149 { 1150 ARKTableauLink link; 1151 ARKTableau t; 1152 PetscInt i, j; 1153 1154 PetscFunctionBegin; 1155 PetscCall(TSARKIMEXInitializePackage()); 1156 for (link = ARKTableauList; link; link = link->next) { 1157 PetscBool match; 1158 1159 PetscCall(PetscStrcmp(link->tab.name, name, &match)); 1160 PetscCheck(!match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Method with name \"%s\" already registered", name); 1161 } 1162 PetscCall(PetscNew(&link)); 1163 t = &link->tab; 1164 PetscCall(PetscStrallocpy(name, &t->name)); 1165 t->order = order; 1166 t->s = s; 1167 PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 1168 PetscCall(PetscArraycpy(t->At, At, s * s)); 1169 if (A) { 1170 PetscCall(PetscArraycpy(t->A, A, s * s)); 1171 t->additive = PETSC_TRUE; 1172 } 1173 1174 if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 1175 else 1176 for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 1177 1178 if (t->additive) { 1179 if (b) PetscCall(PetscArraycpy(t->b, b, s)); 1180 else 1181 for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 1182 } else PetscCall(PetscArrayzero(t->b, s)); 1183 1184 if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 1185 else 1186 for (i = 0; i < s; i++) 1187 for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 1188 1189 if (t->additive) { 1190 if (c) PetscCall(PetscArraycpy(t->c, c, s)); 1191 else 1192 for (i = 0; i < s; i++) 1193 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 1194 } else PetscCall(PetscArrayzero(t->c, s)); 1195 1196 t->stiffly_accurate = PETSC_TRUE; 1197 for (i = 0; i < s; i++) 1198 if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 1199 1200 t->explicit_first_stage = PETSC_TRUE; 1201 for (i = 0; i < s; i++) 1202 if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 1203 1204 /* def of FSAL can be made more precise */ 1205 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 1206 1207 if (bembedt) { 1208 PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 1209 PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 1210 PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 1211 } 1212 1213 t->pinterp = pinterp; 1214 PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 1215 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 1216 PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 1217 1218 link->next = ARKTableauList; 1219 ARKTableauList = link; 1220 PetscFunctionReturn(PETSC_SUCCESS); 1221 } 1222 1223 /*@C 1224 TSDIRKRegister - register a `TSDIRK` scheme by providing the entries in its Butcher tableau and, optionally, embedded approximations and interpolation 1225 1226 Logically Collective. 1227 1228 Input Parameters: 1229 + name - identifier for method 1230 . order - approximation order of method 1231 . s - number of stages, this is the dimension of the matrices below 1232 . At - Butcher table of stage coefficients (dimension `s`*`s`, row-major order) 1233 . bt - Butcher table for completing the step (dimension `s`; pass `NULL` to use the last row of `At`) 1234 . ct - Abscissa of each stage (dimension s, NULL to use row sums of At) 1235 . bembedt - Stiff part of completion table for embedded method (dimension s; `NULL` if not available) 1236 . pinterp - Order of the interpolation scheme, equal to the number of columns of `binterpt` and `binterp` 1237 - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 1238 1239 Level: advanced 1240 1241 Note: 1242 Several `TSDIRK` methods are provided, the use of this function is only needed to create new methods. 1243 1244 .seealso: [](ch_ts), `TSDIRK`, `TSType`, `TS` 1245 @*/ 1246 PetscErrorCode TSDIRKRegister(TSDIRKType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal bembedt[], PetscInt pinterp, const PetscReal binterpt[]) 1247 { 1248 PetscFunctionBegin; 1249 PetscCall(TSARKIMEXRegister(name, order, s, At, bt, ct, NULL, NULL, NULL, bembedt, NULL, pinterp, binterpt, NULL)); 1250 PetscFunctionReturn(PETSC_SUCCESS); 1251 } 1252 1253 /* 1254 The step completion formula is 1255 1256 x1 = x0 - h bt^T YdotI + h b^T YdotRHS 1257 1258 This function can be called before or after ts->vec_sol has been updated. 1259 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 1260 We can write 1261 1262 x1e = x0 - h bet^T YdotI + h be^T YdotRHS 1263 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 1264 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 1265 1266 so we can evaluate the method with different order even after the step has been optimistically completed. 1267 */ 1268 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) 1269 { 1270 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1271 ARKTableau tab = ark->tableau; 1272 PetscScalar *w = ark->work; 1273 PetscReal h; 1274 PetscInt s = tab->s, j; 1275 PetscBool hasE; 1276 1277 PetscFunctionBegin; 1278 switch (ark->status) { 1279 case TS_STEP_INCOMPLETE: 1280 case TS_STEP_PENDING: 1281 h = ts->time_step; 1282 break; 1283 case TS_STEP_COMPLETE: 1284 h = ts->ptime - ts->ptime_prev; 1285 break; 1286 default: 1287 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1288 } 1289 if (order == tab->order) { 1290 if (ark->status == TS_STEP_INCOMPLETE) { 1291 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 1292 PetscCall(VecCopy(ark->Y[s - 1], X)); 1293 } else { /* Use the standard completion formula (bt,b) */ 1294 PetscCall(VecCopy(ts->vec_sol, X)); 1295 for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 1296 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1297 if (tab->additive && ark->imex) { /* Method is IMEX, complete the explicit formula */ 1298 PetscCall(TSHasRHSFunction(ts, &hasE)); 1299 if (hasE) { 1300 for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 1301 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1302 } 1303 } 1304 } 1305 } else PetscCall(VecCopy(ts->vec_sol, X)); 1306 if (done) *done = PETSC_TRUE; 1307 PetscFunctionReturn(PETSC_SUCCESS); 1308 } else if (order == tab->order - 1) { 1309 if (!tab->bembedt) goto unavailable; 1310 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 1311 PetscCall(VecCopy(ts->vec_sol, X)); 1312 for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 1313 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 1314 if (tab->additive) { 1315 PetscCall(TSHasRHSFunction(ts, &hasE)); 1316 if (hasE) { 1317 for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 1318 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1319 } 1320 } 1321 } else { /* Rollback and re-complete using (bet-be,be-b) */ 1322 PetscCall(VecCopy(ts->vec_sol, X)); 1323 for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 1324 PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 1325 if (tab->additive) { 1326 PetscCall(TSHasRHSFunction(ts, &hasE)); 1327 if (hasE) { 1328 for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 1329 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 1330 } 1331 } 1332 } 1333 if (done) *done = PETSC_TRUE; 1334 PetscFunctionReturn(PETSC_SUCCESS); 1335 } 1336 unavailable: 1337 if (done) *done = PETSC_FALSE; 1338 else 1339 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "ARKIMEX '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.", tab->name, 1340 tab->order, order); 1341 PetscFunctionReturn(PETSC_SUCCESS); 1342 } 1343 1344 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 1345 { 1346 Vec Udot, Y1, Y2; 1347 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1348 PetscReal norm; 1349 1350 PetscFunctionBegin; 1351 PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 1352 PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 1353 PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 1354 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 1355 PetscCall(VecSetRandom(Udot, NULL)); 1356 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 1357 PetscCall(VecAXPY(Y2, -1.0, Y1)); 1358 PetscCall(VecAXPY(Y2, -1.0, Udot)); 1359 PetscCall(VecNorm(Y2, NORM_2, &norm)); 1360 if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 1361 *id = PETSC_TRUE; 1362 } else { 1363 *id = PETSC_FALSE; 1364 PetscCall(PetscInfo((PetscObject)ts, "IFunction(Udot = random) - IFunction(Udot = 0) is not near Udot, %g, suspect mass matrix implied in IFunction() is not the identity as required\n", (double)norm)); 1365 } 1366 PetscCall(VecDestroy(&Udot)); 1367 PetscCall(VecDestroy(&Y1)); 1368 PetscCall(VecDestroy(&Y2)); 1369 PetscFunctionReturn(PETSC_SUCCESS); 1370 } 1371 1372 static PetscErrorCode TSRollBack_ARKIMEX(TS ts) 1373 { 1374 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1375 ARKTableau tab = ark->tableau; 1376 const PetscInt s = tab->s; 1377 const PetscReal *bt = tab->bt, *b = tab->b; 1378 PetscScalar *w = ark->work; 1379 Vec *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS; 1380 PetscInt j; 1381 PetscReal h; 1382 1383 PetscFunctionBegin; 1384 switch (ark->status) { 1385 case TS_STEP_INCOMPLETE: 1386 case TS_STEP_PENDING: 1387 h = ts->time_step; 1388 break; 1389 case TS_STEP_COMPLETE: 1390 h = ts->ptime - ts->ptime_prev; 1391 break; 1392 default: 1393 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1394 } 1395 for (j = 0; j < s; j++) w[j] = -h * bt[j]; 1396 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotI)); 1397 if (tab->additive) { 1398 PetscBool hasE; 1399 1400 PetscCall(TSHasRHSFunction(ts, &hasE)); 1401 if (hasE) { 1402 for (j = 0; j < s; j++) w[j] = -h * b[j]; 1403 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS)); 1404 } 1405 } 1406 PetscFunctionReturn(PETSC_SUCCESS); 1407 } 1408 1409 static PetscErrorCode TSStep_ARKIMEX(TS ts) 1410 { 1411 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1412 ARKTableau tab = ark->tableau; 1413 const PetscInt s = tab->s; 1414 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 1415 PetscScalar *w = ark->work; 1416 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 1417 PetscBool extrapolate = ark->extrapolate; 1418 TSAdapt adapt; 1419 SNES snes; 1420 PetscInt i, j, its, lits; 1421 PetscInt rejections = 0; 1422 PetscBool hasE = PETSC_FALSE, dirk = (PetscBool)(!tab->additive), stageok, accept = PETSC_TRUE; 1423 PetscReal next_time_step = ts->time_step; 1424 1425 PetscFunctionBegin; 1426 if (ark->extrapolate && !ark->Y_prev) { 1427 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1428 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1429 if (tab->additive) PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1430 } 1431 1432 if (!dirk) PetscCall(TSHasRHSFunction(ts, &hasE)); 1433 if (!hasE) dirk = PETSC_TRUE; 1434 1435 if (!ts->steprollback) { 1436 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 1437 PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 1438 } 1439 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 1440 for (i = 0; i < s; i++) { 1441 PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 1442 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 1443 if (tab->additive && hasE) PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 1444 } 1445 } 1446 } 1447 1448 /* For fully implicit formulations we can solve the equations 1449 F(tn,xn,xdot) = 0 1450 for the explicit first stage */ 1451 if (dirk && tab->explicit_first_stage && ts->steprestart) { 1452 ark->scoeff = PETSC_MAX_REAL; 1453 PetscCall(VecCopy(ts->vec_sol, Z)); 1454 PetscCall(TSGetSNES(ts, &snes)); 1455 PetscCall(SNESSolve(snes, NULL, Ydot0)); 1456 } 1457 1458 /* For IMEX we compute a step */ 1459 if (!dirk && ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 1460 TS ts_start; 1461 if (PetscDefined(USE_DEBUG) && hasE) { 1462 PetscBool id = PETSC_FALSE; 1463 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1464 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "This scheme requires an identity mass matrix, however the TSIFunction you provided does not utilize an identity mass matrix"); 1465 } 1466 PetscCall(TSClone(ts, &ts_start)); 1467 PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 1468 PetscCall(TSSetTime(ts_start, ts->ptime)); 1469 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 1470 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 1471 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 1472 PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 1473 PetscCall(TSSetType(ts_start, TSARKIMEX)); 1474 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 1475 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 1476 1477 PetscCall(TSRestartStep(ts_start)); 1478 PetscCall(TSSolve(ts_start, ts->vec_sol)); 1479 PetscCall(TSGetTime(ts_start, &ts->ptime)); 1480 PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 1481 1482 { /* Save the initial slope for the next step */ 1483 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 1484 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 1485 } 1486 ts->steps++; 1487 1488 /* Set the correct TS in SNES */ 1489 /* We'll try to bypass this by changing the method on the fly */ 1490 { 1491 PetscCall(TSGetSNES(ts, &snes)); 1492 PetscCall(TSSetSNES(ts, snes)); 1493 } 1494 PetscCall(TSDestroy(&ts_start)); 1495 } 1496 1497 ark->status = TS_STEP_INCOMPLETE; 1498 while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 1499 PetscReal t = ts->ptime; 1500 PetscReal h = ts->time_step; 1501 for (i = 0; i < s; i++) { 1502 ark->stage_time = t + h * ct[i]; 1503 PetscCall(TSPreStage(ts, ark->stage_time)); 1504 if (At[i * s + i] == 0) { /* This stage is explicit */ 1505 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"); 1506 PetscCall(VecCopy(ts->vec_sol, Y[i])); 1507 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1508 PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 1509 if (tab->additive && hasE) { 1510 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1511 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 1512 } 1513 } else { 1514 ark->scoeff = 1. / At[i * s + i]; 1515 /* Ydot = shift*(Y-Z) */ 1516 PetscCall(VecCopy(ts->vec_sol, Z)); 1517 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 1518 PetscCall(VecMAXPY(Z, i, w, YdotI)); 1519 if (tab->additive && hasE) { 1520 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 1521 PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 1522 } 1523 if (extrapolate && !ts->steprestart) { 1524 /* Initial guess extrapolated from previous time step stage values */ 1525 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 1526 } else { 1527 /* Initial guess taken from last stage */ 1528 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 1529 } 1530 PetscCall(TSGetSNES(ts, &snes)); 1531 PetscCall(SNESSolve(snes, NULL, Y[i])); 1532 PetscCall(SNESGetIterationNumber(snes, &its)); 1533 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 1534 ts->snes_its += its; 1535 ts->ksp_its += lits; 1536 PetscCall(TSGetAdapt(ts, &adapt)); 1537 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 1538 if (!stageok) { 1539 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 1540 * use extrapolation to initialize the solves on the next attempt. */ 1541 extrapolate = PETSC_FALSE; 1542 goto reject_step; 1543 } 1544 } 1545 if (dirk || ts->equation_type >= TS_EQ_IMPLICIT) { 1546 if (i == 0 && tab->explicit_first_stage) { 1547 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", 1548 ((PetscObject)ts)->type_name, ark->tableau->name); 1549 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 1550 } else { 1551 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1552 } 1553 } else { 1554 if (i == 0 && tab->explicit_first_stage) { 1555 PetscCall(VecZeroEntries(Ydot)); 1556 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 1557 PetscCall(VecScale(YdotI[i], -1.0)); 1558 } else { 1559 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 1560 } 1561 if (hasE) { 1562 if (ark->imex) { 1563 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 1564 } else { 1565 PetscCall(VecZeroEntries(YdotRHS[i])); 1566 } 1567 } 1568 } 1569 PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 1570 } 1571 1572 ark->status = TS_STEP_INCOMPLETE; 1573 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 1574 ark->status = TS_STEP_PENDING; 1575 PetscCall(TSGetAdapt(ts, &adapt)); 1576 PetscCall(TSAdaptCandidatesClear(adapt)); 1577 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 1578 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1579 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1580 if (!accept) { /* Roll back the current step */ 1581 PetscCall(TSRollBack_ARKIMEX(ts)); 1582 ts->time_step = next_time_step; 1583 goto reject_step; 1584 } 1585 1586 ts->ptime += ts->time_step; 1587 ts->time_step = next_time_step; 1588 break; 1589 1590 reject_step: 1591 ts->reject++; 1592 accept = PETSC_FALSE; 1593 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1594 ts->reason = TS_DIVERGED_STEP_REJECTED; 1595 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1596 } 1597 } 1598 PetscFunctionReturn(PETSC_SUCCESS); 1599 } 1600 1601 /* 1602 This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 1603 Udot = H(t,U) + G(t,U) 1604 This corresponds to F(t,U,Udot) = Udot-H(t,U). 1605 1606 The complete adjoint equations are 1607 (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 1608 dGdU (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 1609 + dHdU (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 1610 lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 1611 mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 1612 + dGdP (b_i lambda_{n+1} + \sum_{j=i+1}^s a[j][i] lambda_s[j]) 1613 + dHdP (bt[i] lambda_{n+1} + \sum_{j=i+1}^s at[j][i] lambda_s[j])), i = s-1,...,0 1614 mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 1615 where shift = 1/(at[i][i]*h) 1616 1617 If at[i][i] is 0, the first equation falls back to 1618 lambda_s[i] = h ( 1619 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 1620 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 1621 1622 */ 1623 static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 1624 { 1625 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1626 ARKTableau tab = ark->tableau; 1627 const PetscInt s = tab->s; 1628 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 1629 PetscScalar *w = ark->work; 1630 Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 1631 Mat Jex, Jim, Jimpre; 1632 PetscInt i, j, nadj; 1633 PetscReal t = ts->ptime, stage_time_ex; 1634 PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 1635 KSP ksp; 1636 1637 PetscFunctionBegin; 1638 ark->status = TS_STEP_INCOMPLETE; 1639 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1640 PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 1641 PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 1642 1643 for (i = s - 1; i >= 0; i--) { 1644 ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 1645 stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 1646 if (At[i * s + i] == 0) { // This stage is explicit 1647 ark->scoeff = 0.; 1648 } else { 1649 ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 1650 } 1651 PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 1652 PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 1653 if (ts->vecs_sensip) { 1654 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 1655 PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 1656 } 1657 /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 1658 for (nadj = 0; nadj < ts->numcost; nadj++) { 1659 /* build implicit part */ 1660 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1661 if (s - i - 1 > 0) { 1662 /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 1663 for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 1664 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1665 } 1666 /* Temp = Temp - bt[i] lambda_{n+1} */ 1667 PetscCall(VecAXPY(VecsSensiTemp[nadj], -bt[i], ts->vecs_sensi[nadj])); 1668 if (bt[i] || s - i - 1 > 0) { 1669 /* (shift I - dHdU) Temp */ 1670 PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 1671 /* cancel out shift Temp where shift=-scoeff/h */ 1672 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 1673 if (ts->vecs_sensip) { 1674 /* - dHdP Temp */ 1675 PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1676 /* mu_n += -h dHdP Temp */ 1677 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1678 } 1679 } else { 1680 PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 1681 } 1682 /* build explicit part */ 1683 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 1684 if (s - i - 1 > 0) { 1685 /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 1686 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 1687 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 1688 } 1689 /* Temp = Temp + b[i] lambda_{n+1} */ 1690 PetscCall(VecAXPY(VecsSensiTemp[nadj], b[i], ts->vecs_sensi[nadj])); 1691 if (b[i] || s - i - 1 > 0) { 1692 /* dGdU Temp */ 1693 PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1694 if (ts->vecs_sensip) { 1695 /* dGdP Temp */ 1696 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 1697 /* mu_n += h dGdP Temp */ 1698 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 1699 } 1700 } 1701 /* Build LHS for first-order adjoint */ 1702 if (At[i * s + i] == 0) { // This stage is explicit 1703 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 1704 } else { 1705 KSP ksp; 1706 KSPConvergedReason kspreason; 1707 PetscCall(SNESGetKSP(ts->snes, &ksp)); 1708 PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 1709 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 1710 PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 1711 PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 1712 if (kspreason < 0) { 1713 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 1714 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj)); 1715 } 1716 if (ts->vecs_sensip) { 1717 /* -dHdP lambda_s[i] */ 1718 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 1719 /* mu_n += h at[i][i] dHdP lambda_s[i] */ 1720 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 1721 } 1722 } 1723 } 1724 } 1725 for (j = 0; j < s; j++) w[j] = 1.0; 1726 for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 1727 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 1728 ark->status = TS_STEP_COMPLETE; 1729 PetscFunctionReturn(PETSC_SUCCESS); 1730 } 1731 1732 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1733 { 1734 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1735 ARKTableau tab = ark->tableau; 1736 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1737 PetscReal h; 1738 PetscReal tt, t; 1739 PetscScalar *bt = ark->work, *b = ark->work + s; 1740 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1741 1742 PetscFunctionBegin; 1743 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "%s %s does not have an interpolation formula", ((PetscObject)ts)->type_name, ark->tableau->name); 1744 switch (ark->status) { 1745 case TS_STEP_INCOMPLETE: 1746 case TS_STEP_PENDING: 1747 h = ts->time_step; 1748 t = (itime - ts->ptime) / h; 1749 break; 1750 case TS_STEP_COMPLETE: 1751 h = ts->ptime - ts->ptime_prev; 1752 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1753 break; 1754 default: 1755 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1756 } 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 PetscCall(VecCopy(ark->Y[0], X)); 1765 PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1766 if (tab->additive) { 1767 PetscBool hasE; 1768 PetscCall(TSHasRHSFunction(ts, &hasE)); 1769 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1770 } 1771 PetscFunctionReturn(PETSC_SUCCESS); 1772 } 1773 1774 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1775 { 1776 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1777 ARKTableau tab = ark->tableau; 1778 PetscInt s = tab->s, pinterp = tab->pinterp, i, j; 1779 PetscReal h, h_prev, t, tt; 1780 PetscScalar *bt = ark->work, *b = ark->work + s; 1781 const PetscReal *Bt = tab->binterpt, *B = tab->binterp; 1782 1783 PetscFunctionBegin; 1784 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 1785 h = ts->time_step; 1786 h_prev = ts->ptime - ts->ptime_prev; 1787 t = 1 + h / h_prev * c; 1788 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1789 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1790 for (i = 0; i < s; i++) { 1791 bt[i] += h * Bt[i * pinterp + j] * tt; 1792 b[i] += h * B[i * pinterp + j] * tt; 1793 } 1794 } 1795 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 1796 PetscCall(VecCopy(ark->Y_prev[0], X)); 1797 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1798 if (tab->additive) { 1799 PetscBool hasE; 1800 PetscCall(TSHasRHSFunction(ts, &hasE)); 1801 if (hasE) PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1802 } 1803 PetscFunctionReturn(PETSC_SUCCESS); 1804 } 1805 1806 static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1807 { 1808 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1809 ARKTableau tab = ark->tableau; 1810 1811 PetscFunctionBegin; 1812 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1813 PetscCall(PetscFree(ark->work)); 1814 PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 1815 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 1816 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 1817 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 1818 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 1819 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 1820 PetscFunctionReturn(PETSC_SUCCESS); 1821 } 1822 1823 static PetscErrorCode TSReset_ARKIMEX(TS ts) 1824 { 1825 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1826 1827 PetscFunctionBegin; 1828 PetscCall(TSARKIMEXTableauReset(ts)); 1829 PetscCall(VecDestroy(&ark->Ydot)); 1830 PetscCall(VecDestroy(&ark->Ydot0)); 1831 PetscCall(VecDestroy(&ark->Z)); 1832 PetscFunctionReturn(PETSC_SUCCESS); 1833 } 1834 1835 static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 1836 { 1837 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1838 ARKTableau tab = ark->tableau; 1839 1840 PetscFunctionBegin; 1841 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 1842 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 1843 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 1844 PetscFunctionReturn(PETSC_SUCCESS); 1845 } 1846 1847 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1848 { 1849 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1850 1851 PetscFunctionBegin; 1852 if (Z) { 1853 if (dm && dm != ts->dm) { 1854 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1855 } else *Z = ax->Z; 1856 } 1857 if (Ydot) { 1858 if (dm && dm != ts->dm) { 1859 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1860 } else *Ydot = ax->Ydot; 1861 } 1862 PetscFunctionReturn(PETSC_SUCCESS); 1863 } 1864 1865 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1866 { 1867 PetscFunctionBegin; 1868 if (Z) { 1869 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1870 } 1871 if (Ydot) { 1872 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1873 } 1874 PetscFunctionReturn(PETSC_SUCCESS); 1875 } 1876 1877 /* 1878 This defines the nonlinear equation that is to be solved with SNES 1879 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1880 */ 1881 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1882 { 1883 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1884 DM dm, dmsave; 1885 Vec Z, Ydot; 1886 1887 PetscFunctionBegin; 1888 PetscCall(SNESGetDM(snes, &dm)); 1889 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1890 dmsave = ts->dm; 1891 ts->dm = dm; 1892 1893 if (ark->scoeff == PETSC_MAX_REAL) { 1894 /* We are solving F(t,x_n,xdot) = 0 to start the method */ 1895 PetscCall(TSComputeIFunction(ts, ark->stage_time, Z, X, F, ark->imex)); 1896 } else { 1897 PetscReal shift = ark->scoeff / ts->time_step; 1898 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1899 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1900 } 1901 1902 ts->dm = dmsave; 1903 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1904 PetscFunctionReturn(PETSC_SUCCESS); 1905 } 1906 1907 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1908 { 1909 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1910 DM dm, dmsave; 1911 Vec Ydot, Z; 1912 PetscReal shift; 1913 1914 PetscFunctionBegin; 1915 PetscCall(SNESGetDM(snes, &dm)); 1916 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1917 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1918 dmsave = ts->dm; 1919 ts->dm = dm; 1920 1921 if (ark->scoeff == PETSC_MAX_REAL) { 1922 /* We are solving F(t,x_n,xdot) = 0 to start the method, we only only dF/dXdot 1923 Jed's proposal is to compute with a very large shift and scale back the matrix */ 1924 shift = 1.0 / PETSC_MACHINE_EPSILON; 1925 PetscCall(TSComputeIJacobian(ts, ark->stage_time, Z, X, shift, A, B, ark->imex)); 1926 PetscCall(MatScale(B, PETSC_MACHINE_EPSILON)); 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 TSIFunction 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->rollback = TSRollBack_ARKIMEX; 2350 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 2351 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 2352 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 2353 ts->ops->getstages = TSGetStages_ARKIMEX; 2354 ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 2355 2356 ts->usessnes = PETSC_TRUE; 2357 2358 PetscCall(PetscNew(&ark)); 2359 ts->data = (void *)ark; 2360 ark->imex = dirk ? PETSC_FALSE : PETSC_TRUE; 2361 2362 ark->VecsDeltaLam = NULL; 2363 ark->VecsSensiTemp = NULL; 2364 ark->VecsSensiPTemp = NULL; 2365 2366 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 2367 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 2368 if (!dirk) { 2369 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 2370 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 2371 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 2372 } 2373 PetscFunctionReturn(PETSC_SUCCESS); 2374 } 2375 2376 /* ------------------------------------------------------------ */ 2377 2378 static PetscErrorCode TSDIRKSetType_DIRK(TS ts, TSDIRKType dirktype) 2379 { 2380 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 2381 2382 PetscFunctionBegin; 2383 PetscCall(TSARKIMEXSetType_ARKIMEX(ts, dirktype)); 2384 PetscCheck(!ark->tableau->additive, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Method \"%s\" is not DIRK", dirktype); 2385 PetscFunctionReturn(PETSC_SUCCESS); 2386 } 2387 2388 /*@C 2389 TSDIRKSetType - Set the type of `TSDIRK` scheme 2390 2391 Logically Collective 2392 2393 Input Parameters: 2394 + ts - timestepping context 2395 - dirktype - type of `TSDIRK` scheme 2396 2397 Options Database Key: 2398 . -ts_dirkimex_type - set `TSDIRK` scheme type 2399 2400 Level: intermediate 2401 2402 .seealso: [](ch_ts), `TSDIRKGetType()`, `TSDIRK`, `TSDIRKType` 2403 @*/ 2404 PetscErrorCode TSDIRKSetType(TS ts, TSDIRKType dirktype) 2405 { 2406 PetscFunctionBegin; 2407 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2408 PetscAssertPointer(dirktype, 2); 2409 PetscTryMethod(ts, "TSDIRKSetType_C", (TS, TSDIRKType), (ts, dirktype)); 2410 PetscFunctionReturn(PETSC_SUCCESS); 2411 } 2412 2413 /*@C 2414 TSDIRKGetType - Get the type of `TSDIRK` scheme 2415 2416 Logically Collective 2417 2418 Input Parameter: 2419 . ts - timestepping context 2420 2421 Output Parameter: 2422 . dirktype - type of `TSDIRK` scheme 2423 2424 Level: intermediate 2425 2426 .seealso: [](ch_ts), `TSDIRKSetType()` 2427 @*/ 2428 PetscErrorCode TSDIRKGetType(TS ts, TSDIRKType *dirktype) 2429 { 2430 PetscFunctionBegin; 2431 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2432 PetscUseMethod(ts, "TSDIRKGetType_C", (TS, TSDIRKType *), (ts, dirktype)); 2433 PetscFunctionReturn(PETSC_SUCCESS); 2434 } 2435 2436 /*MC 2437 TSDIRK - ODE and DAE solver using Diagonally implicit Runge-Kutta schemes. 2438 2439 Level: beginner 2440 2441 Notes: 2442 The default is `TSDIRKES213SAL`, it can be changed with `TSDIRKSetType()` or -ts_dirk_type. 2443 The convention used in PETSc to name the DIRK methods is TSDIRK[E][S]PQS[SA][L][A] with: 2444 + E - whether the method has an explicit first stage 2445 . S - whether the method is single diagonal 2446 . P - order of the advancing method 2447 . Q - order of the embedded method 2448 . S - number of stages 2449 . SA - whether the method is stiffly accurate 2450 . L - whether the method is L-stable 2451 - A - whether the method is A-stable 2452 2453 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSDIRKSetType()`, `TSDIRKGetType()`, `TSDIRKRegister()`. 2454 M*/ 2455 PETSC_EXTERN PetscErrorCode TSCreate_DIRK(TS ts) 2456 { 2457 PetscFunctionBegin; 2458 PetscCall(TSCreate_ARKIMEX(ts)); 2459 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKGetType_C", TSARKIMEXGetType_ARKIMEX)); 2460 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSDIRKSetType_C", TSDIRKSetType_DIRK)); 2461 PetscCall(TSDIRKSetType(ts, TSDIRKDefault)); 2462 PetscFunctionReturn(PETSC_SUCCESS); 2463 } 2464