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 } TS_ARKIMEX; 60 /*MC 61 TSARKIMEXARS122 - Second order ARK IMEX scheme. 62 63 This method has one explicit stage and one implicit stage. 64 65 Options Database: 66 . -ts_arkimex_type ars122 - set arkimex type to ars122 67 68 References: 69 . * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 70 71 Level: advanced 72 73 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 74 M*/ 75 /*MC 76 TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part. 77 78 This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu. 79 80 Options Database: 81 . -ts_arkimex_type a2 - set arkimex type to a2 82 83 Level: advanced 84 85 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 86 M*/ 87 /*MC 88 TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part. 89 90 This method has two implicit stages, and L-stable implicit scheme. 91 92 Options Database: 93 . -ts_arkimex_type l2 - set arkimex type to l2 94 95 References: 96 . * - 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. 97 98 Level: advanced 99 100 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 101 M*/ 102 /*MC 103 TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method. 104 105 This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used. 106 107 Options Database: 108 . -ts_arkimex_type 1bee - set arkimex type to 1bee 109 110 Level: advanced 111 112 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 113 M*/ 114 /*MC 115 TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part. 116 117 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. 118 119 Options Database: 120 . -ts_arkimex_type 2c - set arkimex type to 2c 121 122 Level: advanced 123 124 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 125 M*/ 126 /*MC 127 TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 128 129 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 implict component. This method was provided by Emil Constantinescu. 130 131 Options Database: 132 . -ts_arkimex_type 2d - set arkimex type to 2d 133 134 Level: advanced 135 136 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 137 M*/ 138 /*MC 139 TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 140 141 This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu. 142 143 Options Database: 144 . -ts_arkimex_type 2e - set arkimex type to 2e 145 146 Level: advanced 147 148 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 149 M*/ 150 /*MC 151 TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme. 152 153 This method has three implicit stages. 154 155 References: 156 . * - 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. 157 158 This method is referred to as SSP2-(3,3,2) in https://arxiv.org/abs/1110.4375 159 160 Options Database: 161 . -ts_arkimex_type prssp2 - set arkimex type to prssp2 162 163 Level: advanced 164 165 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 166 M*/ 167 /*MC 168 TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part. 169 170 This method has one explicit stage and three implicit stages. 171 172 Options Database: 173 . -ts_arkimex_type 3 - set arkimex type to 3 174 175 References: 176 . * - Kennedy and Carpenter 2003. 177 178 Level: advanced 179 180 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 181 M*/ 182 /*MC 183 TSARKIMEXARS443 - Third order ARK IMEX scheme. 184 185 This method has one explicit stage and four implicit stages. 186 187 Options Database: 188 . -ts_arkimex_type ars443 - set arkimex type to ars443 189 190 References: 191 + * - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997). 192 - * - This method is referred to as ARS(4,4,3) in https://arxiv.org/abs/1110.4375 193 194 Level: advanced 195 196 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 197 M*/ 198 /*MC 199 TSARKIMEXBPR3 - Third order ARK IMEX scheme. 200 201 This method has one explicit stage and four implicit stages. 202 203 Options Database: 204 . -ts_arkimex_type bpr3 - set arkimex type to bpr3 205 206 References: 207 . * - This method is referred to as ARK3 in https://arxiv.org/abs/1110.4375 208 209 Level: advanced 210 211 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 212 M*/ 213 /*MC 214 TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part. 215 216 This method has one explicit stage and four implicit stages. 217 218 Options Database: 219 . -ts_arkimex_type 4 - set arkimex type to4 220 221 References: 222 . * - Kennedy and Carpenter 2003. 223 224 Level: advanced 225 226 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 227 M*/ 228 /*MC 229 TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part. 230 231 This method has one explicit stage and five implicit stages. 232 233 Options Database: 234 . -ts_arkimex_type 5 - set arkimex type to 5 235 236 References: 237 . * - Kennedy and Carpenter 2003. 238 239 Level: advanced 240 241 .seealso: `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()` 242 M*/ 243 244 /*@C 245 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSARKIMEX 246 247 Not Collective, but should be called by all processes which will need the schemes to be registered 248 249 Level: advanced 250 251 .seealso: `TSARKIMEXRegisterDestroy()` 252 @*/ 253 PetscErrorCode TSARKIMEXRegisterAll(void) { 254 PetscFunctionBegin; 255 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(0); 256 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 257 258 { 259 const PetscReal A[3][3] = 260 { 261 {0.0, 0.0, 0.0}, 262 {0.0, 0.0, 0.0}, 263 {0.0, 0.5, 0.0} 264 }, 265 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}; 266 PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 267 } 268 { 269 const PetscReal A[2][2] = 270 { 271 {0.0, 0.0}, 272 {0.5, 0.0} 273 }, 274 At[2][2] = {{0.0, 0.0}, {0.0, 0.5}}, b[2] = {0.0, 1.0}, bembedt[2] = {0.5, 0.5}; 275 /* 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*/ 276 PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 277 } 278 { 279 const PetscReal A[2][2] = 280 { 281 {0.0, 0.0}, 282 {1.0, 0.0} 283 }, 284 At[2][2] = {{0.0, 0.0}, {0.5, 0.5}}, b[2] = {0.5, 0.5}, bembedt[2] = {0.0, 1.0}; 285 /* 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*/ 286 PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL)); 287 } 288 { 289 /* 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 */ 290 const PetscReal A[2][2] = 291 { 292 {0.0, 0.0}, 293 {1.0, 0.0} 294 }, 295 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}}; 296 PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0])); 297 } 298 { 299 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */ 300 const PetscReal A[3][3] = 301 { 302 {0, 0, 0}, 303 {2 - 1.414213562373095048802, 0, 0}, 304 {0.5, 0.5, 0} 305 }, 306 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}}; 307 PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 308 } 309 { 310 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */ 311 const PetscReal A[3][3] = 312 { 313 {0, 0, 0}, 314 {2 - 1.414213562373095048802, 0, 0}, 315 {0.75, 0.25, 0} 316 }, 317 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}}; 318 PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 319 } 320 { /* Optimal for linear implicit part */ 321 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */ 322 const PetscReal A[3][3] = 323 { 324 {0, 0, 0}, 325 {2 - 1.414213562373095048802, 0, 0}, 326 {(3 - 2 * 1.414213562373095048802) / 6, (3 + 2 * 1.414213562373095048802) / 6, 0} 327 }, 328 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}}; 329 PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 330 } 331 { /* Optimal for linear implicit part */ 332 const PetscReal A[3][3] = 333 { 334 {0, 0, 0}, 335 {0.5, 0, 0}, 336 {0.5, 0.5, 0} 337 }, 338 At[3][3] = {{0.25, 0, 0}, {0, 0.25, 0}, {1. / 3, 1. / 3, 1. / 3}}; 339 PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL)); 340 } 341 { 342 const PetscReal A[4][4] = 343 { 344 {0, 0, 0, 0}, 345 {1767732205903. / 2027836641118., 0, 0, 0}, 346 {5535828885825. / 10492691773637., 788022342437. / 10882634858940., 0, 0}, 347 {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0} 348 }, 349 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.}}; 350 PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL)); 351 } 352 { 353 const PetscReal A[5][5] = 354 { 355 {0, 0, 0, 0, 0}, 356 {1. / 2, 0, 0, 0, 0}, 357 {11. / 18, 1. / 18, 0, 0, 0}, 358 {5. / 6, -5. / 6, .5, 0, 0}, 359 {1. / 4, 7. / 4, 3. / 4, -7. / 4, 0} 360 }, 361 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; 362 PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 0, NULL, NULL)); 363 } 364 { 365 const PetscReal A[5][5] = 366 { 367 {0, 0, 0, 0, 0}, 368 {1, 0, 0, 0, 0}, 369 {4. / 9, 2. / 9, 0, 0, 0}, 370 {1. / 4, 0, 3. / 4, 0, 0}, 371 {1. / 4, 0, 3. / 5, 0, 0} 372 }, 373 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; 374 PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 0, NULL, NULL)); 375 } 376 { 377 const PetscReal A[6][6] = 378 { 379 {0, 0, 0, 0, 0, 0}, 380 {1. / 2, 0, 0, 0, 0, 0}, 381 {13861. / 62500., 6889. / 62500., 0, 0, 0, 0}, 382 {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209., 0, 0, 0}, 383 {-451086348788. / 2902428689909., -2682348792572. / 7519795681897., 12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0, 0}, 384 {647845179188. / 3216320057751., 73281519250. / 8382639484533., 552539513391. / 3454668386233., 3354512671639. / 8306763924573., 4040. / 17871., 0} 385 }, 386 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}, 387 {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035.}, {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.}, 388 {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469.}, {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.}}; 389 PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 390 } 391 { 392 const PetscReal A[8][8] = 393 { 394 {0, 0, 0, 0, 0, 0, 0, 0}, 395 {41. / 100, 0, 0, 0, 0, 0, 0, 0}, 396 {367902744464. / 2072280473677., 677623207551. / 8224143866563., 0, 0, 0, 0, 0, 0}, 397 {1268023523408. / 10340822734521., 0, 1029933939417. / 13636558850479., 0, 0, 0, 0, 0}, 398 {14463281900351. / 6315353703477., 0, 66114435211212. / 5879490589093., -54053170152839. / 4284798021562., 0, 0, 0, 0}, 399 {14090043504691. / 34967701212078., 0, 15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0, 0, 0}, 400 {19230459214898. / 13134317526959., 0, 21275331358303. / 2942455364971., -38145345988419. / 4862620318723., -1. / 8, -1. / 8, 0, 0}, 401 {-19977161125411. / 11928030595625., 0, -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261., -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0} 402 }, 403 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.}}; 404 PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL)); 405 } 406 PetscFunctionReturn(0); 407 } 408 409 /*@C 410 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister(). 411 412 Not Collective 413 414 Level: advanced 415 416 .seealso: `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()` 417 @*/ 418 PetscErrorCode TSARKIMEXRegisterDestroy(void) { 419 ARKTableauLink link; 420 421 PetscFunctionBegin; 422 while ((link = ARKTableauList)) { 423 ARKTableau t = &link->tab; 424 ARKTableauList = link->next; 425 PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c)); 426 PetscCall(PetscFree2(t->bembedt, t->bembed)); 427 PetscCall(PetscFree2(t->binterpt, t->binterp)); 428 PetscCall(PetscFree(t->name)); 429 PetscCall(PetscFree(link)); 430 } 431 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 432 PetscFunctionReturn(0); 433 } 434 435 /*@C 436 TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called 437 from TSInitializePackage(). 438 439 Level: developer 440 441 .seealso: `PetscInitialize()` 442 @*/ 443 PetscErrorCode TSARKIMEXInitializePackage(void) { 444 PetscFunctionBegin; 445 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0); 446 TSARKIMEXPackageInitialized = PETSC_TRUE; 447 PetscCall(TSARKIMEXRegisterAll()); 448 PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage)); 449 PetscFunctionReturn(0); 450 } 451 452 /*@C 453 TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is 454 called from PetscFinalize(). 455 456 Level: developer 457 458 .seealso: `PetscFinalize()` 459 @*/ 460 PetscErrorCode TSARKIMEXFinalizePackage(void) { 461 PetscFunctionBegin; 462 TSARKIMEXPackageInitialized = PETSC_FALSE; 463 PetscCall(TSARKIMEXRegisterDestroy()); 464 PetscFunctionReturn(0); 465 } 466 467 /*@C 468 TSARKIMEXRegister - register an ARK IMEX scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 469 470 Not Collective, but the same schemes should be registered on all processes on which they will be used 471 472 Input Parameters: 473 + name - identifier for method 474 . order - approximation order of method 475 . s - number of stages, this is the dimension of the matrices below 476 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 477 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At) 478 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At) 479 . A - Non-stiff stage coefficients (dimension s*s, row-major) 480 . b - Non-stiff step completion table (dimension s; NULL to use last row of At) 481 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A) 482 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available) 483 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided) 484 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 485 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 486 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt) 487 488 Notes: 489 Several ARK IMEX methods are provided, this function is only needed to create new methods. 490 491 Level: advanced 492 493 .seealso: `TSARKIMEX` 494 @*/ 495 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[]) { 496 ARKTableauLink link; 497 ARKTableau t; 498 PetscInt i, j; 499 500 PetscFunctionBegin; 501 PetscCall(TSARKIMEXInitializePackage()); 502 PetscCall(PetscNew(&link)); 503 t = &link->tab; 504 PetscCall(PetscStrallocpy(name, &t->name)); 505 t->order = order; 506 t->s = s; 507 PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c)); 508 PetscCall(PetscArraycpy(t->At, At, s * s)); 509 PetscCall(PetscArraycpy(t->A, A, s * s)); 510 if (bt) PetscCall(PetscArraycpy(t->bt, bt, s)); 511 else 512 for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i]; 513 if (b) PetscCall(PetscArraycpy(t->b, b, s)); 514 else 515 for (i = 0; i < s; i++) t->b[i] = t->bt[i]; 516 if (ct) PetscCall(PetscArraycpy(t->ct, ct, s)); 517 else 518 for (i = 0; i < s; i++) 519 for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j]; 520 if (c) PetscCall(PetscArraycpy(t->c, c, s)); 521 else 522 for (i = 0; i < s; i++) 523 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; 524 t->stiffly_accurate = PETSC_TRUE; 525 for (i = 0; i < s; i++) 526 if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE; 527 t->explicit_first_stage = PETSC_TRUE; 528 for (i = 0; i < s; i++) 529 if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE; 530 /*def of FSAL can be made more precise*/ 531 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate); 532 if (bembedt) { 533 PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed)); 534 PetscCall(PetscArraycpy(t->bembedt, bembedt, s)); 535 PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s)); 536 } 537 538 t->pinterp = pinterp; 539 PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp)); 540 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 541 PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp)); 542 link->next = ARKTableauList; 543 ARKTableauList = link; 544 PetscFunctionReturn(0); 545 } 546 547 /* 548 The step completion formula is 549 550 x1 = x0 - h bt^T YdotI + h b^T YdotRHS 551 552 This function can be called before or after ts->vec_sol has been updated. 553 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 554 We can write 555 556 x1e = x0 - h bet^T YdotI + h be^T YdotRHS 557 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 558 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 559 560 so we can evaluate the method with different order even after the step has been optimistically completed. 561 */ 562 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done) { 563 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 564 ARKTableau tab = ark->tableau; 565 PetscScalar *w = ark->work; 566 PetscReal h; 567 PetscInt s = tab->s, j; 568 569 PetscFunctionBegin; 570 switch (ark->status) { 571 case TS_STEP_INCOMPLETE: 572 case TS_STEP_PENDING: h = ts->time_step; break; 573 case TS_STEP_COMPLETE: h = ts->ptime - ts->ptime_prev; break; 574 default: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 575 } 576 if (order == tab->order) { 577 if (ark->status == TS_STEP_INCOMPLETE) { 578 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */ 579 PetscCall(VecCopy(ark->Y[s - 1], X)); 580 } else { /* Use the standard completion formula (bt,b) */ 581 PetscCall(VecCopy(ts->vec_sol, X)); 582 for (j = 0; j < s; j++) w[j] = h * tab->bt[j]; 583 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 584 if (ark->imex) { /* Method is IMEX, complete the explicit formula */ 585 for (j = 0; j < s; j++) w[j] = h * tab->b[j]; 586 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 587 } 588 } 589 } else PetscCall(VecCopy(ts->vec_sol, X)); 590 if (done) *done = PETSC_TRUE; 591 PetscFunctionReturn(0); 592 } else if (order == tab->order - 1) { 593 if (!tab->bembedt) goto unavailable; 594 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 595 PetscCall(VecCopy(ts->vec_sol, X)); 596 for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j]; 597 PetscCall(VecMAXPY(X, s, w, ark->YdotI)); 598 for (j = 0; j < s; j++) w[j] = h * tab->bembed[j]; 599 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 600 } else { /* Rollback and re-complete using (bet-be,be-b) */ 601 PetscCall(VecCopy(ts->vec_sol, X)); 602 for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]); 603 PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI)); 604 for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]); 605 PetscCall(VecMAXPY(X, s, w, ark->YdotRHS)); 606 } 607 if (done) *done = PETSC_TRUE; 608 PetscFunctionReturn(0); 609 } 610 unavailable: 611 if (done) *done = PETSC_FALSE; 612 else 613 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, 614 tab->order, order); 615 PetscFunctionReturn(0); 616 } 617 618 static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id) { 619 Vec Udot, Y1, Y2; 620 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 621 PetscReal norm; 622 623 PetscFunctionBegin; 624 PetscCall(VecDuplicate(ts->vec_sol, &Udot)); 625 PetscCall(VecDuplicate(ts->vec_sol, &Y1)); 626 PetscCall(VecDuplicate(ts->vec_sol, &Y2)); 627 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex)); 628 PetscCall(VecSetRandom(Udot, NULL)); 629 PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex)); 630 PetscCall(VecAXPY(Y2, -1.0, Y1)); 631 PetscCall(VecAXPY(Y2, -1.0, Udot)); 632 PetscCall(VecNorm(Y2, NORM_2, &norm)); 633 if (norm < 100.0 * PETSC_MACHINE_EPSILON) { 634 *id = PETSC_TRUE; 635 } else { 636 *id = PETSC_FALSE; 637 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)); 638 } 639 PetscCall(VecDestroy(&Udot)); 640 PetscCall(VecDestroy(&Y1)); 641 PetscCall(VecDestroy(&Y2)); 642 PetscFunctionReturn(0); 643 } 644 645 static PetscErrorCode TSRollBack_ARKIMEX(TS ts) { 646 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 647 ARKTableau tab = ark->tableau; 648 const PetscInt s = tab->s; 649 const PetscReal *bt = tab->bt, *b = tab->b; 650 PetscScalar *w = ark->work; 651 Vec *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS; 652 PetscInt j; 653 PetscReal h; 654 655 PetscFunctionBegin; 656 switch (ark->status) { 657 case TS_STEP_INCOMPLETE: 658 case TS_STEP_PENDING: h = ts->time_step; break; 659 case TS_STEP_COMPLETE: h = ts->ptime - ts->ptime_prev; break; 660 default: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 661 } 662 for (j = 0; j < s; j++) w[j] = -h * bt[j]; 663 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotI)); 664 for (j = 0; j < s; j++) w[j] = -h * b[j]; 665 PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS)); 666 PetscFunctionReturn(0); 667 } 668 669 static PetscErrorCode TSStep_ARKIMEX(TS ts) { 670 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 671 ARKTableau tab = ark->tableau; 672 const PetscInt s = tab->s; 673 const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c; 674 PetscScalar *w = ark->work; 675 Vec *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z; 676 PetscBool extrapolate = ark->extrapolate; 677 TSAdapt adapt; 678 SNES snes; 679 PetscInt i, j, its, lits; 680 PetscInt rejections = 0; 681 PetscBool stageok, accept = PETSC_TRUE; 682 PetscReal next_time_step = ts->time_step; 683 684 PetscFunctionBegin; 685 if (ark->extrapolate && !ark->Y_prev) { 686 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 687 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 688 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 689 } 690 691 if (!ts->steprollback) { 692 if (ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */ 693 PetscCall(VecCopy(YdotI[s - 1], Ydot0)); 694 } 695 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */ 696 for (i = 0; i < s; i++) { 697 PetscCall(VecCopy(Y[i], ark->Y_prev[i])); 698 PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i])); 699 PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i])); 700 } 701 } 702 } 703 704 if (ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) { 705 TS ts_start; 706 if (PetscDefined(USE_DEBUG)) { 707 PetscBool id = PETSC_FALSE; 708 PetscCall(TSARKIMEXTestMassIdentity(ts, &id)); 709 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"); 710 } 711 PetscCall(TSClone(ts, &ts_start)); 712 PetscCall(TSSetSolution(ts_start, ts->vec_sol)); 713 PetscCall(TSSetTime(ts_start, ts->ptime)); 714 PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1)); 715 PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step)); 716 PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER)); 717 PetscCall(TSSetTimeStep(ts_start, ts->time_step)); 718 PetscCall(TSSetType(ts_start, TSARKIMEX)); 719 PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE)); 720 PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE)); 721 722 PetscCall(TSRestartStep(ts_start)); 723 PetscCall(TSSolve(ts_start, ts->vec_sol)); 724 PetscCall(TSGetTime(ts_start, &ts->ptime)); 725 PetscCall(TSGetTimeStep(ts_start, &ts->time_step)); 726 727 { /* Save the initial slope for the next step */ 728 TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data; 729 PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0)); 730 } 731 ts->steps++; 732 733 /* Set the correct TS in SNES */ 734 /* We'll try to bypass this by changing the method on the fly */ 735 { 736 PetscCall(TSGetSNES(ts, &snes)); 737 PetscCall(TSSetSNES(ts, snes)); 738 } 739 PetscCall(TSDestroy(&ts_start)); 740 } 741 742 ark->status = TS_STEP_INCOMPLETE; 743 while (!ts->reason && ark->status != TS_STEP_COMPLETE) { 744 PetscReal t = ts->ptime; 745 PetscReal h = ts->time_step; 746 for (i = 0; i < s; i++) { 747 ark->stage_time = t + h * ct[i]; 748 PetscCall(TSPreStage(ts, ark->stage_time)); 749 if (At[i * s + i] == 0) { /* This stage is explicit */ 750 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"); 751 PetscCall(VecCopy(ts->vec_sol, Y[i])); 752 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 753 PetscCall(VecMAXPY(Y[i], i, w, YdotI)); 754 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 755 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 756 } else { 757 ark->scoeff = 1. / At[i * s + i]; 758 /* Ydot = shift*(Y-Z) */ 759 PetscCall(VecCopy(ts->vec_sol, Z)); 760 for (j = 0; j < i; j++) w[j] = h * At[i * s + j]; 761 PetscCall(VecMAXPY(Z, i, w, YdotI)); 762 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 763 PetscCall(VecMAXPY(Z, i, w, YdotRHS)); 764 if (extrapolate && !ts->steprestart) { 765 /* Initial guess extrapolated from previous time step stage values */ 766 PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i])); 767 } else { 768 /* Initial guess taken from last stage */ 769 PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i])); 770 } 771 PetscCall(TSGetSNES(ts, &snes)); 772 PetscCall(SNESSolve(snes, NULL, Y[i])); 773 PetscCall(SNESGetIterationNumber(snes, &its)); 774 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 775 ts->snes_its += its; 776 ts->ksp_its += lits; 777 PetscCall(TSGetAdapt(ts, &adapt)); 778 PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok)); 779 if (!stageok) { 780 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to 781 * use extrapolation to initialize the solves on the next attempt. */ 782 extrapolate = PETSC_FALSE; 783 goto reject_step; 784 } 785 } 786 if (ts->equation_type >= TS_EQ_IMPLICIT) { 787 if (i == 0 && tab->explicit_first_stage) { 788 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", 789 ark->tableau->name); 790 PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */ 791 } else { 792 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 793 } 794 } else { 795 if (i == 0 && tab->explicit_first_stage) { 796 PetscCall(VecZeroEntries(Ydot)); 797 PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0) */ 798 PetscCall(VecScale(YdotI[i], -1.0)); 799 } else { 800 PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */ 801 } 802 if (ark->imex) { 803 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 804 } else { 805 PetscCall(VecZeroEntries(YdotRHS[i])); 806 } 807 } 808 PetscCall(TSPostStage(ts, ark->stage_time, i, Y)); 809 } 810 811 ark->status = TS_STEP_INCOMPLETE; 812 PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL)); 813 ark->status = TS_STEP_PENDING; 814 PetscCall(TSGetAdapt(ts, &adapt)); 815 PetscCall(TSAdaptCandidatesClear(adapt)); 816 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 817 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 818 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 819 if (!accept) { /* Roll back the current step */ 820 PetscCall(TSRollBack_ARKIMEX(ts)); 821 ts->time_step = next_time_step; 822 goto reject_step; 823 } 824 825 ts->ptime += ts->time_step; 826 ts->time_step = next_time_step; 827 break; 828 829 reject_step: 830 ts->reject++; 831 accept = PETSC_FALSE; 832 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 833 ts->reason = TS_DIVERGED_STEP_REJECTED; 834 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 835 } 836 } 837 PetscFunctionReturn(0); 838 } 839 840 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X) { 841 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 842 PetscInt s = ark->tableau->s, pinterp = ark->tableau->pinterp, i, j; 843 PetscReal h; 844 PetscReal tt, t; 845 PetscScalar *bt, *b; 846 const PetscReal *Bt = ark->tableau->binterpt, *B = ark->tableau->binterp; 847 848 PetscFunctionBegin; 849 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 850 switch (ark->status) { 851 case TS_STEP_INCOMPLETE: 852 case TS_STEP_PENDING: 853 h = ts->time_step; 854 t = (itime - ts->ptime) / h; 855 break; 856 case TS_STEP_COMPLETE: 857 h = ts->ptime - ts->ptime_prev; 858 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 859 break; 860 default: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 861 } 862 PetscCall(PetscMalloc2(s, &bt, s, &b)); 863 for (i = 0; i < s; i++) bt[i] = b[i] = 0; 864 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 865 for (i = 0; i < s; i++) { 866 bt[i] += h * Bt[i * pinterp + j] * tt; 867 b[i] += h * B[i * pinterp + j] * tt; 868 } 869 } 870 PetscCall(VecCopy(ark->Y[0], X)); 871 PetscCall(VecMAXPY(X, s, bt, ark->YdotI)); 872 PetscCall(VecMAXPY(X, s, b, ark->YdotRHS)); 873 PetscCall(PetscFree2(bt, b)); 874 PetscFunctionReturn(0); 875 } 876 877 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X) { 878 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 879 PetscInt s = ark->tableau->s, pinterp = ark->tableau->pinterp, i, j; 880 PetscReal h, h_prev, t, tt; 881 PetscScalar *bt, *b; 882 const PetscReal *Bt = ark->tableau->binterpt, *B = ark->tableau->binterp; 883 884 PetscFunctionBegin; 885 PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name); 886 PetscCall(PetscCalloc2(s, &bt, s, &b)); 887 h = ts->time_step; 888 h_prev = ts->ptime - ts->ptime_prev; 889 t = 1 + h / h_prev * c; 890 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 891 for (i = 0; i < s; i++) { 892 bt[i] += h * Bt[i * pinterp + j] * tt; 893 b[i] += h * B[i * pinterp + j] * tt; 894 } 895 } 896 PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored"); 897 PetscCall(VecCopy(ark->Y_prev[0], X)); 898 PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev)); 899 PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev)); 900 PetscCall(PetscFree2(bt, b)); 901 PetscFunctionReturn(0); 902 } 903 904 /*------------------------------------------------------------*/ 905 906 static PetscErrorCode TSARKIMEXTableauReset(TS ts) { 907 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 908 ARKTableau tab = ark->tableau; 909 910 PetscFunctionBegin; 911 if (!tab) PetscFunctionReturn(0); 912 PetscCall(PetscFree(ark->work)); 913 PetscCall(VecDestroyVecs(tab->s, &ark->Y)); 914 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI)); 915 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS)); 916 PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev)); 917 PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev)); 918 PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev)); 919 PetscFunctionReturn(0); 920 } 921 922 static PetscErrorCode TSReset_ARKIMEX(TS ts) { 923 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 924 925 PetscFunctionBegin; 926 PetscCall(TSARKIMEXTableauReset(ts)); 927 PetscCall(VecDestroy(&ark->Ydot)); 928 PetscCall(VecDestroy(&ark->Ydot0)); 929 PetscCall(VecDestroy(&ark->Z)); 930 PetscFunctionReturn(0); 931 } 932 933 static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) { 934 TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data; 935 936 PetscFunctionBegin; 937 if (Z) { 938 if (dm && dm != ts->dm) { 939 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); 940 } else *Z = ax->Z; 941 } 942 if (Ydot) { 943 if (dm && dm != ts->dm) { 944 PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); 945 } else *Ydot = ax->Ydot; 946 } 947 PetscFunctionReturn(0); 948 } 949 950 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot) { 951 PetscFunctionBegin; 952 if (Z) { 953 if (dm && dm != ts->dm) { PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z)); } 954 } 955 if (Ydot) { 956 if (dm && dm != ts->dm) { PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot)); } 957 } 958 PetscFunctionReturn(0); 959 } 960 961 /* 962 This defines the nonlinear equation that is to be solved with SNES 963 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 964 */ 965 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts) { 966 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 967 DM dm, dmsave; 968 Vec Z, Ydot; 969 PetscReal shift = ark->scoeff / ts->time_step; 970 971 PetscFunctionBegin; 972 PetscCall(SNESGetDM(snes, &dm)); 973 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot)); 974 PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */ 975 dmsave = ts->dm; 976 ts->dm = dm; 977 978 PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex)); 979 980 ts->dm = dmsave; 981 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot)); 982 PetscFunctionReturn(0); 983 } 984 985 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts) { 986 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 987 DM dm, dmsave; 988 Vec Ydot; 989 PetscReal shift = ark->scoeff / ts->time_step; 990 991 PetscFunctionBegin; 992 PetscCall(SNESGetDM(snes, &dm)); 993 PetscCall(TSARKIMEXGetVecs(ts, dm, NULL, &Ydot)); 994 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 995 dmsave = ts->dm; 996 ts->dm = dm; 997 998 PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex)); 999 1000 ts->dm = dmsave; 1001 PetscCall(TSARKIMEXRestoreVecs(ts, dm, NULL, &Ydot)); 1002 PetscFunctionReturn(0); 1003 } 1004 1005 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx) { 1006 PetscFunctionBegin; 1007 PetscFunctionReturn(0); 1008 } 1009 1010 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) { 1011 TS ts = (TS)ctx; 1012 Vec Z, Z_c; 1013 1014 PetscFunctionBegin; 1015 PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL)); 1016 PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL)); 1017 PetscCall(MatRestrict(restrct, Z, Z_c)); 1018 PetscCall(VecPointwiseMult(Z_c, rscale, Z_c)); 1019 PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL)); 1020 PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL)); 1021 PetscFunctionReturn(0); 1022 } 1023 1024 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx) { 1025 PetscFunctionBegin; 1026 PetscFunctionReturn(0); 1027 } 1028 1029 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) { 1030 TS ts = (TS)ctx; 1031 Vec Z, Z_c; 1032 1033 PetscFunctionBegin; 1034 PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL)); 1035 PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL)); 1036 1037 PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1038 PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD)); 1039 1040 PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL)); 1041 PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL)); 1042 PetscFunctionReturn(0); 1043 } 1044 1045 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts) { 1046 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1047 ARKTableau tab = ark->tableau; 1048 1049 PetscFunctionBegin; 1050 PetscCall(PetscMalloc1(tab->s, &ark->work)); 1051 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y)); 1052 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI)); 1053 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS)); 1054 if (ark->extrapolate) { 1055 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev)); 1056 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev)); 1057 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev)); 1058 } 1059 PetscFunctionReturn(0); 1060 } 1061 1062 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) { 1063 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1064 DM dm; 1065 SNES snes; 1066 1067 PetscFunctionBegin; 1068 PetscCall(TSARKIMEXTableauSetUp(ts)); 1069 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot)); 1070 PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0)); 1071 PetscCall(VecDuplicate(ts->vec_sol, &ark->Z)); 1072 PetscCall(TSGetDM(ts, &dm)); 1073 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 1074 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 1075 PetscCall(TSGetSNES(ts, &snes)); 1076 PetscFunctionReturn(0); 1077 } 1078 /*------------------------------------------------------------*/ 1079 1080 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject) { 1081 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1082 1083 PetscFunctionBegin; 1084 PetscOptionsHeadBegin(PetscOptionsObject, "ARKIMEX ODE solver options"); 1085 { 1086 ARKTableauLink link; 1087 PetscInt count, choice; 1088 PetscBool flg; 1089 const char **namelist; 1090 for (link = ARKTableauList, count = 0; link; link = link->next, count++) 1091 ; 1092 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1093 for (link = ARKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1094 PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg)); 1095 if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice])); 1096 PetscCall(PetscFree(namelist)); 1097 1098 flg = (PetscBool)!ark->imex; 1099 PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL)); 1100 ark->imex = (PetscBool)!flg; 1101 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)); 1102 } 1103 PetscOptionsHeadEnd(); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer) { 1108 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1109 PetscBool iascii; 1110 1111 PetscFunctionBegin; 1112 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1113 if (iascii) { 1114 ARKTableau tab = ark->tableau; 1115 TSARKIMEXType arktype; 1116 char buf[512]; 1117 PetscBool flg; 1118 1119 PetscCall(TSARKIMEXGetType(ts, &arktype)); 1120 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 1121 PetscCall(PetscViewerASCIIPrintf(viewer, " ARK IMEX %s\n", arktype)); 1122 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct)); 1123 PetscCall(PetscViewerASCIIPrintf(viewer, " Stiff abscissa ct = %s\n", buf)); 1124 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 1125 PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no")); 1126 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no")); 1127 PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no")); 1128 PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no")); 1129 PetscCall(PetscViewerASCIIPrintf(viewer, " Nonstiff abscissa c = %s\n", buf)); 1130 } 1131 PetscFunctionReturn(0); 1132 } 1133 1134 static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer) { 1135 SNES snes; 1136 TSAdapt adapt; 1137 1138 PetscFunctionBegin; 1139 PetscCall(TSGetAdapt(ts, &adapt)); 1140 PetscCall(TSAdaptLoad(adapt, viewer)); 1141 PetscCall(TSGetSNES(ts, &snes)); 1142 PetscCall(SNESLoad(snes, viewer)); 1143 /* function and Jacobian context for SNES when used with TS is always ts object */ 1144 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 1145 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 1146 PetscFunctionReturn(0); 1147 } 1148 1149 /*@C 1150 TSARKIMEXSetType - Set the type of ARK IMEX scheme 1151 1152 Logically collective 1153 1154 Input Parameters: 1155 + ts - timestepping context 1156 - arktype - type of ARK-IMEX scheme 1157 1158 Options Database: 1159 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set ARK IMEX scheme type 1160 1161 Level: intermediate 1162 1163 .seealso: `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`, 1164 `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5` 1165 @*/ 1166 PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype) { 1167 PetscFunctionBegin; 1168 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1169 PetscValidCharPointer(arktype, 2); 1170 PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype)); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 /*@C 1175 TSARKIMEXGetType - Get the type of ARK IMEX scheme 1176 1177 Logically collective 1178 1179 Input Parameter: 1180 . ts - timestepping context 1181 1182 Output Parameter: 1183 . arktype - type of ARK-IMEX scheme 1184 1185 Level: intermediate 1186 1187 .seealso: `TSARKIMEXGetType()` 1188 @*/ 1189 PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype) { 1190 PetscFunctionBegin; 1191 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1192 PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype)); 1193 PetscFunctionReturn(0); 1194 } 1195 1196 /*@ 1197 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 1198 1199 Logically collective 1200 1201 Input Parameters: 1202 + ts - timestepping context 1203 - flg - PETSC_TRUE for fully implicit 1204 1205 Level: intermediate 1206 1207 .seealso: `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()` 1208 @*/ 1209 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg) { 1210 PetscFunctionBegin; 1211 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1212 PetscValidLogicalCollectiveBool(ts, flg, 2); 1213 PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg)); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /*@ 1218 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly 1219 1220 Logically collective 1221 1222 Input Parameter: 1223 . ts - timestepping context 1224 1225 Output Parameter: 1226 . flg - PETSC_TRUE for fully implicit 1227 1228 Level: intermediate 1229 1230 .seealso: `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()` 1231 @*/ 1232 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg) { 1233 PetscFunctionBegin; 1234 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1235 PetscValidBoolPointer(flg, 2); 1236 PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg)); 1237 PetscFunctionReturn(0); 1238 } 1239 1240 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype) { 1241 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1242 1243 PetscFunctionBegin; 1244 *arktype = ark->tableau->name; 1245 PetscFunctionReturn(0); 1246 } 1247 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype) { 1248 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1249 PetscBool match; 1250 ARKTableauLink link; 1251 1252 PetscFunctionBegin; 1253 if (ark->tableau) { 1254 PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match)); 1255 if (match) PetscFunctionReturn(0); 1256 } 1257 for (link = ARKTableauList; link; link = link->next) { 1258 PetscCall(PetscStrcmp(link->tab.name, arktype, &match)); 1259 if (match) { 1260 if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts)); 1261 ark->tableau = &link->tab; 1262 if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts)); 1263 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1264 PetscFunctionReturn(0); 1265 } 1266 } 1267 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype); 1268 } 1269 1270 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg) { 1271 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1272 1273 PetscFunctionBegin; 1274 ark->imex = (PetscBool)!flg; 1275 PetscFunctionReturn(0); 1276 } 1277 1278 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg) { 1279 TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data; 1280 1281 PetscFunctionBegin; 1282 *flg = (PetscBool)!ark->imex; 1283 PetscFunctionReturn(0); 1284 } 1285 1286 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) { 1287 PetscFunctionBegin; 1288 PetscCall(TSReset_ARKIMEX(ts)); 1289 if (ts->dm) { 1290 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts)); 1291 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts)); 1292 } 1293 PetscCall(PetscFree(ts->data)); 1294 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL)); 1297 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL)); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /* ------------------------------------------------------------ */ 1302 /*MC 1303 TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes 1304 1305 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1306 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1307 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1308 1309 Notes: 1310 The default is TSARKIMEX3, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 1311 1312 If the equation is implicit or a DAE, then TSSetEquationType() needs to be set accordingly. Refer to the manual for further information. 1313 1314 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). 1315 1316 Consider trying TSROSW if the stiff part is linear or weakly nonlinear. 1317 1318 Level: beginner 1319 1320 .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`, 1321 `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`, 1322 `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()` 1323 1324 M*/ 1325 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) { 1326 TS_ARKIMEX *th; 1327 1328 PetscFunctionBegin; 1329 PetscCall(TSARKIMEXInitializePackage()); 1330 1331 ts->ops->reset = TSReset_ARKIMEX; 1332 ts->ops->destroy = TSDestroy_ARKIMEX; 1333 ts->ops->view = TSView_ARKIMEX; 1334 ts->ops->load = TSLoad_ARKIMEX; 1335 ts->ops->setup = TSSetUp_ARKIMEX; 1336 ts->ops->step = TSStep_ARKIMEX; 1337 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1338 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1339 ts->ops->rollback = TSRollBack_ARKIMEX; 1340 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1341 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1342 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1343 1344 ts->usessnes = PETSC_TRUE; 1345 1346 PetscCall(PetscNewLog(ts, &th)); 1347 ts->data = (void *)th; 1348 th->imex = PETSC_TRUE; 1349 1350 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX)); 1351 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX)); 1352 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX)); 1353 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX)); 1354 1355 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault)); 1356 PetscFunctionReturn(0); 1357 } 1358