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__ "TSStep_ARKIMEX" 665 static PetscErrorCode TSStep_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 *At = tab->At,*A = tab->A,*bt = tab->bt,*b = tab->b,*ct = tab->ct,*c = tab->c; 671 PetscScalar *w = ark->work; 672 Vec *Y = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,Ydot0 = ark->Ydot0,W = ark->Work,Z = ark->Z; 673 TSAdapt adapt; 674 SNES snes; 675 PetscInt i,j,its,lits,reject,next_scheme; 676 PetscReal next_time_step; 677 PetscReal t; 678 PetscBool accept; 679 PetscErrorCode ierr; 680 681 PetscFunctionBegin; 682 if (ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage) { 683 PetscReal valid_time; 684 PetscBool isvalid; 685 ierr = PetscObjectComposedDataGetReal((PetscObject)ts->vec_sol, 686 explicit_stage_time_id, 687 valid_time, 688 isvalid); 689 CHKERRQ(ierr); 690 if (!isvalid || valid_time != ts->ptime) { 691 TS ts_start; 692 SNES snes_start; 693 DM dm; 694 PetscReal atol; 695 Vec vatol; 696 PetscReal rtol; 697 Vec vrtol; 698 699 ierr = TSCreate(PETSC_COMM_WORLD,&ts_start);CHKERRQ(ierr); 700 ierr = TSGetSNES(ts,&snes_start);CHKERRQ(ierr); 701 ierr = TSSetSNES(ts_start,snes_start);CHKERRQ(ierr); 702 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 703 ierr = TSSetDM(ts_start,dm);CHKERRQ(ierr); 704 705 ts_start->adapt=ts->adapt; 706 PetscObjectReference((PetscObject)ts_start->adapt); 707 708 ierr = TSSetSolution(ts_start,ts->vec_sol);CHKERRQ(ierr); 709 ierr = TSSetTime(ts_start,ts->ptime);CHKERRQ(ierr); 710 ierr = TSSetDuration(ts_start,1,ts->time_step);CHKERRQ(ierr); 711 ierr = TSSetTimeStep(ts_start,ts->time_step);CHKERRQ(ierr); 712 ierr = TSSetType(ts_start,TSARKIMEX);CHKERRQ(ierr); 713 ierr = TSARKIMEXSetFullyImplicit(ts_start,PETSC_TRUE);CHKERRQ(ierr); 714 ierr = TSARKIMEXSetType(ts_start,TSARKIMEX1BEE);CHKERRQ(ierr); 715 ierr = TSSetEquationType(ts_start,ts->equation_type);CHKERRQ(ierr); 716 ierr = TSGetTolerances(ts,&atol,&vatol,&rtol,&vrtol);CHKERRQ(ierr); 717 ierr = TSSetTolerances(ts_start,atol,vatol,rtol,vrtol);CHKERRQ(ierr); 718 ierr = TSSolve(ts_start,ts->vec_sol);CHKERRQ(ierr); 719 ierr = TSGetTime(ts_start,&ts->ptime);CHKERRQ(ierr); 720 721 ts->time_step = ts_start->time_step; 722 ts->steps++; 723 ierr = VecCopy(((TS_ARKIMEX*)ts_start->data)->Ydot0,Ydot0);CHKERRQ(ierr); 724 ierr = TSDestroy(&ts_start);CHKERRQ(ierr); 725 ierr = TSSetSNES(ts,snes_start);CHKERRQ(ierr); 726 } 727 } 728 729 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 730 next_time_step = ts->time_step; 731 t = ts->ptime; 732 accept = PETSC_TRUE; 733 ark->status = TS_STEP_INCOMPLETE; 734 735 736 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 737 PetscReal h = ts->time_step; 738 ierr = TSPreStep(ts);CHKERRQ(ierr); 739 for (i=0; i<s; i++) { 740 if (At[i*s+i] == 0) { /* This stage is explicit */ 741 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 742 for (j=0; j<i; j++) w[j] = h*At[i*s+j]; 743 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 744 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 745 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 746 } else { 747 ark->stage_time = t + h*ct[i]; 748 ark->scoeff = 1./At[i*s+i]; 749 ierr = TSPreStage(ts,ark->stage_time);CHKERRQ(ierr); 750 /* Affine part */ 751 ierr = VecZeroEntries(W);CHKERRQ(ierr); 752 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 753 ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 754 ierr = VecScale(W, ark->scoeff/h);CHKERRQ(ierr); 755 756 /* Ydot = shift*(Y-Z) */ 757 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 758 for (j=0; j<i; j++) w[j] = h*At[i*s+j]; 759 ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr); 760 761 /* Initial guess taken from last stage */ 762 ierr = VecCopy(i>0 ? Y[i-1] : ts->vec_sol,Y[i]);CHKERRQ(ierr); 763 ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 764 ierr = (ts->ops->snesfunction)(snes,Y[i],W,ts);CHKERRQ(ierr); 765 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 766 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 767 ts->snes_its += its; ts->ksp_its += lits; 768 ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 769 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 770 if (!accept) goto reject_step; 771 } 772 if (ts->equation_type>=TS_EQ_IMPLICIT) { 773 if (i==0 && tab->explicit_first_stage) { 774 ierr = VecCopy(Ydot0,YdotI[0]);CHKERRQ(ierr); 775 } else { 776 ierr = VecAXPBYPCZ(YdotI[i],-ark->scoeff/h,ark->scoeff/h,0,Z,Y[i]);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 777 } 778 } else { 779 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 780 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr); 781 ierr = VecScale(YdotI[i], -1.0);CHKERRQ(ierr); 782 if (ark->imex) { 783 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 784 } else { 785 ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr); 786 } 787 } 788 } 789 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 790 ark->status = TS_STEP_PENDING; 791 792 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 793 ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 794 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 795 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 796 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 797 if (accept) { 798 /* ignore next_scheme for now */ 799 ts->ptime += ts->time_step; 800 ts->time_step = next_time_step; 801 ts->steps++; 802 if (ts->equation_type>=TS_EQ_IMPLICIT) { /* save the initial slope for the next step*/ 803 ierr = VecCopy(YdotI[s-1],Ydot0);CHKERRQ(ierr); 804 } 805 ark->status = TS_STEP_COMPLETE; 806 if (tab->explicit_first_stage) { 807 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 808 } 809 810 break; 811 } else { /* Roll back the current step */ 812 for (j=0; j<s; j++) w[j] = -h*bt[j]; 813 ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotI);CHKERRQ(ierr); 814 for (j=0; j<s; j++) w[j] = -h*b[j]; 815 ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotRHS);CHKERRQ(ierr); 816 ts->time_step = next_time_step; 817 ark->status = TS_STEP_INCOMPLETE; 818 } 819 reject_step: continue; 820 } 821 if (ark->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 822 PetscFunctionReturn(0); 823 } 824 825 #undef __FUNCT__ 826 #define __FUNCT__ "TSInterpolate_ARKIMEX" 827 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X) 828 { 829 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 830 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 831 PetscReal h; 832 PetscReal tt,t; 833 PetscScalar *bt,*b; 834 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 835 PetscErrorCode ierr; 836 837 PetscFunctionBegin; 838 if (!Bt || !B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name); 839 switch (ark->status) { 840 case TS_STEP_INCOMPLETE: 841 case TS_STEP_PENDING: 842 h = ts->time_step; 843 t = (itime - ts->ptime)/h; 844 break; 845 case TS_STEP_COMPLETE: 846 h = ts->time_step_prev; 847 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 848 break; 849 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 850 } 851 ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 852 for (i=0; i<s; i++) bt[i] = b[i] = 0; 853 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 854 for (i=0; i<s; i++) { 855 bt[i] += h * Bt[i*pinterp+j] * tt * -1.0; 856 b[i] += h * B[i*pinterp+j] * tt; 857 } 858 } 859 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"); 860 ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 861 ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 862 ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 863 ierr = PetscFree2(bt,b);CHKERRQ(ierr); 864 PetscFunctionReturn(0); 865 } 866 867 /*------------------------------------------------------------*/ 868 #undef __FUNCT__ 869 #define __FUNCT__ "TSReset_ARKIMEX" 870 static PetscErrorCode TSReset_ARKIMEX(TS ts) 871 { 872 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 873 PetscInt s; 874 PetscErrorCode ierr; 875 876 PetscFunctionBegin; 877 if (!ark->tableau) PetscFunctionReturn(0); 878 s = ark->tableau->s; 879 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 880 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 881 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 882 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 883 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 884 ierr = VecDestroy(&ark->Ydot0);CHKERRQ(ierr); 885 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 886 ierr = PetscFree(ark->work);CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 #undef __FUNCT__ 891 #define __FUNCT__ "TSDestroy_ARKIMEX" 892 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 893 { 894 PetscErrorCode ierr; 895 896 PetscFunctionBegin; 897 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 898 ierr = PetscFree(ts->data);CHKERRQ(ierr); 899 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetType_C",NULL);CHKERRQ(ierr); 900 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetType_C",NULL);CHKERRQ(ierr); 901 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",NULL);CHKERRQ(ierr); 902 PetscFunctionReturn(0); 903 } 904 905 906 #undef __FUNCT__ 907 #define __FUNCT__ "TSARKIMEXGetVecs" 908 static PetscErrorCode TSARKIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot) 909 { 910 TS_ARKIMEX *ax = (TS_ARKIMEX*)ts->data; 911 PetscErrorCode ierr; 912 913 PetscFunctionBegin; 914 if (Z) { 915 if (dm && dm != ts->dm) { 916 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr); 917 } else *Z = ax->Z; 918 } 919 if (Ydot) { 920 if (dm && dm != ts->dm) { 921 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr); 922 } else *Ydot = ax->Ydot; 923 } 924 PetscFunctionReturn(0); 925 } 926 927 928 #undef __FUNCT__ 929 #define __FUNCT__ "TSARKIMEXRestoreVecs" 930 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot) 931 { 932 PetscErrorCode ierr; 933 934 PetscFunctionBegin; 935 if (Z) { 936 if (dm && dm != ts->dm) { 937 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr); 938 } 939 } 940 if (Ydot) { 941 if (dm && dm != ts->dm) { 942 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr); 943 } 944 } 945 PetscFunctionReturn(0); 946 } 947 948 /* 949 This defines the nonlinear equation that is to be solved with SNES 950 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 951 */ 952 #undef __FUNCT__ 953 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 954 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 955 { 956 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 957 DM dm,dmsave; 958 Vec Z,Ydot; 959 PetscReal shift = ark->scoeff / ts->time_step; 960 PetscErrorCode ierr; 961 962 PetscFunctionBegin; 963 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 964 ierr = TSARKIMEXGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 965 ierr = VecAXPBYPCZ(Ydot,-shift,shift,0,Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 966 dmsave = ts->dm; 967 ts->dm = dm; 968 969 ierr = TSComputeIFunction(ts,ark->stage_time,X,Ydot,F,ark->imex);CHKERRQ(ierr); 970 971 ts->dm = dmsave; 972 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 973 PetscFunctionReturn(0); 974 } 975 976 #undef __FUNCT__ 977 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 978 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 979 { 980 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 981 DM dm,dmsave; 982 Vec Ydot; 983 PetscReal shift = ark->scoeff / ts->time_step; 984 PetscErrorCode ierr; 985 986 PetscFunctionBegin; 987 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 988 ierr = TSARKIMEXGetVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); 989 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 990 dmsave = ts->dm; 991 ts->dm = dm; 992 993 ierr = TSComputeIJacobian(ts,ark->stage_time,X,Ydot,shift,A,B,str,ark->imex);CHKERRQ(ierr); 994 995 ts->dm = dmsave; 996 ierr = TSARKIMEXRestoreVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); 997 PetscFunctionReturn(0); 998 } 999 1000 #undef __FUNCT__ 1001 #define __FUNCT__ "DMCoarsenHook_TSARKIMEX" 1002 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine,DM coarse,void *ctx) 1003 { 1004 PetscFunctionBegin; 1005 PetscFunctionReturn(0); 1006 } 1007 1008 #undef __FUNCT__ 1009 #define __FUNCT__ "DMRestrictHook_TSARKIMEX" 1010 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1011 { 1012 TS ts = (TS)ctx; 1013 PetscErrorCode ierr; 1014 Vec Z,Z_c; 1015 1016 PetscFunctionBegin; 1017 ierr = TSARKIMEXGetVecs(ts,fine,&Z,NULL);CHKERRQ(ierr); 1018 ierr = TSARKIMEXGetVecs(ts,coarse,&Z_c,NULL);CHKERRQ(ierr); 1019 ierr = MatRestrict(restrct,Z,Z_c);CHKERRQ(ierr); 1020 ierr = VecPointwiseMult(Z_c,rscale,Z_c);CHKERRQ(ierr); 1021 ierr = TSARKIMEXRestoreVecs(ts,fine,&Z,NULL);CHKERRQ(ierr); 1022 ierr = TSARKIMEXRestoreVecs(ts,coarse,&Z_c,NULL);CHKERRQ(ierr); 1023 PetscFunctionReturn(0); 1024 } 1025 1026 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "DMSubDomainHook_TSARKIMEX" 1029 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm,DM subdm,void *ctx) 1030 { 1031 PetscFunctionBegin; 1032 PetscFunctionReturn(0); 1033 } 1034 1035 #undef __FUNCT__ 1036 #define __FUNCT__ "DMSubDomainRestrictHook_TSARKIMEX" 1037 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1038 { 1039 TS ts = (TS)ctx; 1040 PetscErrorCode ierr; 1041 Vec Z,Z_c; 1042 1043 PetscFunctionBegin; 1044 ierr = TSARKIMEXGetVecs(ts,dm,&Z,NULL);CHKERRQ(ierr); 1045 ierr = TSARKIMEXGetVecs(ts,subdm,&Z_c,NULL);CHKERRQ(ierr); 1046 1047 ierr = VecScatterBegin(gscat,Z,Z_c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1048 ierr = VecScatterEnd(gscat,Z,Z_c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1049 1050 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,NULL);CHKERRQ(ierr); 1051 ierr = TSARKIMEXRestoreVecs(ts,subdm,&Z_c,NULL);CHKERRQ(ierr); 1052 PetscFunctionReturn(0); 1053 } 1054 1055 #undef __FUNCT__ 1056 #define __FUNCT__ "TSSetUp_ARKIMEX" 1057 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 1058 { 1059 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1060 ARKTableau tab; 1061 PetscInt s; 1062 PetscErrorCode ierr; 1063 DM dm; 1064 1065 PetscFunctionBegin; 1066 if (!ark->tableau) { 1067 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 1068 } 1069 tab = ark->tableau; 1070 s = tab->s; 1071 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 1072 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 1073 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 1074 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 1075 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 1076 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot0);CHKERRQ(ierr); 1077 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 1078 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 1079 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1080 if (dm) { 1081 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSARKIMEX,DMRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr); 1082 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSARKIMEX,DMSubDomainRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr); 1083 } 1084 PetscFunctionReturn(0); 1085 } 1086 /*------------------------------------------------------------*/ 1087 1088 #undef __FUNCT__ 1089 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 1090 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 1091 { 1092 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1093 PetscErrorCode ierr; 1094 char arktype[256]; 1095 1096 PetscFunctionBegin; 1097 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 1098 { 1099 ARKTableauLink link; 1100 PetscInt count,choice; 1101 PetscBool flg; 1102 const char **namelist; 1103 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof(arktype));CHKERRQ(ierr); 1104 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 1105 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 1106 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1107 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 1108 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 1109 ierr = PetscFree(namelist);CHKERRQ(ierr); 1110 flg = (PetscBool) !ark->imex; 1111 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,NULL);CHKERRQ(ierr); 1112 ark->imex = (PetscBool) !flg; 1113 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 1114 } 1115 ierr = PetscOptionsTail();CHKERRQ(ierr); 1116 PetscFunctionReturn(0); 1117 } 1118 1119 #undef __FUNCT__ 1120 #define __FUNCT__ "PetscFormatRealArray" 1121 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1122 { 1123 PetscErrorCode ierr; 1124 PetscInt i; 1125 size_t left,count; 1126 char *p; 1127 1128 PetscFunctionBegin; 1129 for (i=0,p=buf,left=len; i<n; i++) { 1130 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1131 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1132 left -= count; 1133 p += count; 1134 *p++ = ' '; 1135 } 1136 p[i ? 0 : -1] = 0; 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "TSView_ARKIMEX" 1142 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 1143 { 1144 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1145 ARKTableau tab = ark->tableau; 1146 PetscBool iascii; 1147 PetscErrorCode ierr; 1148 TSAdapt adapt; 1149 1150 PetscFunctionBegin; 1151 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1152 if (iascii) { 1153 TSARKIMEXType arktype; 1154 char buf[512]; 1155 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 1156 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 1157 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 1158 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 1159 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1160 ierr = PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s\n",tab->stiffly_accurate ? "yes" : "no");CHKERRQ(ierr); 1161 ierr = PetscViewerASCIIPrintf(viewer,"Explicit first stage: %s\n",tab->explicit_first_stage ? "yes" : "no");CHKERRQ(ierr); 1162 ierr = PetscViewerASCIIPrintf(viewer,"FSAL property: %s\n",tab->FSAL_implicit ? "yes" : "no");CHKERRQ(ierr); 1163 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 1164 } 1165 ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 1166 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1167 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "TSLoad_ARKIMEX" 1173 static PetscErrorCode TSLoad_ARKIMEX(TS ts,PetscViewer viewer) 1174 { 1175 PetscErrorCode ierr; 1176 SNES snes; 1177 TSAdapt tsadapt; 1178 1179 PetscFunctionBegin; 1180 ierr = TSGetTSAdapt(ts,&tsadapt);CHKERRQ(ierr); 1181 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 1182 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1183 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 1184 /* function and Jacobian context for SNES when used with TS is always ts object */ 1185 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 1186 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 1187 PetscFunctionReturn(0); 1188 } 1189 1190 #undef __FUNCT__ 1191 #define __FUNCT__ "TSARKIMEXSetType" 1192 /*@C 1193 TSARKIMEXSetType - Set the type of ARK IMEX scheme 1194 1195 Logically collective 1196 1197 Input Parameter: 1198 + ts - timestepping context 1199 - arktype - type of ARK-IMEX scheme 1200 1201 Level: intermediate 1202 1203 .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 1204 @*/ 1205 PetscErrorCode TSARKIMEXSetType(TS ts,TSARKIMEXType arktype) 1206 { 1207 PetscErrorCode ierr; 1208 1209 PetscFunctionBegin; 1210 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1211 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 1212 PetscFunctionReturn(0); 1213 } 1214 1215 #undef __FUNCT__ 1216 #define __FUNCT__ "TSARKIMEXGetType" 1217 /*@C 1218 TSARKIMEXGetType - Get the type of ARK IMEX scheme 1219 1220 Logically collective 1221 1222 Input Parameter: 1223 . ts - timestepping context 1224 1225 Output Parameter: 1226 . arktype - type of ARK-IMEX scheme 1227 1228 Level: intermediate 1229 1230 .seealso: TSARKIMEXGetType() 1231 @*/ 1232 PetscErrorCode TSARKIMEXGetType(TS ts,TSARKIMEXType *arktype) 1233 { 1234 PetscErrorCode ierr; 1235 1236 PetscFunctionBegin; 1237 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1238 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 1239 PetscFunctionReturn(0); 1240 } 1241 1242 #undef __FUNCT__ 1243 #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 1244 /*@C 1245 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 1246 1247 Logically collective 1248 1249 Input Parameter: 1250 + ts - timestepping context 1251 - flg - PETSC_TRUE for fully implicit 1252 1253 Level: intermediate 1254 1255 .seealso: TSARKIMEXGetType() 1256 @*/ 1257 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 1258 { 1259 PetscErrorCode ierr; 1260 1261 PetscFunctionBegin; 1262 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1263 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 1267 #undef __FUNCT__ 1268 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 1269 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType *arktype) 1270 { 1271 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1272 PetscErrorCode ierr; 1273 1274 PetscFunctionBegin; 1275 if (!ark->tableau) { 1276 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 1277 } 1278 *arktype = ark->tableau->name; 1279 PetscFunctionReturn(0); 1280 } 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 1283 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype) 1284 { 1285 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1286 PetscErrorCode ierr; 1287 PetscBool match; 1288 ARKTableauLink link; 1289 1290 PetscFunctionBegin; 1291 if (ark->tableau) { 1292 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 1293 if (match) PetscFunctionReturn(0); 1294 } 1295 for (link = ARKTableauList; link; link=link->next) { 1296 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 1297 if (match) { 1298 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 1299 ark->tableau = &link->tab; 1300 PetscFunctionReturn(0); 1301 } 1302 } 1303 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 1304 PetscFunctionReturn(0); 1305 } 1306 #undef __FUNCT__ 1307 #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 1308 PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 1309 { 1310 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1311 1312 PetscFunctionBegin; 1313 ark->imex = (PetscBool)!flg; 1314 PetscFunctionReturn(0); 1315 } 1316 1317 /* ------------------------------------------------------------ */ 1318 /*MC 1319 TSARKIMEX - ODE and DAE solver using Additive Runge-Kutta IMEX schemes 1320 1321 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1322 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1323 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1324 1325 Notes: 1326 The default is TSARKIMEX3, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 1327 1328 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). 1329 1330 Level: beginner 1331 1332 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 1333 TSARKIMEX4, TSARKIMEX5, TSARKIMEXPRSSP2, TSARKIMEXBPR3, TSARKIMEXType, TSARKIMEXRegister() 1334 1335 M*/ 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "TSCreate_ARKIMEX" 1338 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts) 1339 { 1340 TS_ARKIMEX *th; 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1345 ierr = TSARKIMEXInitializePackage();CHKERRQ(ierr); 1346 #endif 1347 1348 ts->ops->reset = TSReset_ARKIMEX; 1349 ts->ops->destroy = TSDestroy_ARKIMEX; 1350 ts->ops->view = TSView_ARKIMEX; 1351 ts->ops->load = TSLoad_ARKIMEX; 1352 ts->ops->setup = TSSetUp_ARKIMEX; 1353 ts->ops->step = TSStep_ARKIMEX; 1354 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1355 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1356 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1357 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1358 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1359 1360 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 1361 ts->data = (void*)th; 1362 th->imex = PETSC_TRUE; 1363 1364 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetType_C",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 1365 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetType_C",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 1366 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 1367 PetscFunctionReturn(0); 1368 } 1369