1 /* 2 Code for timestepping with additive Runge-Kutta IMEX method 3 4 Notes: 5 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 PetscBool TSARKIMEXRegisterAllCalled; 17 static PetscBool TSARKIMEXPackageInitialized; 18 static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec); 19 20 typedef struct _ARKTableau *ARKTableau; 21 struct _ARKTableau { 22 char *name; 23 PetscInt order; /* Classical approximation order of the method */ 24 PetscInt s; /* Number of stages */ 25 PetscBool stiffly_accurate; /* The implicit part is stiffly accurate*/ 26 PetscBool FSAL_implicit; /* The implicit part is FSAL*/ 27 PetscBool explicit_first_stage; /* The implicit part has an explicit first stage*/ 28 PetscInt pinterp; /* Interpolation order */ 29 PetscReal *At, *bt, *ct; /* Stiff tableau */ 30 PetscReal *A, *b, *c; /* Non-stiff tableau */ 31 PetscReal *bembedt, *bembed; /* Embedded formula of order one less (order-1) */ 32 PetscReal *binterpt, *binterp; /* Dense output formula */ 33 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 34 }; 35 typedef struct _ARKTableauLink *ARKTableauLink; 36 struct _ARKTableauLink { 37 struct _ARKTableau tab; 38 ARKTableauLink next; 39 }; 40 static ARKTableauLink ARKTableauList; 41 42 typedef struct { 43 ARKTableau tableau; 44 Vec *Y; /* States computed during the step */ 45 Vec *YdotI; /* Time derivatives for the stiff part */ 46 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 47 Vec *Y_prev; /* States computed during the previous time step */ 48 Vec *YdotI_prev; /* Time derivatives for the stiff part for the previous time step*/ 49 Vec *YdotRHS_prev; /* Function evaluations for the non-stiff part for the previous time step*/ 50 Vec Ydot0; /* Holds the slope from the previous step in FSAL case */ 51 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 52 Vec Z; /* Ydot = shift(Y-Z) */ 53 PetscScalar *work; /* Scalar work */ 54 PetscReal scoeff; /* shift = scoeff/dt */ 55 PetscReal stage_time; 56 PetscBool imex; 57 PetscBool extrapolate; /* Extrapolate initial guess from previous time-step stage values */ 58 TSStepStatus status; 59 60 /* context for sensitivity analysis */ 61 Vec *VecsDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */ 62 Vec *VecsSensiTemp; /* Vectors to be multiplied with Jacobian transpose */ 63 Vec *VecsSensiPTemp; /* Temporary Vectors to store JacobianP-transpose-vector product */ 64 } TS_ARKIMEX; 65 /*MC 66 TSARKIMEXARS122 - Second order ARK IMEX scheme. 67 68 This method has one explicit stage and one implicit stage. 69 70 Options Database Key: 71 . -ts_arkimex_type ars122 - set arkimex type to ars122 72 73 Level: advanced 74 75 References: 76 . * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 77 78 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 79 M*/ 80 /*MC 81 TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part. 82 83 This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu. 84 85 Options Database Key: 86 . -ts_arkimex_type a2 - set arkimex type to a2 87 88 Level: advanced 89 90 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 91 M*/ 92 /*MC 93 TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part. 94 95 This method has two implicit stages, and L-stable implicit scheme. 96 97 Options Database Key: 98 . -ts_arkimex_type l2 - set arkimex type to l2 99 100 Level: advanced 101 102 References: 103 . * - L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005. 104 105 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 106 M*/ 107 /*MC 108 TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method. 109 110 This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used. 111 112 Options Database Key: 113 . -ts_arkimex_type 1bee - set arkimex type to 1bee 114 115 Level: advanced 116 117 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 118 M*/ 119 /*MC 120 TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part. 121 122 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. 123 124 Options Database Key: 125 . -ts_arkimex_type 2c - set arkimex type to 2c 126 127 Level: advanced 128 129 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 130 M*/ 131 /*MC 132 TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 133 134 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. 135 136 Options Database Key: 137 . -ts_arkimex_type 2d - set arkimex type to 2d 138 139 Level: advanced 140 141 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 142 M*/ 143 /*MC 144 TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 145 146 This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu. 147 148 Options Database Key: 149 . -ts_arkimex_type 2e - set arkimex type to 2e 150 151 Level: advanced 152 153 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 154 M*/ 155 /*MC 156 TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme. 157 158 This method has three implicit stages. 159 160 References: 161 . * - L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005. 162 163 This method is referred to as SSP2-(3,3,2) in https://arxiv.org/abs/1110.4375 164 165 Options Database Key: 166 . -ts_arkimex_type prssp2 - set arkimex type to prssp2 167 168 Level: advanced 169 170 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 171 M*/ 172 /*MC 173 TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part. 174 175 This method has one explicit stage and three implicit stages. 176 177 Options Database Key: 178 . -ts_arkimex_type 3 - set arkimex type to 3 179 180 Level: advanced 181 182 References: 183 . * - Kennedy and Carpenter 2003. 184 185 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 186 M*/ 187 /*MC 188 TSARKIMEXARS443 - Third order ARK IMEX scheme. 189 190 This method has one explicit stage and four implicit stages. 191 192 Options Database Key: 193 . -ts_arkimex_type ars443 - set arkimex type to ars443 194 195 Level: advanced 196 197 References: 198 + * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 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 /*MC 204 TSARKIMEXBPR3 - Third order ARK IMEX scheme. 205 206 This method has one explicit stage and four implicit stages. 207 208 Options Database Key: 209 . -ts_arkimex_type bpr3 - set arkimex type to bpr3 210 211 Level: advanced 212 213 References: 214 . * - This method is referred to as ARK3 in https://arxiv.org/abs/1110.4375 215 216 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 217 M*/ 218 /*MC 219 TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part. 220 221 This method has one explicit stage and four implicit stages. 222 223 Options Database Key: 224 . -ts_arkimex_type 4 - set arkimex type to4 225 226 Level: advanced 227 228 References: 229 . * - Kennedy and Carpenter 2003. 230 231 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 232 M*/ 233 /*MC 234 TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part. 235 236 This method has one explicit stage and five implicit stages. 237 238 Options Database Key: 239 . -ts_arkimex_type 5 - set arkimex type to 5 240 241 Level: advanced 242 243 References: 244 . * - Kennedy and Carpenter 2003. 245 246 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 247 M*/ 248 249 /*@C 250 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX` 251 252 Not Collective, but should be called by all processes which will need the schemes to be registered 253 254 Level: advanced 255 256 .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()` 257 @*/ 258 PetscErrorCode TSARKIMEXRegisterAll(void) 259 { 260 PetscFunctionBegin; 261 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 262 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 263 264 { 265 const PetscReal A[3][3] = 266 { 267 {0.0, 0.0, 0.0}, 268 {0.0, 0.0, 0.0}, 269 {0.0, 0.5, 0.0} 270 }, 271 At[3][3] = {{1.0, 0.0, 0.0}, {0.0, 0.5, 0.0}, {0.0, 0.5, 0.5}}, b[3] = {0.0, 0.5, 0.5}, bembedt[3] = {1.0, 0.0, 0.0}; 272 PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 273 } 274 { 275 const PetscReal A[2][2] = 276 { 277 {0.0, 0.0}, 278 {0.5, 0.0} 279 }, 280 At[2][2] = {{0.0, 0.0}, {0.0, 0.5}}, b[2] = {0.0, 1.0}, bembedt[2] = {0.5, 0.5}; 281 /* 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*/ 282 PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 283 } 284 { 285 const PetscReal A[2][2] = 286 { 287 {0.0, 0.0}, 288 {1.0, 0.0} 289 }, 290 At[2][2] = {{0.0, 0.0}, {0.5, 0.5}}, b[2] = {0.5, 0.5}, bembedt[2] = {0.0, 1.0}; 291 /* 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*/ 292 PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 293 } 294 { 295 /* const PetscReal us2 = 1.0-1.0/PetscSqrtReal((PetscReal)2.0); Direct evaluation: 0.2928932188134524755992. Used below to ensure all values are available at compile time */ 296 const PetscReal A[2][2] = 297 { 298 {0.0, 0.0}, 299 {1.0, 0.0} 300 }, 301 At[2][2] = {{0.2928932188134524755992, 0.0}, {1.0 - 2.0 * 0.2928932188134524755992, 0.2928932188134524755992}}, b[2] = {0.5, 0.5}, bembedt[2] = {0.0, 1.0}, binterpt[2][2] = {{(0.2928932188134524755992 - 1.0) / (2.0 * 0.2928932188134524755992 - 1.0), -1 / (2.0 * (1.0 - 2.0 * 0.2928932188134524755992))}, {1 - (0.2928932188134524755992 - 1.0) / (2.0 * 0.2928932188134524755992 - 1.0), -1 / (2.0 * (1.0 - 2.0 * 0.2928932188134524755992))}}, binterp[2][2] = {{1.0, -0.5}, {0.0, 0.5}}; 302 PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 303 } 304 { 305 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */ 306 const PetscReal A[3][3] = 307 { 308 {0, 0, 0}, 309 {2 - 1.414213562373095048802, 0, 0}, 310 {0.5, 0.5, 0} 311 }, 312 At[3][3] = {{0, 0, 0}, {1 - 1 / 1.414213562373095048802, 1 - 1 / 1.414213562373095048802, 0}, {1 / (2 * 1.414213562373095048802), 1 / (2 * 1.414213562373095048802), 1 - 1 / 1.414213562373095048802}}, bembedt[3] = {(4. - 1.414213562373095048802) / 8., (4. - 1.414213562373095048802) / 8., 1 / (2. * 1.414213562373095048802)}, binterpt[3][2] = {{1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 - 1.414213562373095048802, 1.0 / 1.414213562373095048802}}; 313 PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 314 } 315 { 316 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */ 317 const PetscReal A[3][3] = 318 { 319 {0, 0, 0}, 320 {2 - 1.414213562373095048802, 0, 0}, 321 {0.75, 0.25, 0} 322 }, 323 At[3][3] = {{0, 0, 0}, {1 - 1 / 1.414213562373095048802, 1 - 1 / 1.414213562373095048802, 0}, {1 / (2 * 1.414213562373095048802), 1 / (2 * 1.414213562373095048802), 1 - 1 / 1.414213562373095048802}}, bembedt[3] = {(4. - 1.414213562373095048802) / 8., (4. - 1.414213562373095048802) / 8., 1 / (2. * 1.414213562373095048802)}, binterpt[3][2] = {{1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 - 1.414213562373095048802, 1.0 / 1.414213562373095048802}}; 324 PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 325 } 326 { /* Optimal for linear implicit part */ 327 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */ 328 const PetscReal A[3][3] = 329 { 330 {0, 0, 0}, 331 {2 - 1.414213562373095048802, 0, 0}, 332 {(3 - 2 * 1.414213562373095048802) / 6, (3 + 2 * 1.414213562373095048802) / 6, 0} 333 }, 334 At[3][3] = {{0, 0, 0}, {1 - 1 / 1.414213562373095048802, 1 - 1 / 1.414213562373095048802, 0}, {1 / (2 * 1.414213562373095048802), 1 / (2 * 1.414213562373095048802), 1 - 1 / 1.414213562373095048802}}, bembedt[3] = {(4. - 1.414213562373095048802) / 8., (4. - 1.414213562373095048802) / 8., 1 / (2. * 1.414213562373095048802)}, binterpt[3][2] = {{1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 - 1.414213562373095048802, 1.0 / 1.414213562373095048802}}; 335 PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 336 } 337 { /* Optimal for linear implicit part */ 338 const PetscReal A[3][3] = 339 { 340 {0, 0, 0}, 341 {0.5, 0, 0}, 342 {0.5, 0.5, 0} 343 }, 344 At[3][3] = {{0.25, 0, 0}, {0, 0.25, 0}, {1. / 3, 1. / 3, 1. / 3}}; 345 PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 346 } 347 { 348 const PetscReal A[4][4] = 349 { 350 {0, 0, 0, 0}, 351 {1767732205903. / 2027836641118., 0, 0, 0}, 352 {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 353 {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 354 }, 355 At[4][4] = {{0, 0, 0, 0}, {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0}, {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0}, {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.}}, bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}, binterpt[4][2] = {{4655552711362. / 22874653954995., -215264564351. / 13552729205753.}, {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119.}, {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, {584795268549. / 6622622206610., 2508943948391. / 7218656332882.}}; 356 PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 357 } 358 { 359 const PetscReal A[5][5] = 360 { 361 {0, 0, 0, 0, 0}, 362 {1. / 2, 0, 0, 0, 0}, 363 {11. / 18, 1. / 18, 0, 0, 0}, 364 {5. / 6, -5. / 6, .5, 0, 0}, 365 {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 366 }, 367 At[5][5] = {{0, 0, 0, 0, 0}, {0, 1. / 2, 0, 0, 0}, {0, 1. / 6, 1. / 2, 0, 0}, {0, -1. / 2, 1. / 2, 1. / 2, 0}, {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2}}, *bembedt = NULL; 368 PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 0, NULL, NULL)); 369 } 370 { 371 const PetscReal A[5][5] = 372 { 373 {0, 0, 0, 0, 0}, 374 {1, 0, 0, 0, 0}, 375 {4. / 9, 2. / 9, 0, 0, 0}, 376 {1. / 4, 0, 3. / 4, 0, 0}, 377 {1. / 4, 0, 3. / 5, 0, 0} 378 }, 379 At[5][5] = {{0, 0, 0, 0, 0}, {.5, .5, 0, 0, 0}, {5. / 18, -1. / 9, .5, 0, 0}, {.5, 0, 0, .5, 0}, {.25, 0, .75, -.5, .5}}, *bembedt = NULL; 380 PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 0, NULL, NULL)); 381 } 382 { 383 const PetscReal A[6][6] = 384 { 385 {0, 0, 0, 0, 0, 0}, 386 {1. / 2, 0, 0, 0, 0, 0}, 387 {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 388 {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 389 {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 390 {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 391 }, 392 At[6][6] = {{0, 0, 0, 0, 0, 0}, {1. / 4, 1. / 4, 0, 0, 0, 0}, {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0}, {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0}, {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0}, {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4}}, bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}, binterpt[6][3] = {{6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025.}, {0, 0, 0}, 393 {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035.}, {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 394 {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469.}, {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.}}; 395 PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 396 } 397 { 398 const PetscReal A[8][8] = 399 { 400 {0, 0, 0, 0, 0, 0, 0, 0}, 401 {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 402 {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 403 {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 404 {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 405 {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 406 {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 407 {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 408 }, 409 At[8][8] = {{0, 0, 0, 0, 0, 0, 0, 0}, {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0}, {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0}, {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0}, {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0}, {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0}, {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0}, {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.}}, bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}, binterpt[8][3] = {{-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439.}, {0, 0, 0}, {0, 0, 0}, {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460.}}; 410 PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 411 } 412 PetscFunctionReturn(PETSC_SUCCESS); 413 } 414 415 /*@C 416 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`. 417 418 Not Collective 419 420 Level: advanced 421 422 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 423 @*/ 424 PetscErrorCode TSARKIMEXRegisterDestroy(void) 425 { 426 ARKTableauLink link; 427 428 PetscFunctionBegin; 429 while ((link = ARKTableauList)) { 430 ARKTableau t = &link->tab; 431 ARKTableauList = link->next; 432 PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 433 PetscCall(PetscFree2(t->bembedt, t->bembed)); 434 PetscCall(PetscFree2(t->binterpt, t->binterp)); 435 PetscCall(PetscFree(t->name)); 436 PetscCall(PetscFree(link)); 437 } 438 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 439 PetscFunctionReturn(PETSC_SUCCESS); 440 } 441 442 /*@C 443 TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called 444 from `TSInitializePackage()`. 445 446 Level: developer 447 448 .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()` 449 @*/ 450 PetscErrorCode TSARKIMEXInitializePackage(void) 451 { 452 PetscFunctionBegin; 453 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 454 TSARKIMEXPackageInitialized = PETSC_TRUE; 455 PetscCall(TSARKIMEXRegisterAll()); 456 PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 457 PetscFunctionReturn(PETSC_SUCCESS); 458 } 459 460 /*@C 461 TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is 462 called from `PetscFinalize()`. 463 464 Level: developer 465 466 .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()` 467 @*/ 468 PetscErrorCode TSARKIMEXFinalizePackage(void) 469 { 470 PetscFunctionBegin; 471 TSARKIMEXPackageInitialized = PETSC_FALSE; 472 PetscCall(TSARKIMEXRegisterDestroy()); 473 PetscFunctionReturn(PETSC_SUCCESS); 474 } 475 476 /*@C 477 TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 478 479 Not Collective, but the same schemes should be registered on all processes on which they will be used 480 481 Input Parameters: 482 + name - identifier for method 483 . order - approximation order of method 484 . s - number of stages, this is the dimension of the matrices below 485 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 486 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 487 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 488 . A - Non-stiff stage coefficients (dimension s*s, row-major) 489 . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 490 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 491 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 492 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 493 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 494 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 495 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 496 497 Level: advanced 498 499 Note: 500 Several `TSARKIMEX` methods are provided, this function is only needed to create new methods. 501 502 .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS` 503 @*/ 504 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[]) 505 { 506 ARKTableauLink link; 507 ARKTableau t; 508 PetscInt i, j; 509 510 PetscFunctionBegin; 511 PetscCall(TSARKIMEXInitializePackage()); 512 PetscCall(PetscNew(&link)); 513 t = &link->tab; 514 PetscCall(PetscStrallocpy(name, &t->name)); 515 t->order = order; 516 t->s = s; 517 PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 518 PetscCall(PetscArraycpy(t->At, At, s * s)); 519 PetscCall(PetscArraycpy(t->A, A, s * s)); 520 if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 521 else 522 for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 523 if (b) PetscCall(PetscArraycpy(t->b, b, s)); 524 else 525 for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 526 if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 527 else 528 for (i = 0; i < s; i++) 529 for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 530 if (c) PetscCall(PetscArraycpy(t->c, c, s)); 531 else 532 for (i = 0; i < s; i++) 533 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 534 t->stiffly_accurate = PETSC_TRUE; 535 for (i = 0; i < s; i++) 536 if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 537 t->explicit_first_stage = PETSC_TRUE; 538 for (i = 0; i < s; i++) 539 if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 540 /*def of FSAL can be made more precise*/ 541 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 542 if (bembedt) { 543 PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 544 PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 545 PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 546 } 547 548 t->pinterp = pinterp; 549 PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 550 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 551 PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 552 link->next = ARKTableauList; 553 ARKTableauList = link; 554 PetscFunctionReturn(PETSC_SUCCESS); 555 } 556 557 /* 558 The step completion formula is 559 560 x1 = x0 - h bt^T YdotI + h b^T YdotRHS 561 562 This function can be called before or after ts->vec_sol has been updated. 563 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 564 We can write 565 566 x1e = x0 - h bet^T YdotI + h be^T YdotRHS 567 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 568 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 569 570 so we can evaluate the method with different order even after the step has been optimistically completed. 571 */ 572 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) 573 { 574 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 575 ARKTableau tab = ark->tableau; 576 PetscScalar *w = ark->work; 577 PetscReal h; 578 PetscInt s = tab->s, j; 579 580 PetscFunctionBegin; 581 switch (ark->status) { 582 case TS_STEP_INCOMPLETE: 583 case TS_STEP_PENDING: 584 h = ts->time_step; 585 break; 586 case TS_STEP_COMPLETE: 587 h = ts->ptime - ts->ptime_prev; 588 break; 589 default: 590 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 591 } 592 if (order == tab->order) { 593 if (ark->status == TS_STEP_INCOMPLETE) { 594 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 595 PetscCall(VecCopy(ark->Y[s - 1], X)); 596 } else { /* Use the standard completion formula (bt,b) */ 597 PetscCall(VecCopy(ts->vec_sol, X)); 598 for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 599 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 600 if (ark->imex) { /* Method is IMEX, complete the explicit formula */ 601 for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 602 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 603 } 604 } 605 } else PetscCall(VecCopy(ts->vec_sol, X)); 606 if (done) *done = PETSC_TRUE; 607 PetscFunctionReturn(PETSC_SUCCESS); 608 } else if (order == tab->order - 1) { 609 if (!tab->bembedt) goto unavailable; 610 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 611 PetscCall(VecCopy(ts->vec_sol, X)); 612 for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 613 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 614 for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 615 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 616 } else { /* Rollback and re-complete using (bet-be,be-b) */ 617 PetscCall(VecCopy(ts->vec_sol, X)); 618 for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 619 PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 620 for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 621 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 622 } 623 if (done) *done = PETSC_TRUE; 624 PetscFunctionReturn(PETSC_SUCCESS); 625 } 626 unavailable: 627 if (done) *done = PETSC_FALSE; 628 else 629 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, 630 tab->order, order); 631 PetscFunctionReturn(PETSC_SUCCESS); 632 } 633 634 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) 635 { 636 Vec Udot, Y1, Y2; 637 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 638 PetscReal norm; 639 640 PetscFunctionBegin; 641 PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 642 PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 643 PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 644 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 645 PetscCall(VecSetRandom(Udot, NULL)); 646 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 647 PetscCall(VecAXPY(Y2, -1.0, Y1)); 648 PetscCall(VecAXPY(Y2, -1.0, Udot)); 649 PetscCall(VecNorm(Y2, NORM_2, &norm)); 650 if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 651 *id = PETSC_TRUE; 652 } else { 653 *id = PETSC_FALSE; 654 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)); 655 } 656 PetscCall(VecDestroy(&Udot)); 657 PetscCall(VecDestroy(&Y1)); 658 PetscCall(VecDestroy(&Y2)); 659 PetscFunctionReturn(PETSC_SUCCESS); 660 } 661 662 static PetscErrorCode TSRollBack_ARKIMEX(TS ts) 663 { 664 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 665 ARKTableau tab = ark->tableau; 666 const PetscInt s = tab->s; 667 const PetscReal *bt = tab->bt, *b = tab->b; 668 PetscScalar *w = ark->work; 669 Vec *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS; 670 PetscInt j; 671 PetscReal h; 672 673 PetscFunctionBegin; 674 switch (ark->status) { 675 case TS_STEP_INCOMPLETE: 676 case TS_STEP_PENDING: 677 h = ts->time_step; 678 break; 679 case TS_STEP_COMPLETE: 680 h = ts->ptime - ts->ptime_prev; 681 break; 682 default: 683 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 684 } 685 for (j = 0; j < s; j++) w[j] = -h * bt[j]; 686 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotI)); 687 for (j = 0; j < s; j++) w[j] = -h * b[j]; 688 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS)); 689 PetscFunctionReturn(PETSC_SUCCESS); 690 } 691 692 static PetscErrorCode TSStep_ARKIMEX(TS ts) 693 { 694 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 695 ARKTableau tab = ark->tableau; 696 const PetscInt s = tab->s; 697 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 698 PetscScalar *w = ark->work; 699 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 700 PetscBool extrapolate = ark->extrapolate; 701 TSAdapt adapt; 702 SNES snes; 703 PetscInt i, j, its, lits; 704 PetscInt rejections = 0; 705 PetscBool stageok, accept = PETSC_TRUE; 706 PetscReal next_time_step = ts->time_step; 707 708 PetscFunctionBegin; 709 if (ark->extrapolate && !ark->Y_prev) { 710 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 711 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 712 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 713 } 714 715 if (!ts->steprollback) { 716 if (ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 717 PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 718 } 719 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 720 for (i = 0; i < s; i++) { 721 PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 722 PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 723 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 724 } 725 } 726 } 727 728 if (ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 729 TS ts_start; 730 if (PetscDefined(USE_DEBUG)) { 731 PetscBool id = PETSC_FALSE; 732 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 733 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "This scheme requires an identity mass matrix, however the TSIFunction you provide does not utilize an identity mass matrix"); 734 } 735 PetscCall(TSClone(ts, &ts_start)); 736 PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 737 PetscCall(TSSetTime(ts_start, ts->ptime)); 738 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 739 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 740 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 741 PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 742 PetscCall(TSSetType(ts_start, TSARKIMEX)); 743 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 744 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 745 746 PetscCall(TSRestartStep(ts_start)); 747 PetscCall(TSSolve(ts_start, ts->vec_sol)); 748 PetscCall(TSGetTime(ts_start, &ts->ptime)); 749 PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 750 751 { /* Save the initial slope for the next step */ 752 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 753 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 754 } 755 ts->steps++; 756 757 /* Set the correct TS in SNES */ 758 /* We'll try to bypass this by changing the method on the fly */ 759 { 760 PetscCall(TSGetSNES(ts, &snes)); 761 PetscCall(TSSetSNES(ts, snes)); 762 } 763 PetscCall(TSDestroy(&ts_start)); 764 } 765 766 ark->status = TS_STEP_INCOMPLETE; 767 while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 768 PetscReal t = ts->ptime; 769 PetscReal h = ts->time_step; 770 for (i = 0; i < s; i++) { 771 ark->stage_time = t + h * ct[i]; 772 PetscCall(TSPreStage(ts, ark->stage_time)); 773 if (At[i * s + i] == 0) { /* This stage is explicit */ 774 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"); 775 PetscCall(VecCopy(ts->vec_sol, Y[i])); 776 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 777 PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 778 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 779 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 780 } else { 781 ark->scoeff = 1. / At[i * s + i]; 782 /* Ydot = shift*(Y-Z) */ 783 PetscCall(VecCopy(ts->vec_sol, Z)); 784 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 785 PetscCall(VecMAXPY(Z, i, w, YdotI)); 786 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 787 PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 788 if (extrapolate && !ts->steprestart) { 789 /* Initial guess extrapolated from previous time step stage values */ 790 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 791 } else { 792 /* Initial guess taken from last stage */ 793 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 794 } 795 PetscCall(TSGetSNES(ts, &snes)); 796 PetscCall(SNESSolve(snes, NULL, Y[i])); 797 PetscCall(SNESGetIterationNumber(snes, &its)); 798 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 799 ts->snes_its += its; 800 ts->ksp_its += lits; 801 PetscCall(TSGetAdapt(ts, &adapt)); 802 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 803 if (!stageok) { 804 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 805 * use extrapolation to initialize the solves on the next attempt. */ 806 extrapolate = PETSC_FALSE; 807 goto reject_step; 808 } 809 } 810 if (ts->equation_type >= TS_EQ_IMPLICIT) { 811 if (i == 0 && tab->explicit_first_stage) { 812 PetscCheck(tab->stiffly_accurate, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %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", 813 ark->tableau->name); 814 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 815 } else { 816 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 817 } 818 } else { 819 if (i == 0 && tab->explicit_first_stage) { 820 PetscCall(VecZeroEntries(Ydot)); 821 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 822 PetscCall(VecScale(YdotI[i], -1.0)); 823 } else { 824 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 825 } 826 if (ark->imex) { 827 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 828 } else { 829 PetscCall(VecZeroEntries(YdotRHS[i])); 830 } 831 } 832 PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 833 } 834 835 ark->status = TS_STEP_INCOMPLETE; 836 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 837 ark->status = TS_STEP_PENDING; 838 PetscCall(TSGetAdapt(ts, &adapt)); 839 PetscCall(TSAdaptCandidatesClear(adapt)); 840 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 841 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 842 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 843 if (!accept) { /* Roll back the current step */ 844 PetscCall(TSRollBack_ARKIMEX(ts)); 845 ts->time_step = next_time_step; 846 goto reject_step; 847 } 848 849 ts->ptime += ts->time_step; 850 ts->time_step = next_time_step; 851 break; 852 853 reject_step: 854 ts->reject++; 855 accept = PETSC_FALSE; 856 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 857 ts->reason = TS_DIVERGED_STEP_REJECTED; 858 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 859 } 860 } 861 PetscFunctionReturn(PETSC_SUCCESS); 862 } 863 /* 864 This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form 865 Udot = H(t,U) + G(t,U) 866 This corresponds to F(t,U,Udot) = Udot-H(t,U). 867 868 The complete adjoint equations are 869 (shift*I - dHdu) lambda_s[i] = 1/at[i][i] ( 870 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 871 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0 872 lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j] 873 mu_{n+1}[i] = h (at[i][i] dHdP lambda_s[i] 874 + (b_i dGdP + bt[i] dHdP) lambda_{n+1} + dGdP \sum_{j=i+1}^s a[j][i] lambda_s[j] 875 + dHdP \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0 876 mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j] 877 where shift = 1/(at[i][i]*h) 878 879 If at[i][i] is 0, the first equation falls back to 880 lambda_s[i] = h ( 881 (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j] 882 + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]) 883 884 */ 885 static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts) 886 { 887 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 888 ARKTableau tab = ark->tableau; 889 const PetscInt s = tab->s; 890 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt; 891 PetscScalar *w = ark->work; 892 Vec *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp; 893 Mat Jex, Jim, Jimpre; 894 PetscInt i, j, nadj; 895 PetscReal t = ts->ptime, stage_time_ex; 896 PetscReal adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */ 897 KSP ksp; 898 899 PetscFunctionBegin; 900 ark->status = TS_STEP_INCOMPLETE; 901 PetscCall(SNESGetKSP(ts->snes, &ksp)); 902 PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL)); 903 PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL)); 904 905 for (i = s - 1; i >= 0; i--) { 906 ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]); 907 stage_time_ex = t - adjoint_time_step * (1.0 - c[i]); 908 if (At[i * s + i] == 0) { // This stage is explicit 909 ark->scoeff = 0.; 910 } else { 911 ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative 912 } 913 PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre)); 914 PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex)); 915 if (ts->vecs_sensip) { 916 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 917 PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs)); // get dGdP 918 } 919 for (nadj = 0; nadj < ts->numcost; nadj++) { 920 /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */ 921 if (s - i - 1 > 0) { 922 /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */ 923 for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i]; 924 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 925 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 926 /* (shift I - dHdU) Temp */ 927 PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 928 /* cancel out shift Temp where shift=-scoeff/h */ 929 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj])); 930 if (ts->vecs_sensip) { 931 /* - dHdP Temp */ 932 PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 933 /* mu_n += h dHdP Temp */ 934 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsSensiPTemp[nadj])); 935 } 936 /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */ 937 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; 938 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 939 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 940 /* dGdU Temp */ 941 PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 942 if (ts->vecs_sensip) { 943 /* dGdP Temp */ 944 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj])); 945 /* mu_n += h dGdP Temp */ 946 PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj])); 947 } 948 } else { 949 PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized 950 } 951 if (bt[i]) { // bt[i] dHdU lambda_{n+1} 952 /* (shift I - dHdU)^T lambda_{n+1} */ 953 PetscCall(MatMultTranspose(Jim, ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); 954 /* -bt[i] (shift I - dHdU)^T lambda_{n+1} */ 955 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i], VecsSensiTemp[nadj])); 956 /* cancel out -bt[i] shift lambda_{n+1} */ 957 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i] * ark->scoeff / adjoint_time_step, ts->vecs_sensi[nadj])); 958 if (ts->vecs_sensip) { 959 /* -dHdP lambda_{n+1} */ 960 PetscCall(MatMultTranspose(ts->Jacp, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj])); 961 /* mu_n += h bt[i] dHdP lambda_{n+1} */ 962 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -bt[i] * adjoint_time_step, VecsSensiPTemp[nadj])); 963 } 964 } 965 if (b[i]) { // b[i] dGdU lambda_{n+1} 966 /* dGdU lambda_{n+1} */ 967 PetscCall(MatMultTranspose(Jex, ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); 968 /* b[i] dGdU lambda_{n+1} */ 969 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], b[i], VecsSensiTemp[nadj])); 970 if (ts->vecs_sensip) { 971 /* dGdP lambda_{n+1} */ 972 PetscCall(MatMultTranspose(ts->Jacprhs, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj])); 973 /* mu_n += h b[i] dGdP lambda_{n+1} */ 974 PetscCall(VecAXPY(ts->vecs_sensip[nadj], b[i] * adjoint_time_step, VecsSensiPTemp[nadj])); 975 } 976 } 977 /* Build LHS for first-order adjoint */ 978 if (At[i * s + i] == 0) { // This stage is explicit 979 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step)); 980 } else { 981 KSP ksp; 982 KSPConvergedReason kspreason; 983 PetscCall(SNESGetKSP(ts->snes, &ksp)); 984 PetscCall(KSPSetOperators(ksp, Jim, Jimpre)); 985 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i])); 986 PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i])); 987 PetscCall(KSPGetConvergedReason(ksp, &kspreason)); 988 if (kspreason < 0) { 989 ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE; 990 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj)); 991 } 992 if (ts->vecs_sensip) { 993 /* -dHdP lambda_s[i] */ 994 PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj])); 995 /* mu_n += h at[i][i] dHdP lambda_s[i] */ 996 PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj])); 997 } 998 } 999 } 1000 } 1001 for (j = 0; j < s; j++) w[j] = 1.0; 1002 for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's 1003 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 1004 ark->status = TS_STEP_COMPLETE; 1005 PetscFunctionReturn(PETSC_SUCCESS); 1006 } 1007 1008 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) 1009 { 1010 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1011 PetscInt s = ark->tableau->s, pinterp = ark->tableau->pinterp, i, j; 1012 PetscReal h; 1013 PetscReal tt, t; 1014 PetscScalar *bt, *b; 1015 const PetscReal *Bt = ark->tableau->binterpt, *B = ark->tableau->binterp; 1016 1017 PetscFunctionBegin; 1018 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 1019 switch (ark->status) { 1020 case TS_STEP_INCOMPLETE: 1021 case TS_STEP_PENDING: 1022 h = ts->time_step; 1023 t = (itime - ts->ptime) / h; 1024 break; 1025 case TS_STEP_COMPLETE: 1026 h = ts->ptime - ts->ptime_prev; 1027 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1028 break; 1029 default: 1030 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1031 } 1032 PetscCall(PetscMalloc2(s, &bt, s, &b)); 1033 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 1034 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1035 for (i = 0; i < s; i++) { 1036 bt[i] += h * Bt[i * pinterp + j] * tt; 1037 b[i] += h * B[i * pinterp + j] * tt; 1038 } 1039 } 1040 PetscCall(VecCopy(ark->Y[0], X)); 1041 PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 1042 PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 1043 PetscCall(PetscFree2(bt, b)); 1044 PetscFunctionReturn(PETSC_SUCCESS); 1045 } 1046 1047 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) 1048 { 1049 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1050 PetscInt s = ark->tableau->s, pinterp = ark->tableau->pinterp, i, j; 1051 PetscReal h, h_prev, t, tt; 1052 PetscScalar *bt, *b; 1053 const PetscReal *Bt = ark->tableau->binterpt, *B = ark->tableau->binterp; 1054 1055 PetscFunctionBegin; 1056 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 1057 PetscCall(PetscCalloc2(s, &bt, s, &b)); 1058 h = ts->time_step; 1059 h_prev = ts->ptime - ts->ptime_prev; 1060 t = 1 + h / h_prev * c; 1061 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1062 for (i = 0; i < s; i++) { 1063 bt[i] += h * Bt[i * pinterp + j] * tt; 1064 b[i] += h * B[i * pinterp + j] * tt; 1065 } 1066 } 1067 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 1068 PetscCall(VecCopy(ark->Y_prev[0], X)); 1069 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 1070 PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 1071 PetscCall(PetscFree2(bt, b)); 1072 PetscFunctionReturn(PETSC_SUCCESS); 1073 } 1074 1075 /*------------------------------------------------------------*/ 1076 1077 static PetscErrorCode TSARKIMEXTableauReset(TS ts) 1078 { 1079 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1080 ARKTableau tab = ark->tableau; 1081 1082 PetscFunctionBegin; 1083 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1084 PetscCall(PetscFree(ark->work)); 1085 PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 1086 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 1087 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 1088 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 1089 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 1090 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 1091 PetscFunctionReturn(PETSC_SUCCESS); 1092 } 1093 1094 static PetscErrorCode TSReset_ARKIMEX(TS ts) 1095 { 1096 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1097 1098 PetscFunctionBegin; 1099 PetscCall(TSARKIMEXTableauReset(ts)); 1100 PetscCall(VecDestroy(&ark->Ydot)); 1101 PetscCall(VecDestroy(&ark->Ydot0)); 1102 PetscCall(VecDestroy(&ark->Z)); 1103 PetscFunctionReturn(PETSC_SUCCESS); 1104 } 1105 1106 static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts) 1107 { 1108 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1109 ARKTableau tab = ark->tableau; 1110 1111 PetscFunctionBegin; 1112 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam)); 1113 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp)); 1114 PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp)); 1115 PetscFunctionReturn(PETSC_SUCCESS); 1116 } 1117 1118 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1119 { 1120 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 1121 1122 PetscFunctionBegin; 1123 if (Z) { 1124 if (dm && dm != ts->dm) { 1125 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1126 } else *Z = ax->Z; 1127 } 1128 if (Ydot) { 1129 if (dm && dm != ts->dm) { 1130 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1131 } else *Ydot = ax->Ydot; 1132 } 1133 PetscFunctionReturn(PETSC_SUCCESS); 1134 } 1135 1136 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) 1137 { 1138 PetscFunctionBegin; 1139 if (Z) { 1140 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 1141 } 1142 if (Ydot) { 1143 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 1144 } 1145 PetscFunctionReturn(PETSC_SUCCESS); 1146 } 1147 1148 /* 1149 This defines the nonlinear equation that is to be solved with SNES 1150 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1151 */ 1152 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) 1153 { 1154 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1155 DM dm, dmsave; 1156 Vec Z, Ydot; 1157 PetscReal shift = ark->scoeff / ts->time_step; 1158 1159 PetscFunctionBegin; 1160 PetscCall(SNESGetDM(snes, &dm)); 1161 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 1162 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 1163 dmsave = ts->dm; 1164 ts->dm = dm; 1165 1166 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 1167 1168 ts->dm = dmsave; 1169 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 1170 PetscFunctionReturn(PETSC_SUCCESS); 1171 } 1172 1173 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) 1174 { 1175 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1176 DM dm, dmsave; 1177 Vec Ydot; 1178 PetscReal shift = ark->scoeff / ts->time_step; 1179 1180 PetscFunctionBegin; 1181 PetscCall(SNESGetDM(snes, &dm)); 1182 PetscCall(TSARKIMEXGetVecs(ts, dm, NULL, &Ydot)); 1183 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1184 dmsave = ts->dm; 1185 ts->dm = dm; 1186 1187 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 1188 1189 ts->dm = dmsave; 1190 PetscCall(TSARKIMEXRestoreVecs(ts, dm, NULL, &Ydot)); 1191 PetscFunctionReturn(PETSC_SUCCESS); 1192 } 1193 1194 static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[]) 1195 { 1196 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1197 1198 PetscFunctionBegin; 1199 if (ns) *ns = ark->tableau->s; 1200 if (Y) *Y = ark->Y; 1201 PetscFunctionReturn(PETSC_SUCCESS); 1202 } 1203 1204 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) 1205 { 1206 PetscFunctionBegin; 1207 PetscFunctionReturn(PETSC_SUCCESS); 1208 } 1209 1210 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1211 { 1212 TS ts = (TS)ctx; 1213 Vec Z, Z_c; 1214 1215 PetscFunctionBegin; 1216 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 1217 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 1218 PetscCall(MatRestrict(restrct, Z, Z_c)); 1219 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 1220 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 1221 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 1222 PetscFunctionReturn(PETSC_SUCCESS); 1223 } 1224 1225 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) 1226 { 1227 PetscFunctionBegin; 1228 PetscFunctionReturn(PETSC_SUCCESS); 1229 } 1230 1231 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1232 { 1233 TS ts = (TS)ctx; 1234 Vec Z, Z_c; 1235 1236 PetscFunctionBegin; 1237 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 1238 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 1239 1240 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1241 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1242 1243 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 1244 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 1245 PetscFunctionReturn(PETSC_SUCCESS); 1246 } 1247 1248 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) 1249 { 1250 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1251 ARKTableau tab = ark->tableau; 1252 1253 PetscFunctionBegin; 1254 PetscCall(PetscMalloc1(tab->s, &ark->work)); 1255 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 1256 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 1257 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 1258 if (ark->extrapolate) { 1259 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1260 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1261 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1262 } 1263 PetscFunctionReturn(PETSC_SUCCESS); 1264 } 1265 1266 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 1267 { 1268 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1269 DM dm; 1270 SNES snes; 1271 1272 PetscFunctionBegin; 1273 PetscCall(TSARKIMEXTableauSetUp(ts)); 1274 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 1275 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 1276 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 1277 PetscCall(TSGetDM(ts, &dm)); 1278 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 1279 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 1280 PetscCall(TSGetSNES(ts, &snes)); 1281 PetscFunctionReturn(PETSC_SUCCESS); 1282 } 1283 /*------------------------------------------------------------*/ 1284 1285 static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts) 1286 { 1287 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1288 ARKTableau tab = ark->tableau; 1289 1290 PetscFunctionBegin; 1291 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam)); 1292 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp)); 1293 if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); } 1294 if (PetscDefined(USE_DEBUG)) { 1295 PetscBool id = PETSC_FALSE; 1296 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 1297 PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Adjoint ARKIMEX requires an identity mass matrix, however the TSIFunction you provide does not utilize an identity mass matrix"); 1298 } 1299 PetscFunctionReturn(PETSC_SUCCESS); 1300 } 1301 1302 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject) 1303 { 1304 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1305 1306 PetscFunctionBegin; 1307 PetscOptionsHeadBegin(PetscOptionsObject, "ARKIMEX ODE solver options"); 1308 { 1309 ARKTableauLink link; 1310 PetscInt count, choice; 1311 PetscBool flg; 1312 const char **namelist; 1313 for (link = ARKTableauList, count = 0; link; link = link->next, count++) 1314 ; 1315 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1316 for (link = ARKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1317 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 1318 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 1319 PetscCall(PetscFree(namelist)); 1320 1321 flg = (PetscBool)!ark->imex; 1322 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 1323 ark->imex = (PetscBool)!flg; 1324 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)); 1325 } 1326 PetscOptionsHeadEnd(); 1327 PetscFunctionReturn(PETSC_SUCCESS); 1328 } 1329 1330 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) 1331 { 1332 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1333 PetscBool iascii; 1334 1335 PetscFunctionBegin; 1336 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1337 if (iascii) { 1338 ARKTableau tab = ark->tableau; 1339 TSARKIMEXType arktype; 1340 char buf[512]; 1341 PetscBool flg; 1342 1343 PetscCall(TSARKIMEXGetType(ts, &arktype)); 1344 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 1345 PetscCall(PetscViewerASCIIPrintf(viewer, " ARK IMEX %s\n", arktype)); 1346 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 1347 PetscCall(PetscViewerASCIIPrintf(viewer, " Stiff abscissa ct = %s\n", buf)); 1348 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 1349 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 1350 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 1351 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 1352 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 1353 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 1354 } 1355 PetscFunctionReturn(PETSC_SUCCESS); 1356 } 1357 1358 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) 1359 { 1360 SNES snes; 1361 TSAdapt adapt; 1362 1363 PetscFunctionBegin; 1364 PetscCall(TSGetAdapt(ts, &adapt)); 1365 PetscCall(TSAdaptLoad(adapt, viewer)); 1366 PetscCall(TSGetSNES(ts, &snes)); 1367 PetscCall(SNESLoad(snes, viewer)); 1368 /* function and Jacobian context for SNES when used with TS is always ts object */ 1369 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 1370 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 1371 PetscFunctionReturn(PETSC_SUCCESS); 1372 } 1373 1374 /*@C 1375 TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme 1376 1377 Logically Collective 1378 1379 Input Parameters: 1380 + ts - timestepping context 1381 - arktype - type of `TSARKIMEX` scheme 1382 1383 Options Database Key: 1384 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type 1385 1386 Level: intermediate 1387 1388 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 1389 `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 1390 @*/ 1391 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) 1392 { 1393 PetscFunctionBegin; 1394 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1395 PetscValidCharPointer(arktype, 2); 1396 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 1397 PetscFunctionReturn(PETSC_SUCCESS); 1398 } 1399 1400 /*@C 1401 TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme 1402 1403 Logically Collective 1404 1405 Input Parameter: 1406 . ts - timestepping context 1407 1408 Output Parameter: 1409 . arktype - type of `TSARKIMEX` scheme 1410 1411 Level: intermediate 1412 1413 .seealso: [](ch_ts), `TSARKIMEX`c, `TSARKIMEXGetType()` 1414 @*/ 1415 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) 1416 { 1417 PetscFunctionBegin; 1418 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1419 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 1420 PetscFunctionReturn(PETSC_SUCCESS); 1421 } 1422 1423 /*@ 1424 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly 1425 1426 Logically Collective 1427 1428 Input Parameters: 1429 + ts - timestepping context 1430 - flg - `PETSC_TRUE` for fully implicit 1431 1432 Level: intermediate 1433 1434 .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 1435 @*/ 1436 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) 1437 { 1438 PetscFunctionBegin; 1439 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1440 PetscValidLogicalCollectiveBool(ts, flg, 2); 1441 PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 1442 PetscFunctionReturn(PETSC_SUCCESS); 1443 } 1444 1445 /*@ 1446 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 1447 1448 Logically Collective 1449 1450 Input Parameter: 1451 . ts - timestepping context 1452 1453 Output Parameter: 1454 . flg - `PETSC_TRUE` for fully implicit 1455 1456 Level: intermediate 1457 1458 .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 1459 @*/ 1460 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) 1461 { 1462 PetscFunctionBegin; 1463 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1464 PetscValidBoolPointer(flg, 2); 1465 PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 1466 PetscFunctionReturn(PETSC_SUCCESS); 1467 } 1468 1469 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) 1470 { 1471 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1472 1473 PetscFunctionBegin; 1474 *arktype = ark->tableau->name; 1475 PetscFunctionReturn(PETSC_SUCCESS); 1476 } 1477 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) 1478 { 1479 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1480 PetscBool match; 1481 ARKTableauLink link; 1482 1483 PetscFunctionBegin; 1484 if (ark->tableau) { 1485 PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 1486 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1487 } 1488 for (link = ARKTableauList; link; link = link->next) { 1489 PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 1490 if (match) { 1491 if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 1492 ark->tableau = &link->tab; 1493 if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 1494 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1495 PetscFunctionReturn(PETSC_SUCCESS); 1496 } 1497 } 1498 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 1499 } 1500 1501 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) 1502 { 1503 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1504 1505 PetscFunctionBegin; 1506 ark->imex = (PetscBool)!flg; 1507 PetscFunctionReturn(PETSC_SUCCESS); 1508 } 1509 1510 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) 1511 { 1512 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1513 1514 PetscFunctionBegin; 1515 *flg = (PetscBool)!ark->imex; 1516 PetscFunctionReturn(PETSC_SUCCESS); 1517 } 1518 1519 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 1520 { 1521 PetscFunctionBegin; 1522 PetscCall(TSReset_ARKIMEX(ts)); 1523 if (ts->dm) { 1524 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 1525 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 1526 } 1527 PetscCall(PetscFree(ts->data)); 1528 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 1529 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 1530 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 1531 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 1532 PetscFunctionReturn(PETSC_SUCCESS); 1533 } 1534 1535 /* ------------------------------------------------------------ */ 1536 /*MC 1537 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 1538 1539 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1540 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1541 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 1542 1543 Level: beginner 1544 1545 Notes: 1546 The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type 1547 1548 If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information. 1549 1550 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). 1551 1552 Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear. 1553 1554 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 1555 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 1556 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType` 1557 M*/ 1558 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 1559 { 1560 TS_ARKIMEX *ark; 1561 1562 PetscFunctionBegin; 1563 PetscCall(TSARKIMEXInitializePackage()); 1564 1565 ts->ops->reset = TSReset_ARKIMEX; 1566 ts->ops->adjointreset = TSAdjointReset_ARKIMEX; 1567 ts->ops->destroy = TSDestroy_ARKIMEX; 1568 ts->ops->view = TSView_ARKIMEX; 1569 ts->ops->load = TSLoad_ARKIMEX; 1570 ts->ops->setup = TSSetUp_ARKIMEX; 1571 ts->ops->adjointsetup = TSAdjointSetUp_ARKIMEX; 1572 ts->ops->step = TSStep_ARKIMEX; 1573 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1574 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1575 ts->ops->rollback = TSRollBack_ARKIMEX; 1576 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1577 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1578 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1579 ts->ops->getstages = TSGetStages_ARKIMEX; 1580 ts->ops->adjointstep = TSAdjointStep_ARKIMEX; 1581 1582 ts->usessnes = PETSC_TRUE; 1583 1584 PetscCall(PetscNew(&ark)); 1585 ts->data = (void *)ark; 1586 ark->imex = PETSC_TRUE; 1587 1588 ark->VecsDeltaLam = NULL; 1589 ark->VecsSensiTemp = NULL; 1590 ark->VecsSensiPTemp = NULL; 1591 1592 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 1593 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 1594 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 1595 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 1596 1597 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 1598 PetscFunctionReturn(PETSC_SUCCESS); 1599 } 1600