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 for (i=0; i<s; i++) { 644 if (At[i*s+i] == 0) { /* This stage is explicit */ 645 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 646 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 647 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 648 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 649 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 650 } else { 651 ark->stage_time = t + h*ct[i]; 652 ark->shift = 1./(h*At[i*s+i]); 653 /* Affine part */ 654 ierr = VecZeroEntries(W);CHKERRQ(ierr); 655 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 656 ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 657 ierr = VecScale(W, ark->shift);CHKERRQ(ierr); 658 659 /* Ydot = shift*(Y-Z) */ 660 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 661 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 662 ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr); 663 664 /* Initial guess taken from last stage */ 665 ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 666 ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 667 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 668 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 669 ts->snes_its += its; ts->ksp_its += lits; 670 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 671 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 672 if (!accept) goto reject_step; 673 } 674 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 675 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr); 676 if (ark->imex) { 677 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 678 } else { 679 ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr); 680 } 681 } 682 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 683 ark->status = TS_STEP_PENDING; 684 685 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 686 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 687 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 688 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 689 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 690 if (accept) { 691 /* ignore next_scheme for now */ 692 ts->ptime += ts->time_step; 693 ts->time_step = next_time_step; 694 ts->steps++; 695 ark->status = TS_STEP_COMPLETE; 696 break; 697 } else { /* Roll back the current step */ 698 for (j=0; j<s; j++) w[j] = h*bt[j]; 699 ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotI);CHKERRQ(ierr); 700 for (j=0; j<s; j++) w[j] = -h*b[j]; 701 ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotRHS);CHKERRQ(ierr); 702 ts->time_step = next_time_step; 703 ark->status = TS_STEP_INCOMPLETE; 704 } 705 reject_step: continue; 706 } 707 if (ark->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 708 PetscFunctionReturn(0); 709 } 710 711 #undef __FUNCT__ 712 #define __FUNCT__ "TSInterpolate_ARKIMEX" 713 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X) 714 { 715 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 716 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 717 PetscReal h; 718 PetscReal tt,t; 719 PetscScalar *bt,*b; 720 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 721 PetscErrorCode ierr; 722 723 PetscFunctionBegin; 724 if (!Bt || !B) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name); 725 switch (ark->status) { 726 case TS_STEP_INCOMPLETE: 727 case TS_STEP_PENDING: 728 h = ts->time_step; 729 t = (itime - ts->ptime)/h; 730 break; 731 case TS_STEP_COMPLETE: 732 h = ts->time_step_prev; 733 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 734 break; 735 default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 736 } 737 ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 738 for (i=0; i<s; i++) bt[i] = b[i] = 0; 739 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 740 for (i=0; i<s; i++) { 741 bt[i] += h * Bt[i*pinterp+j] * tt * -1.0; 742 b[i] += h * B[i*pinterp+j] * tt; 743 } 744 } 745 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"); 746 ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 747 ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 748 ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 749 ierr = PetscFree2(bt,b);CHKERRQ(ierr); 750 PetscFunctionReturn(0); 751 } 752 753 /*------------------------------------------------------------*/ 754 #undef __FUNCT__ 755 #define __FUNCT__ "TSReset_ARKIMEX" 756 static PetscErrorCode TSReset_ARKIMEX(TS ts) 757 { 758 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 759 PetscInt s; 760 PetscErrorCode ierr; 761 762 PetscFunctionBegin; 763 if (!ark->tableau) PetscFunctionReturn(0); 764 s = ark->tableau->s; 765 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 766 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 767 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 768 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 769 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 770 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 771 ierr = PetscFree(ark->work);CHKERRQ(ierr); 772 PetscFunctionReturn(0); 773 } 774 775 #undef __FUNCT__ 776 #define __FUNCT__ "TSDestroy_ARKIMEX" 777 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 778 { 779 PetscErrorCode ierr; 780 781 PetscFunctionBegin; 782 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 783 ierr = PetscFree(ts->data);CHKERRQ(ierr); 784 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr); 785 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr); 786 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","",PETSC_NULL);CHKERRQ(ierr); 787 PetscFunctionReturn(0); 788 } 789 790 /* 791 This defines the nonlinear equation that is to be solved with SNES 792 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 793 */ 794 #undef __FUNCT__ 795 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 796 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 797 { 798 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 799 PetscErrorCode ierr; 800 801 PetscFunctionBegin; 802 ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 803 ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,ark->imex);CHKERRQ(ierr); 804 PetscFunctionReturn(0); 805 } 806 807 #undef __FUNCT__ 808 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 809 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 810 { 811 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 812 PetscErrorCode ierr; 813 814 PetscFunctionBegin; 815 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 816 ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 817 PetscFunctionReturn(0); 818 } 819 820 #undef __FUNCT__ 821 #define __FUNCT__ "TSSetUp_ARKIMEX" 822 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 823 { 824 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 825 ARKTableau tab = ark->tableau; 826 PetscInt s = tab->s; 827 PetscErrorCode ierr; 828 829 PetscFunctionBegin; 830 if (!ark->tableau) { 831 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 832 } 833 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 834 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 835 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 836 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 837 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 838 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 839 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 /*------------------------------------------------------------*/ 843 844 #undef __FUNCT__ 845 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 846 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 847 { 848 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 849 PetscErrorCode ierr; 850 char arktype[256]; 851 852 PetscFunctionBegin; 853 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 854 { 855 ARKTableauLink link; 856 PetscInt count,choice; 857 PetscBool flg; 858 const char **namelist; 859 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr); 860 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 861 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 862 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 863 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 864 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 865 ierr = PetscFree(namelist);CHKERRQ(ierr); 866 flg = (PetscBool)!ark->imex; 867 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 868 ark->imex = (PetscBool)!flg; 869 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 870 } 871 ierr = PetscOptionsTail();CHKERRQ(ierr); 872 PetscFunctionReturn(0); 873 } 874 875 #undef __FUNCT__ 876 #define __FUNCT__ "PetscFormatRealArray" 877 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 878 { 879 PetscErrorCode ierr; 880 PetscInt i; 881 size_t left,count; 882 char *p; 883 884 PetscFunctionBegin; 885 for (i=0,p=buf,left=len; i<n; i++) { 886 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 887 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 888 left -= count; 889 p += count; 890 *p++ = ' '; 891 } 892 p[i ? 0 : -1] = 0; 893 PetscFunctionReturn(0); 894 } 895 896 #undef __FUNCT__ 897 #define __FUNCT__ "TSView_ARKIMEX" 898 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 899 { 900 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 901 ARKTableau tab = ark->tableau; 902 PetscBool iascii; 903 PetscErrorCode ierr; 904 905 PetscFunctionBegin; 906 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 907 if (iascii) { 908 const TSARKIMEXType arktype; 909 char buf[512]; 910 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 911 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 912 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 913 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 914 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 915 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 916 } 917 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 918 PetscFunctionReturn(0); 919 } 920 921 #undef __FUNCT__ 922 #define __FUNCT__ "TSARKIMEXSetType" 923 /*@C 924 TSARKIMEXSetType - Set the type of ARK IMEX scheme 925 926 Logically collective 927 928 Input Parameter: 929 + ts - timestepping context 930 - arktype - type of ARK-IMEX scheme 931 932 Level: intermediate 933 934 .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 935 @*/ 936 PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 937 { 938 PetscErrorCode ierr; 939 940 PetscFunctionBegin; 941 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 942 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 943 PetscFunctionReturn(0); 944 } 945 946 #undef __FUNCT__ 947 #define __FUNCT__ "TSARKIMEXGetType" 948 /*@C 949 TSARKIMEXGetType - Get the type of ARK IMEX scheme 950 951 Logically collective 952 953 Input Parameter: 954 . ts - timestepping context 955 956 Output Parameter: 957 . arktype - type of ARK-IMEX scheme 958 959 Level: intermediate 960 961 .seealso: TSARKIMEXGetType() 962 @*/ 963 PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 964 { 965 PetscErrorCode ierr; 966 967 PetscFunctionBegin; 968 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 969 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 970 PetscFunctionReturn(0); 971 } 972 973 #undef __FUNCT__ 974 #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 975 /*@C 976 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 977 978 Logically collective 979 980 Input Parameter: 981 + ts - timestepping context 982 - flg - PETSC_TRUE for fully implicit 983 984 Level: intermediate 985 986 .seealso: TSARKIMEXGetType() 987 @*/ 988 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 989 { 990 PetscErrorCode ierr; 991 992 PetscFunctionBegin; 993 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 994 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 995 PetscFunctionReturn(0); 996 } 997 998 EXTERN_C_BEGIN 999 #undef __FUNCT__ 1000 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 1001 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 1002 { 1003 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1004 PetscErrorCode ierr; 1005 1006 PetscFunctionBegin; 1007 if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 1008 *arktype = ark->tableau->name; 1009 PetscFunctionReturn(0); 1010 } 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 1013 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 1014 { 1015 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1016 PetscErrorCode ierr; 1017 PetscBool match; 1018 ARKTableauLink link; 1019 1020 PetscFunctionBegin; 1021 if (ark->tableau) { 1022 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 1023 if (match) PetscFunctionReturn(0); 1024 } 1025 for (link = ARKTableauList; link; link=link->next) { 1026 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 1027 if (match) { 1028 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 1029 ark->tableau = &link->tab; 1030 PetscFunctionReturn(0); 1031 } 1032 } 1033 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 1034 PetscFunctionReturn(0); 1035 } 1036 #undef __FUNCT__ 1037 #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 1038 PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 1039 { 1040 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1041 1042 PetscFunctionBegin; 1043 ark->imex = (PetscBool)!flg; 1044 PetscFunctionReturn(0); 1045 } 1046 EXTERN_C_END 1047 1048 /* ------------------------------------------------------------ */ 1049 /*MC 1050 TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 1051 1052 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1053 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1054 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1055 1056 Notes: 1057 The default is TSARKIMEX2E, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 1058 1059 This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 1060 1061 Level: beginner 1062 1063 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 1064 TSARKIMEX4, TSARKIMEX5, TSARKIMEXType, TSARKIMEXRegister() 1065 1066 M*/ 1067 EXTERN_C_BEGIN 1068 #undef __FUNCT__ 1069 #define __FUNCT__ "TSCreate_ARKIMEX" 1070 PetscErrorCode TSCreate_ARKIMEX(TS ts) 1071 { 1072 TS_ARKIMEX *th; 1073 PetscErrorCode ierr; 1074 1075 PetscFunctionBegin; 1076 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1077 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1078 #endif 1079 1080 ts->ops->reset = TSReset_ARKIMEX; 1081 ts->ops->destroy = TSDestroy_ARKIMEX; 1082 ts->ops->view = TSView_ARKIMEX; 1083 ts->ops->setup = TSSetUp_ARKIMEX; 1084 ts->ops->step = TSStep_ARKIMEX; 1085 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1086 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1087 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1088 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1089 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1090 1091 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 1092 ts->data = (void*)th; 1093 th->imex = PETSC_TRUE; 1094 1095 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 1096 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 1097 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 1098 PetscFunctionReturn(0); 1099 } 1100 EXTERN_C_END 1101