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