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