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