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