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