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 ark->status = TS_STEP_INCOMPLETE; 839 ierr = TSRollBack(ts);CHKERRQ(ierr); 840 } 841 reject_step: continue; 842 } 843 if (ark->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "TSInterpolate_ARKIMEX" 849 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X) 850 { 851 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 852 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 853 PetscReal h; 854 PetscReal tt,t; 855 PetscScalar *bt,*b; 856 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 857 PetscErrorCode ierr; 858 859 PetscFunctionBegin; 860 if (!Bt || !B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name); 861 switch (ark->status) { 862 case TS_STEP_INCOMPLETE: 863 case TS_STEP_PENDING: 864 h = ts->time_step; 865 t = (itime - ts->ptime)/h; 866 break; 867 case TS_STEP_COMPLETE: 868 h = ts->time_step_prev; 869 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 870 break; 871 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 872 } 873 ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 874 for (i=0; i<s; i++) bt[i] = b[i] = 0; 875 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 876 for (i=0; i<s; i++) { 877 bt[i] += h * Bt[i*pinterp+j] * tt * -1.0; 878 b[i] += h * B[i*pinterp+j] * tt; 879 } 880 } 881 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"); 882 ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 883 ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 884 ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 885 ierr = PetscFree2(bt,b);CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 /*------------------------------------------------------------*/ 890 #undef __FUNCT__ 891 #define __FUNCT__ "TSReset_ARKIMEX" 892 static PetscErrorCode TSReset_ARKIMEX(TS ts) 893 { 894 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 895 PetscInt s; 896 PetscErrorCode ierr; 897 898 PetscFunctionBegin; 899 if (!ark->tableau) PetscFunctionReturn(0); 900 s = ark->tableau->s; 901 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 902 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 903 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 904 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 905 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 906 ierr = VecDestroy(&ark->Ydot0);CHKERRQ(ierr); 907 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 908 ierr = PetscFree(ark->work);CHKERRQ(ierr); 909 PetscFunctionReturn(0); 910 } 911 912 #undef __FUNCT__ 913 #define __FUNCT__ "TSDestroy_ARKIMEX" 914 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 915 { 916 PetscErrorCode ierr; 917 918 PetscFunctionBegin; 919 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 920 ierr = PetscFree(ts->data);CHKERRQ(ierr); 921 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetType_C",NULL);CHKERRQ(ierr); 922 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetType_C",NULL);CHKERRQ(ierr); 923 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",NULL);CHKERRQ(ierr); 924 PetscFunctionReturn(0); 925 } 926 927 928 #undef __FUNCT__ 929 #define __FUNCT__ "TSARKIMEXGetVecs" 930 static PetscErrorCode TSARKIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot) 931 { 932 TS_ARKIMEX *ax = (TS_ARKIMEX*)ts->data; 933 PetscErrorCode ierr; 934 935 PetscFunctionBegin; 936 if (Z) { 937 if (dm && dm != ts->dm) { 938 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr); 939 } else *Z = ax->Z; 940 } 941 if (Ydot) { 942 if (dm && dm != ts->dm) { 943 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr); 944 } else *Ydot = ax->Ydot; 945 } 946 PetscFunctionReturn(0); 947 } 948 949 950 #undef __FUNCT__ 951 #define __FUNCT__ "TSARKIMEXRestoreVecs" 952 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot) 953 { 954 PetscErrorCode ierr; 955 956 PetscFunctionBegin; 957 if (Z) { 958 if (dm && dm != ts->dm) { 959 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr); 960 } 961 } 962 if (Ydot) { 963 if (dm && dm != ts->dm) { 964 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr); 965 } 966 } 967 PetscFunctionReturn(0); 968 } 969 970 /* 971 This defines the nonlinear equation that is to be solved with SNES 972 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 973 */ 974 #undef __FUNCT__ 975 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 976 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 977 { 978 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 979 DM dm,dmsave; 980 Vec Z,Ydot; 981 PetscReal shift = ark->scoeff / ts->time_step; 982 PetscErrorCode ierr; 983 984 PetscFunctionBegin; 985 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 986 ierr = TSARKIMEXGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 987 ierr = VecAXPBYPCZ(Ydot,-shift,shift,0,Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 988 dmsave = ts->dm; 989 ts->dm = dm; 990 991 ierr = TSComputeIFunction(ts,ark->stage_time,X,Ydot,F,ark->imex);CHKERRQ(ierr); 992 993 ts->dm = dmsave; 994 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 995 PetscFunctionReturn(0); 996 } 997 998 #undef __FUNCT__ 999 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 1000 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 1001 { 1002 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1003 DM dm,dmsave; 1004 Vec Ydot; 1005 PetscReal shift = ark->scoeff / ts->time_step; 1006 PetscErrorCode ierr; 1007 1008 PetscFunctionBegin; 1009 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1010 ierr = TSARKIMEXGetVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); 1011 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 1012 dmsave = ts->dm; 1013 ts->dm = dm; 1014 1015 ierr = TSComputeIJacobian(ts,ark->stage_time,X,Ydot,shift,A,B,str,ark->imex);CHKERRQ(ierr); 1016 1017 ts->dm = dmsave; 1018 ierr = TSARKIMEXRestoreVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); 1019 PetscFunctionReturn(0); 1020 } 1021 1022 #undef __FUNCT__ 1023 #define __FUNCT__ "DMCoarsenHook_TSARKIMEX" 1024 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine,DM coarse,void *ctx) 1025 { 1026 PetscFunctionBegin; 1027 PetscFunctionReturn(0); 1028 } 1029 1030 #undef __FUNCT__ 1031 #define __FUNCT__ "DMRestrictHook_TSARKIMEX" 1032 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1033 { 1034 TS ts = (TS)ctx; 1035 PetscErrorCode ierr; 1036 Vec Z,Z_c; 1037 1038 PetscFunctionBegin; 1039 ierr = TSARKIMEXGetVecs(ts,fine,&Z,NULL);CHKERRQ(ierr); 1040 ierr = TSARKIMEXGetVecs(ts,coarse,&Z_c,NULL);CHKERRQ(ierr); 1041 ierr = MatRestrict(restrct,Z,Z_c);CHKERRQ(ierr); 1042 ierr = VecPointwiseMult(Z_c,rscale,Z_c);CHKERRQ(ierr); 1043 ierr = TSARKIMEXRestoreVecs(ts,fine,&Z,NULL);CHKERRQ(ierr); 1044 ierr = TSARKIMEXRestoreVecs(ts,coarse,&Z_c,NULL);CHKERRQ(ierr); 1045 PetscFunctionReturn(0); 1046 } 1047 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "DMSubDomainHook_TSARKIMEX" 1051 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm,DM subdm,void *ctx) 1052 { 1053 PetscFunctionBegin; 1054 PetscFunctionReturn(0); 1055 } 1056 1057 #undef __FUNCT__ 1058 #define __FUNCT__ "DMSubDomainRestrictHook_TSARKIMEX" 1059 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1060 { 1061 TS ts = (TS)ctx; 1062 PetscErrorCode ierr; 1063 Vec Z,Z_c; 1064 1065 PetscFunctionBegin; 1066 ierr = TSARKIMEXGetVecs(ts,dm,&Z,NULL);CHKERRQ(ierr); 1067 ierr = TSARKIMEXGetVecs(ts,subdm,&Z_c,NULL);CHKERRQ(ierr); 1068 1069 ierr = VecScatterBegin(gscat,Z,Z_c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1070 ierr = VecScatterEnd(gscat,Z,Z_c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1071 1072 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,NULL);CHKERRQ(ierr); 1073 ierr = TSARKIMEXRestoreVecs(ts,subdm,&Z_c,NULL);CHKERRQ(ierr); 1074 PetscFunctionReturn(0); 1075 } 1076 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "TSSetUp_ARKIMEX" 1079 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 1080 { 1081 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1082 ARKTableau tab; 1083 PetscInt s; 1084 PetscErrorCode ierr; 1085 DM dm; 1086 1087 PetscFunctionBegin; 1088 if (!ark->tableau) { 1089 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 1090 } 1091 tab = ark->tableau; 1092 s = tab->s; 1093 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 1094 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 1095 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 1096 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 1097 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 1098 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot0);CHKERRQ(ierr); 1099 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 1100 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 1101 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1102 if (dm) { 1103 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSARKIMEX,DMRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr); 1104 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSARKIMEX,DMSubDomainRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr); 1105 } 1106 PetscFunctionReturn(0); 1107 } 1108 /*------------------------------------------------------------*/ 1109 1110 #undef __FUNCT__ 1111 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 1112 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 1113 { 1114 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1115 PetscErrorCode ierr; 1116 char arktype[256]; 1117 1118 PetscFunctionBegin; 1119 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 1120 { 1121 ARKTableauLink link; 1122 PetscInt count,choice; 1123 PetscBool flg; 1124 const char **namelist; 1125 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof(arktype));CHKERRQ(ierr); 1126 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 1127 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 1128 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1129 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 1130 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 1131 ierr = PetscFree(namelist);CHKERRQ(ierr); 1132 flg = (PetscBool) !ark->imex; 1133 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,NULL);CHKERRQ(ierr); 1134 ark->imex = (PetscBool) !flg; 1135 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 1136 } 1137 ierr = PetscOptionsTail();CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 #undef __FUNCT__ 1142 #define __FUNCT__ "PetscFormatRealArray" 1143 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1144 { 1145 PetscErrorCode ierr; 1146 PetscInt i; 1147 size_t left,count; 1148 char *p; 1149 1150 PetscFunctionBegin; 1151 for (i=0,p=buf,left=len; i<n; i++) { 1152 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1153 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1154 left -= count; 1155 p += count; 1156 *p++ = ' '; 1157 } 1158 p[i ? 0 : -1] = 0; 1159 PetscFunctionReturn(0); 1160 } 1161 1162 #undef __FUNCT__ 1163 #define __FUNCT__ "TSView_ARKIMEX" 1164 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 1165 { 1166 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1167 ARKTableau tab = ark->tableau; 1168 PetscBool iascii; 1169 PetscErrorCode ierr; 1170 TSAdapt adapt; 1171 1172 PetscFunctionBegin; 1173 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1174 if (iascii) { 1175 TSARKIMEXType arktype; 1176 char buf[512]; 1177 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 1178 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 1179 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 1180 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 1181 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1182 ierr = PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s\n",tab->stiffly_accurate ? "yes" : "no");CHKERRQ(ierr); 1183 ierr = PetscViewerASCIIPrintf(viewer,"Explicit first stage: %s\n",tab->explicit_first_stage ? "yes" : "no");CHKERRQ(ierr); 1184 ierr = PetscViewerASCIIPrintf(viewer,"FSAL property: %s\n",tab->FSAL_implicit ? "yes" : "no");CHKERRQ(ierr); 1185 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 1186 } 1187 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1188 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1189 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 #undef __FUNCT__ 1194 #define __FUNCT__ "TSLoad_ARKIMEX" 1195 static PetscErrorCode TSLoad_ARKIMEX(TS ts,PetscViewer viewer) 1196 { 1197 PetscErrorCode ierr; 1198 SNES snes; 1199 TSAdapt tsadapt; 1200 1201 PetscFunctionBegin; 1202 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 1203 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 1204 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1205 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 1206 /* function and Jacobian context for SNES when used with TS is always ts object */ 1207 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 1208 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 1209 PetscFunctionReturn(0); 1210 } 1211 1212 #undef __FUNCT__ 1213 #define __FUNCT__ "TSARKIMEXSetType" 1214 /*@C 1215 TSARKIMEXSetType - Set the type of ARK IMEX scheme 1216 1217 Logically collective 1218 1219 Input Parameter: 1220 + ts - timestepping context 1221 - arktype - type of ARK-IMEX scheme 1222 1223 Level: intermediate 1224 1225 .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 1226 @*/ 1227 PetscErrorCode TSARKIMEXSetType(TS ts,TSARKIMEXType arktype) 1228 { 1229 PetscErrorCode ierr; 1230 1231 PetscFunctionBegin; 1232 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1233 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 1234 PetscFunctionReturn(0); 1235 } 1236 1237 #undef __FUNCT__ 1238 #define __FUNCT__ "TSARKIMEXGetType" 1239 /*@C 1240 TSARKIMEXGetType - Get the type of ARK IMEX scheme 1241 1242 Logically collective 1243 1244 Input Parameter: 1245 . ts - timestepping context 1246 1247 Output Parameter: 1248 . arktype - type of ARK-IMEX scheme 1249 1250 Level: intermediate 1251 1252 .seealso: TSARKIMEXGetType() 1253 @*/ 1254 PetscErrorCode TSARKIMEXGetType(TS ts,TSARKIMEXType *arktype) 1255 { 1256 PetscErrorCode ierr; 1257 1258 PetscFunctionBegin; 1259 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1260 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 #undef __FUNCT__ 1265 #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 1266 /*@C 1267 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 1268 1269 Logically collective 1270 1271 Input Parameter: 1272 + ts - timestepping context 1273 - flg - PETSC_TRUE for fully implicit 1274 1275 Level: intermediate 1276 1277 .seealso: TSARKIMEXGetType() 1278 @*/ 1279 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 1280 { 1281 PetscErrorCode ierr; 1282 1283 PetscFunctionBegin; 1284 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1285 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1286 PetscFunctionReturn(0); 1287 } 1288 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 1291 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType *arktype) 1292 { 1293 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1294 PetscErrorCode ierr; 1295 1296 PetscFunctionBegin; 1297 if (!ark->tableau) { 1298 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 1299 } 1300 *arktype = ark->tableau->name; 1301 PetscFunctionReturn(0); 1302 } 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 1305 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype) 1306 { 1307 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1308 PetscErrorCode ierr; 1309 PetscBool match; 1310 ARKTableauLink link; 1311 1312 PetscFunctionBegin; 1313 if (ark->tableau) { 1314 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 1315 if (match) PetscFunctionReturn(0); 1316 } 1317 for (link = ARKTableauList; link; link=link->next) { 1318 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 1319 if (match) { 1320 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 1321 ark->tableau = &link->tab; 1322 PetscFunctionReturn(0); 1323 } 1324 } 1325 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 1326 PetscFunctionReturn(0); 1327 } 1328 #undef __FUNCT__ 1329 #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 1330 PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 1331 { 1332 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1333 1334 PetscFunctionBegin; 1335 ark->imex = (PetscBool)!flg; 1336 PetscFunctionReturn(0); 1337 } 1338 1339 /* ------------------------------------------------------------ */ 1340 /*MC 1341 TSARKIMEX - ODE and DAE solver using Additive Runge-Kutta IMEX schemes 1342 1343 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1344 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1345 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1346 1347 Notes: 1348 The default is TSARKIMEX3, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 1349 1350 Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 1351 1352 Level: beginner 1353 1354 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 1355 TSARKIMEX4, TSARKIMEX5, TSARKIMEXPRSSP2, TSARKIMEXBPR3, TSARKIMEXType, TSARKIMEXRegister() 1356 1357 M*/ 1358 #undef __FUNCT__ 1359 #define __FUNCT__ "TSCreate_ARKIMEX" 1360 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 1361 { 1362 TS_ARKIMEX *th; 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1367 ierr = TSARKIMEXInitializePackage();CHKERRQ(ierr); 1368 #endif 1369 1370 ts->ops->reset = TSReset_ARKIMEX; 1371 ts->ops->destroy = TSDestroy_ARKIMEX; 1372 ts->ops->view = TSView_ARKIMEX; 1373 ts->ops->load = TSLoad_ARKIMEX; 1374 ts->ops->setup = TSSetUp_ARKIMEX; 1375 ts->ops->step = TSStep_ARKIMEX; 1376 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1377 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1378 ts->ops->rollback = TSRollBack_ARKIMEX; 1379 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1380 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1381 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1382 1383 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 1384 ts->data = (void*)th; 1385 th->imex = PETSC_TRUE; 1386 1387 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetType_C",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 1388 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetType_C",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 1389 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 1390 PetscFunctionReturn(0); 1391 } 1392