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 14 static 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 scoeff; /* shift = scoeff/dt */ 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(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->scoeff = 1./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->scoeff/h);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 #undef __FUNCT__ 794 #define __FUNCT__ "TSARKIMEXGetVecs" 795 static PetscErrorCode TSARKIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot) 796 { 797 TS_ARKIMEX *ax = (TS_ARKIMEX*)ts->data; 798 PetscErrorCode ierr; 799 800 PetscFunctionBegin; 801 if (Z) { 802 if (dm && dm != ts->dm) { 803 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr); 804 } else *Z = ax->Z; 805 } 806 if (Ydot) { 807 if (dm && dm != ts->dm) { 808 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr); 809 } else *Ydot = ax->Ydot; 810 } 811 PetscFunctionReturn(0); 812 } 813 814 815 #undef __FUNCT__ 816 #define __FUNCT__ "TSARKIMEXRestoreVecs" 817 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot) 818 { 819 PetscErrorCode ierr; 820 821 PetscFunctionBegin; 822 if (Z) { 823 if (dm && dm != ts->dm) { 824 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr); 825 } 826 } 827 if (Ydot) { 828 if (dm && dm != ts->dm) { 829 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr); 830 } 831 } 832 PetscFunctionReturn(0); 833 } 834 835 /* 836 This defines the nonlinear equation that is to be solved with SNES 837 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 838 */ 839 #undef __FUNCT__ 840 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 841 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 842 { 843 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 844 DM dm,dmsave; 845 Vec Z,Ydot; 846 PetscReal shift = ark->scoeff / ts->time_step; 847 PetscErrorCode ierr; 848 849 PetscFunctionBegin; 850 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 851 ierr = TSARKIMEXGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 852 ierr = VecAXPBYPCZ(Ydot,-shift,shift,0,Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 853 dmsave = ts->dm; 854 ts->dm = dm; 855 ierr = TSComputeIFunction(ts,ark->stage_time,X,Ydot,F,ark->imex);CHKERRQ(ierr); 856 ts->dm = dmsave; 857 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 858 PetscFunctionReturn(0); 859 } 860 861 #undef __FUNCT__ 862 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 863 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 864 { 865 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 866 DM dm,dmsave; 867 Vec Ydot; 868 PetscReal shift = ark->scoeff / ts->time_step; 869 PetscErrorCode ierr; 870 871 PetscFunctionBegin; 872 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 873 ierr = TSARKIMEXGetVecs(ts,dm,PETSC_NULL,&Ydot);CHKERRQ(ierr); 874 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 875 dmsave = ts->dm; 876 ts->dm = dm; 877 ierr = TSComputeIJacobian(ts,ark->stage_time,X,Ydot,shift,A,B,str,ark->imex);CHKERRQ(ierr); 878 ts->dm = dmsave; 879 ierr = TSARKIMEXRestoreVecs(ts,dm,PETSC_NULL,&Ydot);CHKERRQ(ierr); 880 PetscFunctionReturn(0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "DMCoarsenHook_TSARKIMEX" 885 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine,DM coarse,void *ctx) 886 { 887 888 PetscFunctionBegin; 889 PetscFunctionReturn(0); 890 } 891 892 #undef __FUNCT__ 893 #define __FUNCT__ "DMRestrictHook_TSARKIMEX" 894 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 895 { 896 TS ts = (TS)ctx; 897 PetscErrorCode ierr; 898 Vec Z,Z_c; 899 900 PetscFunctionBegin; 901 ierr = TSARKIMEXGetVecs(ts,fine,&Z,PETSC_NULL);CHKERRQ(ierr); 902 ierr = TSARKIMEXGetVecs(ts,coarse,&Z_c,PETSC_NULL);CHKERRQ(ierr); 903 ierr = MatRestrict(restrct,Z,Z_c);CHKERRQ(ierr); 904 ierr = VecPointwiseMult(Z_c,rscale,Z_c);CHKERRQ(ierr); 905 ierr = TSARKIMEXRestoreVecs(ts,fine,&Z,PETSC_NULL);CHKERRQ(ierr); 906 ierr = TSARKIMEXRestoreVecs(ts,coarse,&Z_c,PETSC_NULL);CHKERRQ(ierr); 907 PetscFunctionReturn(0); 908 } 909 910 #undef __FUNCT__ 911 #define __FUNCT__ "TSSetUp_ARKIMEX" 912 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 913 { 914 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 915 ARKTableau tab = ark->tableau; 916 PetscInt s = tab->s; 917 PetscErrorCode ierr; 918 DM dm; 919 920 PetscFunctionBegin; 921 if (!ark->tableau) { 922 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 923 } 924 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 925 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 926 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 927 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 928 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 929 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 930 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 931 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 932 if (dm) { 933 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSARKIMEX,DMRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr); 934 } 935 PetscFunctionReturn(0); 936 } 937 /*------------------------------------------------------------*/ 938 939 #undef __FUNCT__ 940 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 941 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 942 { 943 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 944 PetscErrorCode ierr; 945 char arktype[256]; 946 947 PetscFunctionBegin; 948 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 949 { 950 ARKTableauLink link; 951 PetscInt count,choice; 952 PetscBool flg; 953 const char **namelist; 954 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof(arktype));CHKERRQ(ierr); 955 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 956 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 957 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 958 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 959 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 960 ierr = PetscFree(namelist);CHKERRQ(ierr); 961 flg = (PetscBool)!ark->imex; 962 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 963 ark->imex = (PetscBool)!flg; 964 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 965 } 966 ierr = PetscOptionsTail();CHKERRQ(ierr); 967 PetscFunctionReturn(0); 968 } 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "PetscFormatRealArray" 972 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 973 { 974 PetscErrorCode ierr; 975 PetscInt i; 976 size_t left,count; 977 char *p; 978 979 PetscFunctionBegin; 980 for (i=0,p=buf,left=len; i<n; i++) { 981 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 982 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 983 left -= count; 984 p += count; 985 *p++ = ' '; 986 } 987 p[i ? 0 : -1] = 0; 988 PetscFunctionReturn(0); 989 } 990 991 #undef __FUNCT__ 992 #define __FUNCT__ "TSView_ARKIMEX" 993 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 994 { 995 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 996 ARKTableau tab = ark->tableau; 997 PetscBool iascii; 998 PetscErrorCode ierr; 999 TSAdapt adapt; 1000 1001 PetscFunctionBegin; 1002 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1003 if (iascii) { 1004 TSARKIMEXType arktype; 1005 char buf[512]; 1006 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 1007 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 1008 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 1009 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 1010 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1011 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 1012 } 1013 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1014 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1015 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "TSARKIMEXSetType" 1021 /*@C 1022 TSARKIMEXSetType - Set the type of ARK IMEX scheme 1023 1024 Logically collective 1025 1026 Input Parameter: 1027 + ts - timestepping context 1028 - arktype - type of ARK-IMEX scheme 1029 1030 Level: intermediate 1031 1032 .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 1033 @*/ 1034 PetscErrorCode TSARKIMEXSetType(TS ts,TSARKIMEXType arktype) 1035 { 1036 PetscErrorCode ierr; 1037 1038 PetscFunctionBegin; 1039 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1040 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 1041 PetscFunctionReturn(0); 1042 } 1043 1044 #undef __FUNCT__ 1045 #define __FUNCT__ "TSARKIMEXGetType" 1046 /*@C 1047 TSARKIMEXGetType - Get the type of ARK IMEX scheme 1048 1049 Logically collective 1050 1051 Input Parameter: 1052 . ts - timestepping context 1053 1054 Output Parameter: 1055 . arktype - type of ARK-IMEX scheme 1056 1057 Level: intermediate 1058 1059 .seealso: TSARKIMEXGetType() 1060 @*/ 1061 PetscErrorCode TSARKIMEXGetType(TS ts,TSARKIMEXType *arktype) 1062 { 1063 PetscErrorCode ierr; 1064 1065 PetscFunctionBegin; 1066 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1067 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 1073 /*@C 1074 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 1075 1076 Logically collective 1077 1078 Input Parameter: 1079 + ts - timestepping context 1080 - flg - PETSC_TRUE for fully implicit 1081 1082 Level: intermediate 1083 1084 .seealso: TSARKIMEXGetType() 1085 @*/ 1086 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 1087 { 1088 PetscErrorCode ierr; 1089 1090 PetscFunctionBegin; 1091 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1092 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 EXTERN_C_BEGIN 1097 #undef __FUNCT__ 1098 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 1099 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType *arktype) 1100 { 1101 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1102 PetscErrorCode ierr; 1103 1104 PetscFunctionBegin; 1105 if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 1106 *arktype = ark->tableau->name; 1107 PetscFunctionReturn(0); 1108 } 1109 #undef __FUNCT__ 1110 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 1111 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype) 1112 { 1113 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1114 PetscErrorCode ierr; 1115 PetscBool match; 1116 ARKTableauLink link; 1117 1118 PetscFunctionBegin; 1119 if (ark->tableau) { 1120 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 1121 if (match) PetscFunctionReturn(0); 1122 } 1123 for (link = ARKTableauList; link; link=link->next) { 1124 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 1125 if (match) { 1126 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 1127 ark->tableau = &link->tab; 1128 PetscFunctionReturn(0); 1129 } 1130 } 1131 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 1132 PetscFunctionReturn(0); 1133 } 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 1136 PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 1137 { 1138 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1139 1140 PetscFunctionBegin; 1141 ark->imex = (PetscBool)!flg; 1142 PetscFunctionReturn(0); 1143 } 1144 EXTERN_C_END 1145 1146 /* ------------------------------------------------------------ */ 1147 /*MC 1148 TSARKIMEX - ODE and DAE solver using Additive Runge-Kutta IMEX schemes 1149 1150 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1151 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1152 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1153 1154 Notes: 1155 The default is TSARKIMEX3, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 1156 1157 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). 1158 1159 Level: beginner 1160 1161 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 1162 TSARKIMEX4, TSARKIMEX5, TSARKIMEXPRSSP2, TSARKIMEXBPR3, TSARKIMEXType, TSARKIMEXRegister() 1163 1164 M*/ 1165 EXTERN_C_BEGIN 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "TSCreate_ARKIMEX" 1168 PetscErrorCode TSCreate_ARKIMEX(TS ts) 1169 { 1170 TS_ARKIMEX *th; 1171 PetscErrorCode ierr; 1172 1173 PetscFunctionBegin; 1174 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1175 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1176 #endif 1177 1178 ts->ops->reset = TSReset_ARKIMEX; 1179 ts->ops->destroy = TSDestroy_ARKIMEX; 1180 ts->ops->view = TSView_ARKIMEX; 1181 ts->ops->setup = TSSetUp_ARKIMEX; 1182 ts->ops->step = TSStep_ARKIMEX; 1183 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1184 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1185 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1186 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1187 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1188 1189 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 1190 ts->data = (void*)th; 1191 th->imex = PETSC_TRUE; 1192 1193 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 1194 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 1195 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 1196 PetscFunctionReturn(0); 1197 } 1198 EXTERN_C_END 1199