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