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