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