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