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