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