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