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