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