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