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