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 = TSGetTSAdapt(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 = TSGetTSAdapt(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; 916 PetscInt s; 917 PetscErrorCode ierr; 918 DM dm; 919 920 PetscFunctionBegin; 921 if (!ark->tableau) { 922 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 923 } 924 tab = ark->tableau; 925 s = tab->s; 926 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 927 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 928 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 929 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 930 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 931 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 932 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 933 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 934 if (dm) { 935 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSARKIMEX,DMRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr); 936 } 937 PetscFunctionReturn(0); 938 } 939 /*------------------------------------------------------------*/ 940 941 #undef __FUNCT__ 942 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 943 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 944 { 945 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 946 PetscErrorCode ierr; 947 char arktype[256]; 948 949 PetscFunctionBegin; 950 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 951 { 952 ARKTableauLink link; 953 PetscInt count,choice; 954 PetscBool flg; 955 const char **namelist; 956 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof(arktype));CHKERRQ(ierr); 957 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 958 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 959 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 960 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 961 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 962 ierr = PetscFree(namelist);CHKERRQ(ierr); 963 flg = (PetscBool)!ark->imex; 964 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 965 ark->imex = (PetscBool)!flg; 966 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 967 } 968 ierr = PetscOptionsTail();CHKERRQ(ierr); 969 PetscFunctionReturn(0); 970 } 971 972 #undef __FUNCT__ 973 #define __FUNCT__ "PetscFormatRealArray" 974 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 975 { 976 PetscErrorCode ierr; 977 PetscInt i; 978 size_t left,count; 979 char *p; 980 981 PetscFunctionBegin; 982 for (i=0,p=buf,left=len; i<n; i++) { 983 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 984 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 985 left -= count; 986 p += count; 987 *p++ = ' '; 988 } 989 p[i ? 0 : -1] = 0; 990 PetscFunctionReturn(0); 991 } 992 993 #undef __FUNCT__ 994 #define __FUNCT__ "TSView_ARKIMEX" 995 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 996 { 997 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 998 ARKTableau tab = ark->tableau; 999 PetscBool iascii; 1000 PetscErrorCode ierr; 1001 TSAdapt adapt; 1002 1003 PetscFunctionBegin; 1004 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1005 if (iascii) { 1006 TSARKIMEXType arktype; 1007 char buf[512]; 1008 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 1009 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 1010 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 1011 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 1012 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1013 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 1014 } 1015 ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 1016 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1017 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1018 PetscFunctionReturn(0); 1019 } 1020 1021 #undef __FUNCT__ 1022 #define __FUNCT__ "TSLoad_ARKIMEX" 1023 static PetscErrorCode TSLoad_ARKIMEX(TS ts,PetscViewer viewer) 1024 { 1025 PetscErrorCode ierr; 1026 SNES snes; 1027 TSAdapt tsadapt; 1028 1029 PetscFunctionBegin; 1030 ierr = TSGetTSAdapt(ts,&tsadapt);CHKERRQ(ierr); 1031 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 1032 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1033 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 1034 /* function and Jacobian context for SNES when used with TS is always ts object */ 1035 ierr = SNESSetFunction(snes,PETSC_NULL,PETSC_NULL,ts);CHKERRQ(ierr); 1036 ierr = SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,PETSC_NULL,ts);CHKERRQ(ierr); 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "TSARKIMEXSetType" 1042 /*@C 1043 TSARKIMEXSetType - Set the type of ARK IMEX scheme 1044 1045 Logically collective 1046 1047 Input Parameter: 1048 + ts - timestepping context 1049 - arktype - type of ARK-IMEX scheme 1050 1051 Level: intermediate 1052 1053 .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 1054 @*/ 1055 PetscErrorCode TSARKIMEXSetType(TS ts,TSARKIMEXType arktype) 1056 { 1057 PetscErrorCode ierr; 1058 1059 PetscFunctionBegin; 1060 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1061 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 #undef __FUNCT__ 1066 #define __FUNCT__ "TSARKIMEXGetType" 1067 /*@C 1068 TSARKIMEXGetType - Get the type of ARK IMEX scheme 1069 1070 Logically collective 1071 1072 Input Parameter: 1073 . ts - timestepping context 1074 1075 Output Parameter: 1076 . arktype - type of ARK-IMEX scheme 1077 1078 Level: intermediate 1079 1080 .seealso: TSARKIMEXGetType() 1081 @*/ 1082 PetscErrorCode TSARKIMEXGetType(TS ts,TSARKIMEXType *arktype) 1083 { 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1088 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 1094 /*@C 1095 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 1096 1097 Logically collective 1098 1099 Input Parameter: 1100 + ts - timestepping context 1101 - flg - PETSC_TRUE for fully implicit 1102 1103 Level: intermediate 1104 1105 .seealso: TSARKIMEXGetType() 1106 @*/ 1107 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 1108 { 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1113 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1114 PetscFunctionReturn(0); 1115 } 1116 1117 EXTERN_C_BEGIN 1118 #undef __FUNCT__ 1119 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 1120 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType *arktype) 1121 { 1122 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1123 PetscErrorCode ierr; 1124 1125 PetscFunctionBegin; 1126 if (!ark->tableau) { 1127 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 1128 } 1129 *arktype = ark->tableau->name; 1130 PetscFunctionReturn(0); 1131 } 1132 #undef __FUNCT__ 1133 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 1134 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype) 1135 { 1136 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1137 PetscErrorCode ierr; 1138 PetscBool match; 1139 ARKTableauLink link; 1140 1141 PetscFunctionBegin; 1142 if (ark->tableau) { 1143 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 1144 if (match) PetscFunctionReturn(0); 1145 } 1146 for (link = ARKTableauList; link; link=link->next) { 1147 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 1148 if (match) { 1149 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 1150 ark->tableau = &link->tab; 1151 PetscFunctionReturn(0); 1152 } 1153 } 1154 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 1155 PetscFunctionReturn(0); 1156 } 1157 #undef __FUNCT__ 1158 #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 1159 PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 1160 { 1161 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 1162 1163 PetscFunctionBegin; 1164 ark->imex = (PetscBool)!flg; 1165 PetscFunctionReturn(0); 1166 } 1167 EXTERN_C_END 1168 1169 /* ------------------------------------------------------------ */ 1170 /*MC 1171 TSARKIMEX - ODE and DAE solver using Additive Runge-Kutta IMEX schemes 1172 1173 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1174 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1175 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1176 1177 Notes: 1178 The default is TSARKIMEX3, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 1179 1180 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). 1181 1182 Level: beginner 1183 1184 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 1185 TSARKIMEX4, TSARKIMEX5, TSARKIMEXPRSSP2, TSARKIMEXBPR3, TSARKIMEXType, TSARKIMEXRegister() 1186 1187 M*/ 1188 EXTERN_C_BEGIN 1189 #undef __FUNCT__ 1190 #define __FUNCT__ "TSCreate_ARKIMEX" 1191 PetscErrorCode TSCreate_ARKIMEX(TS ts) 1192 { 1193 TS_ARKIMEX *th; 1194 PetscErrorCode ierr; 1195 1196 PetscFunctionBegin; 1197 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1198 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1199 #endif 1200 1201 ts->ops->reset = TSReset_ARKIMEX; 1202 ts->ops->destroy = TSDestroy_ARKIMEX; 1203 ts->ops->view = TSView_ARKIMEX; 1204 ts->ops->load = TSLoad_ARKIMEX; 1205 ts->ops->setup = TSSetUp_ARKIMEX; 1206 ts->ops->step = TSStep_ARKIMEX; 1207 ts->ops->interpolate = TSInterpolate_ARKIMEX; 1208 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 1209 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 1210 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 1211 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 1212 1213 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 1214 ts->data = (void*)th; 1215 th->imex = PETSC_TRUE; 1216 1217 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 1218 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 1219 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 1220 PetscFunctionReturn(0); 1221 } 1222 EXTERN_C_END 1223