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