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