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