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 = TSARKIMEX2E; 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 *binterpt,*binterp; /* Dense output formula */ 27 }; 28 typedef struct _ARKTableauLink *ARKTableauLink; 29 struct _ARKTableauLink { 30 struct _ARKTableau tab; 31 ARKTableauLink next; 32 }; 33 static ARKTableauLink ARKTableauList; 34 35 typedef struct { 36 ARKTableau tableau; 37 Vec *Y; /* States computed during the step */ 38 Vec *YdotI; /* Time derivatives for the stiff part */ 39 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 40 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 41 Vec Work; /* Generic work vector */ 42 Vec Z; /* Ydot = shift(Y-Z) */ 43 PetscScalar *work; /* Scalar work */ 44 PetscReal shift; 45 PetscReal stage_time; 46 PetscBool imex; 47 } TS_ARKIMEX; 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "TSARKIMEXRegisterAll" 51 /*@C 52 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSARKIMEX 53 54 Not Collective, but should be called by all processes which will need the schemes to be registered 55 56 Level: advanced 57 58 .keywords: TS, TSARKIMEX, register, all 59 60 .seealso: TSARKIMEXRegisterDestroy() 61 @*/ 62 PetscErrorCode TSARKIMEXRegisterAll(void) 63 { 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(0); 68 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 69 70 { 71 const PetscReal 72 A[3][3] = {{0,0,0}, 73 {0.41421356237309504880,0,0}, 74 {0.75,0.25,0}}, 75 At[3][3] = {{0,0,0}, 76 {0.12132034355964257320,0.29289321881345247560,0}, 77 {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}}, 78 binterpt[3][2] = {{1,-0.5},{0,0},{0,0.5}}; 79 ierr = TSARKIMEXRegister(TSARKIMEX2D,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 80 } 81 { /* Optimal for linear implicit part */ 82 const PetscReal s2 = sqrt(2), 83 A[3][3] = {{0,0,0}, 84 {2-s2,0,0}, 85 {(3-2*s2)/6,(3+2*s2)/6,0}}, 86 At[3][3] = {{0,0,0}, 87 {1-1/s2,1-1/s2,0}, 88 {1/(2*s2),1/(2*s2),1-1/s2}}, 89 binterpt[3][2] = {{1,-0.5},{0,0},{0,0.5}}; 90 ierr = TSARKIMEXRegister(TSARKIMEX2E,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 91 } 92 { 93 const PetscReal 94 A[4][4] = {{0,0,0,0}, 95 {1767732205903./2027836641118.,0,0,0}, 96 {5535828885825./10492691773637.,788022342437./10882634858940.,0,0}, 97 {6485989280629./16251701735622.,-4246266847089./9704473918619.,10755448449292./10357097424841.,0}}, 98 At[4][4] = {{0,0,0,0}, 99 {1767732205903./4055673282236.,1767732205903./4055673282236.,0,0}, 100 {2746238789719./10658868560708.,-640167445237./6845629431997.,1767732205903./4055673282236.,0}, 101 {1471266399579./7840856788654.,-4482444167858./7529755066697.,11266239266428./11593286722821.,1767732205903./4055673282236.}}, 102 binterpt[4][2] = {{4655552711362./22874653954995., -215264564351./13552729205753.}, 103 {-18682724506714./9892148508045.,17870216137069./13817060693119.}, 104 {34259539580243./13192909600954.,-28141676662227./17317692491321.}, 105 {584795268549./6622622206610., 2508943948391./7218656332882.}}; 106 ierr = TSARKIMEXRegister(TSARKIMEX3,3,4,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 107 } 108 { 109 const PetscReal 110 A[6][6] = {{0,0,0,0,0,0}, 111 {1./2,0,0,0,0,0}, 112 {13861./62500.,6889./62500.,0,0,0,0}, 113 {-116923316275./2393684061468.,-2731218467317./15368042101831.,9408046702089./11113171139209.,0,0,0}, 114 {-451086348788./2902428689909.,-2682348792572./7519795681897.,12662868775082./11960479115383.,3355817975965./11060851509271.,0,0}, 115 {647845179188./3216320057751.,73281519250./8382639484533.,552539513391./3454668386233.,3354512671639./8306763924573.,4040./17871.,0}}, 116 At[6][6] = {{0,0,0,0,0,0}, 117 {1./4,1./4,0,0,0,0}, 118 {8611./62500.,-1743./31250.,1./4,0,0,0}, 119 {5012029./34652500.,-654441./2922500.,174375./388108.,1./4,0,0}, 120 {15267082809./155376265600.,-71443401./120774400.,730878875./902184768.,2285395./8070912.,1./4,0}, 121 {82889./524892.,0,15625./83664.,69875./102672.,-2260./8211,1./4}}, 122 binterpt[6][3] = {{6943876665148./7220017795957.,-54480133./30881146.,6818779379841./7100303317025.}, 123 {0,0,0}, 124 {7640104374378./9702883013639.,-11436875./14766696.,2173542590792./12501825683035.}, 125 {-20649996744609./7521556579894.,174696575./18121608.,-31592104683404./5083833661969.}, 126 {8854892464581./2390941311638.,-12120380./966161.,61146701046299./7138195549469.}, 127 {-11397109935349./6675773540249.,3843./706.,-17219254887155./4939391667607.}}; 128 ierr = TSARKIMEXRegister(TSARKIMEX4,4,6,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,3,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 129 } 130 { 131 const PetscReal 132 A[8][8] = {{0,0,0,0,0,0,0,0}, 133 {41./100,0,0,0,0,0,0,0}, 134 {367902744464./2072280473677.,677623207551./8224143866563.,0,0,0,0,0,0}, 135 {1268023523408./10340822734521.,0,1029933939417./13636558850479.,0,0,0,0,0}, 136 {14463281900351./6315353703477.,0,66114435211212./5879490589093.,-54053170152839./4284798021562.,0,0,0,0}, 137 {14090043504691./34967701212078.,0,15191511035443./11219624916014.,-18461159152457./12425892160975.,-281667163811./9011619295870.,0,0,0}, 138 {19230459214898./13134317526959.,0,21275331358303./2942455364971.,-38145345988419./4862620318723.,-1./8,-1./8,0,0}, 139 {-19977161125411./11928030595625.,0,-40795976796054./6384907823539.,177454434618887./12078138498510.,782672205425./8267701900261.,-69563011059811./9646580694205.,7356628210526./4942186776405.,0}}, 140 At[8][8] = {{0,0,0,0,0,0,0,0}, 141 {41./200.,41./200.,0,0,0,0,0,0}, 142 {41./400.,-567603406766./11931857230679.,41./200.,0,0,0,0,0}, 143 {683785636431./9252920307686.,0,-110385047103./1367015193373.,41./200.,0,0,0,0}, 144 {3016520224154./10081342136671.,0,30586259806659./12414158314087.,-22760509404356./11113319521817.,41./200.,0,0,0}, 145 {218866479029./1489978393911.,0,638256894668./5436446318841.,-1179710474555./5321154724896.,-60928119172./8023461067671.,41./200.,0,0}, 146 {1020004230633./5715676835656.,0,25762820946817./25263940353407.,-2161375909145./9755907335909.,-211217309593./5846859502534.,-4269925059573./7827059040749.,41./200,0}, 147 {-872700587467./9133579230613.,0,0,22348218063261./9555858737531.,-1143369518992./8141816002931.,-39379526789629./19018526304540.,32727382324388./42900044865799.,41./200.}}, 148 binterpt[8][3] = {{-17674230611817./10670229744614. , 43486358583215./12773830924787. , -9257016797708./5021505065439.}, 149 {0 , 0 , 0 }, 150 {0 , 0 , 0 }, 151 {65168852399939./7868540260826. , -91478233927265./11067650958493., 26096422576131./11239449250142.}, 152 {15494834004392./5936557850923. , -79368583304911./10890268929626., 92396832856987./20362823103730.}, 153 {-99329723586156./26959484932159., -12239297817655./9152339842473. , 30029262896817./10175596800299.}, 154 {-19024464361622./5461577185407. , 115839755401235./10719374521269., -26136350496073./3983972220547.}, 155 {-6511271360970./6095937251113. , 5843115559534./2180450260947. , -5289405421727./3760307252460. }}; 156 ierr = TSARKIMEXRegister(TSARKIMEX5,5,8,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,3,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 157 } 158 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "TSARKIMEXRegisterDestroy" 164 /*@C 165 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister(). 166 167 Not Collective 168 169 Level: advanced 170 171 .keywords: TSARKIMEX, register, destroy 172 .seealso: TSARKIMEXRegister(), TSARKIMEXRegisterAll(), TSARKIMEXRegisterDynamic() 173 @*/ 174 PetscErrorCode TSARKIMEXRegisterDestroy(void) 175 { 176 PetscErrorCode ierr; 177 ARKTableauLink link; 178 179 PetscFunctionBegin; 180 while ((link = ARKTableauList)) { 181 ARKTableau t = &link->tab; 182 ARKTableauList = link->next; 183 ierr = PetscFree6(t->At,t->bt,t->ct,t->A,t->b,t->c);CHKERRQ(ierr); 184 ierr = PetscFree2(t->binterpt,t->binterp);CHKERRQ(ierr); 185 ierr = PetscFree(t->name);CHKERRQ(ierr); 186 ierr = PetscFree(link);CHKERRQ(ierr); 187 } 188 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 189 PetscFunctionReturn(0); 190 } 191 192 #undef __FUNCT__ 193 #define __FUNCT__ "TSARKIMEXInitializePackage" 194 /*@C 195 TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called 196 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_ARKIMEX() 197 when using static libraries. 198 199 Input Parameter: 200 path - The dynamic library path, or PETSC_NULL 201 202 Level: developer 203 204 .keywords: TS, TSARKIMEX, initialize, package 205 .seealso: PetscInitialize() 206 @*/ 207 PetscErrorCode TSARKIMEXInitializePackage(const char path[]) 208 { 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0); 213 TSARKIMEXPackageInitialized = PETSC_TRUE; 214 ierr = TSARKIMEXRegisterAll();CHKERRQ(ierr); 215 ierr = PetscRegisterFinalize(TSARKIMEXFinalizePackage);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "TSARKIMEXFinalizePackage" 221 /*@C 222 TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is 223 called from PetscFinalize(). 224 225 Level: developer 226 227 .keywords: Petsc, destroy, package 228 .seealso: PetscFinalize() 229 @*/ 230 PetscErrorCode TSARKIMEXFinalizePackage(void) 231 { 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 TSARKIMEXPackageInitialized = PETSC_FALSE; 236 ierr = TSARKIMEXRegisterDestroy();CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "TSARKIMEXRegister" 242 /*@C 243 TSARKIMEXRegister - register an ARK IMEX scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 244 245 Not Collective, but the same schemes should be registered on all processes on which they will be used 246 247 Input Parameters: 248 + name - identifier for method 249 . order - approximation order of method 250 . s - number of stages, this is the dimension of the matrices below 251 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 252 . bt - Butcher table for completing the stiff part of the step (dimension s; PETSC_NULL to use the last row of At) 253 . ct - Abscissa of each stiff stage (dimension s, PETSC_NULL to use row sums of At) 254 . A - Non-stiff stage coefficients (dimension s*s, row-major) 255 . b - Non-stiff step completion table (dimension s; PETSC_NULL to use last row of At) 256 . c - Non-stiff abscissa (dimension s; PETSC_NULL to use row sums of A) 257 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 258 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 259 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; PETSC_NULL to reuse binterpt) 260 261 Notes: 262 Several ARK IMEX methods are provided, this function is only needed to create new methods. 263 264 Level: advanced 265 266 .keywords: TS, register 267 268 .seealso: TSARKIMEX 269 @*/ 270 PetscErrorCode TSARKIMEXRegister(const TSARKIMEXType name,PetscInt order,PetscInt s, 271 const PetscReal At[],const PetscReal bt[],const PetscReal ct[], 272 const PetscReal A[],const PetscReal b[],const PetscReal c[], 273 PetscInt pinterp,const PetscReal binterpt[],const PetscReal binterp[]) 274 { 275 PetscErrorCode ierr; 276 ARKTableauLink link; 277 ARKTableau t; 278 PetscInt i,j; 279 280 PetscFunctionBegin; 281 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 282 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 283 t = &link->tab; 284 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 285 t->order = order; 286 t->s = s; 287 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); 288 ierr = PetscMemcpy(t->At,At,s*s*sizeof(At[0]));CHKERRQ(ierr); 289 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 290 if (bt) {ierr = PetscMemcpy(t->bt,bt,s*sizeof(bt[0]));CHKERRQ(ierr);} 291 else for (i=0; i<s; i++) t->bt[i] = At[(s-1)*s+i]; 292 if (b) {ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr);} 293 else for (i=0; i<s; i++) t->b[i] = At[(s-1)*s+i]; 294 if (ct) {ierr = PetscMemcpy(t->ct,ct,s*sizeof(ct[0]));CHKERRQ(ierr);} 295 else for (i=0; i<s; i++) for (j=0,t->ct[i]=0; j<s; j++) t->ct[i] += At[i*s+j]; 296 if (c) {ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr);} 297 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 298 t->pinterp = pinterp; 299 ierr = PetscMalloc2(s*pinterp,PetscReal,&t->binterpt,s*pinterp,PetscReal,&t->binterp);CHKERRQ(ierr); 300 ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 301 ierr = PetscMemcpy(t->binterp,binterp?binterp:binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 302 link->next = ARKTableauList; 303 ARKTableauList = link; 304 PetscFunctionReturn(0); 305 } 306 307 #undef __FUNCT__ 308 #define __FUNCT__ "TSStep_ARKIMEX" 309 static PetscErrorCode TSStep_ARKIMEX(TS ts) 310 { 311 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 312 ARKTableau tab = ark->tableau; 313 const PetscInt s = tab->s; 314 const PetscReal *At = tab->At,*A = tab->A,*bt = tab->bt,*b = tab->b,*ct = tab->ct,*c = tab->c; 315 PetscScalar *w = ark->work; 316 Vec *Y = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,W = ark->Work,Z = ark->Z; 317 SNES snes; 318 PetscInt i,j,its,lits; 319 PetscReal h,t; 320 PetscErrorCode ierr; 321 322 PetscFunctionBegin; 323 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 324 h = ts->time_step = ts->next_time_step; 325 t = ts->ptime; 326 327 for (i=0; i<s; i++) { 328 if (At[i*s+i] == 0) { /* This stage is explicit */ 329 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 330 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 331 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 332 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 333 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 334 } else { 335 ark->stage_time = t + h*ct[i]; 336 ark->shift = 1./(h*At[i*s+i]); 337 /* Affine part */ 338 ierr = VecZeroEntries(W);CHKERRQ(ierr); 339 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 340 ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 341 /* Ydot = shift*(Y-Z) */ 342 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 343 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 344 ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr); 345 /* Initial guess taken from last stage */ 346 ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 347 ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 348 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 349 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 350 ts->nonlinear_its += its; ts->linear_its += lits; 351 } 352 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 353 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr); 354 if (ark->imex) { 355 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 356 } else { 357 ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr); 358 } 359 } 360 for (j=0; j<s; j++) w[j] = -h*bt[j]; 361 ierr = VecMAXPY(ts->vec_sol,s,w,YdotI);CHKERRQ(ierr); 362 for (j=0; j<s; j++) w[j] = h*b[j]; 363 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 364 365 ts->ptime += ts->time_step; 366 ts->next_time_step = ts->time_step; 367 ts->steps++; 368 PetscFunctionReturn(0); 369 } 370 371 #undef __FUNCT__ 372 #define __FUNCT__ "TSInterpolate_ARKIMEX" 373 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X) 374 { 375 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 376 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 377 PetscReal tt,t = (itime - ts->ptime)/ts->time_step + 1; /* In the interval [0,1] */ 378 PetscScalar *bt,*b; 379 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 380 PetscErrorCode ierr; 381 382 PetscFunctionBegin; 383 if (!Bt || !B) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name); 384 ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 385 for (i=0; i<s; i++) bt[i] = b[i] = 0; 386 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 387 for (i=0; i<s; i++) { 388 bt[i] += ts->time_step * Bt[i*pinterp+j] * tt * -1.0; 389 b[i] += ts->time_step * B[i*pinterp+j] * tt; 390 } 391 } 392 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"); 393 ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 394 ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 395 ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 396 ierr = PetscFree2(bt,b);CHKERRQ(ierr); 397 PetscFunctionReturn(0); 398 } 399 400 /*------------------------------------------------------------*/ 401 #undef __FUNCT__ 402 #define __FUNCT__ "TSReset_ARKIMEX" 403 static PetscErrorCode TSReset_ARKIMEX(TS ts) 404 { 405 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 406 PetscInt s; 407 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 if (!ark->tableau) PetscFunctionReturn(0); 411 s = ark->tableau->s; 412 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 413 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 414 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 415 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 416 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 417 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 418 ierr = PetscFree(ark->work);CHKERRQ(ierr); 419 PetscFunctionReturn(0); 420 } 421 422 #undef __FUNCT__ 423 #define __FUNCT__ "TSDestroy_ARKIMEX" 424 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 425 { 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 430 ierr = PetscFree(ts->data);CHKERRQ(ierr); 431 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr); 432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr); 433 PetscFunctionReturn(0); 434 } 435 436 /* 437 This defines the nonlinear equation that is to be solved with SNES 438 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 439 */ 440 #undef __FUNCT__ 441 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 442 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 443 { 444 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 445 PetscErrorCode ierr; 446 447 PetscFunctionBegin; 448 ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 449 ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,ark->imex);CHKERRQ(ierr); 450 PetscFunctionReturn(0); 451 } 452 453 #undef __FUNCT__ 454 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 455 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 456 { 457 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 462 ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 463 PetscFunctionReturn(0); 464 } 465 466 #undef __FUNCT__ 467 #define __FUNCT__ "TSSetUp_ARKIMEX" 468 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 469 { 470 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 471 ARKTableau tab = ark->tableau; 472 PetscInt s = tab->s; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 if (!ark->tableau) { 477 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 478 } 479 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 480 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 481 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 482 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 483 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 484 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 485 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 /*------------------------------------------------------------*/ 489 490 #undef __FUNCT__ 491 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 492 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 493 { 494 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 495 PetscErrorCode ierr; 496 char arktype[256]; 497 498 PetscFunctionBegin; 499 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 500 { 501 ARKTableauLink link; 502 PetscInt count,choice; 503 PetscBool flg; 504 const char **namelist; 505 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr); 506 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 507 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 508 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 509 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 510 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 511 ierr = PetscFree(namelist);CHKERRQ(ierr); 512 flg = (PetscBool)!ark->imex; 513 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 514 ark->imex = (PetscBool)!flg; 515 } 516 ierr = PetscOptionsTail();CHKERRQ(ierr); 517 PetscFunctionReturn(0); 518 } 519 520 #undef __FUNCT__ 521 #define __FUNCT__ "PetscFormatRealArray" 522 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 523 { 524 int i,left,count; 525 char *p; 526 527 PetscFunctionBegin; 528 for (i=0,p=buf,left=(int)len; i<n; i++) { 529 count = snprintf(p,left,fmt,x[i]); 530 if (count < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"snprintf()"); 531 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 532 left -= count; 533 p += count; 534 *p++ = ' '; 535 } 536 p[i ? 0 : -1] = 0; 537 PetscFunctionReturn(0); 538 } 539 540 #undef __FUNCT__ 541 #define __FUNCT__ "TSView_ARKIMEX" 542 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 543 { 544 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 545 ARKTableau tab = ark->tableau; 546 PetscBool iascii; 547 PetscErrorCode ierr; 548 549 PetscFunctionBegin; 550 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 551 if (iascii) { 552 const TSARKIMEXType arktype; 553 char buf[512]; 554 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 555 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 556 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 557 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 558 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 559 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 560 } 561 PetscFunctionReturn(0); 562 } 563 564 #undef __FUNCT__ 565 #define __FUNCT__ "TSARKIMEXSetType" 566 /*@C 567 TSARKIMEXSetType - Set the type of ARK IMEX scheme 568 569 Logically collective 570 571 Input Parameter: 572 + ts - timestepping context 573 - arktype - type of ARK-IMEX scheme 574 575 Level: intermediate 576 577 .seealso: TSARKIMEXGetType() 578 @*/ 579 PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 580 { 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 585 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 586 PetscFunctionReturn(0); 587 } 588 589 #undef __FUNCT__ 590 #define __FUNCT__ "TSARKIMEXGetType" 591 /*@C 592 TSARKIMEXGetType - Get the type of ARK IMEX scheme 593 594 Logically collective 595 596 Input Parameter: 597 . ts - timestepping context 598 599 Output Parameter: 600 . arktype - type of ARK-IMEX scheme 601 602 Level: intermediate 603 604 .seealso: TSARKIMEXGetType() 605 @*/ 606 PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 607 { 608 PetscErrorCode ierr; 609 610 PetscFunctionBegin; 611 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 612 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 613 PetscFunctionReturn(0); 614 } 615 616 #undef __FUNCT__ 617 #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 618 /*@C 619 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 620 621 Logically collective 622 623 Input Parameter: 624 + ts - timestepping context 625 - flg - PETSC_TRUE for fully implicit 626 627 Level: intermediate 628 629 .seealso: TSARKIMEXGetType() 630 @*/ 631 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 632 { 633 PetscErrorCode ierr; 634 635 PetscFunctionBegin; 636 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 637 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 638 PetscFunctionReturn(0); 639 } 640 641 EXTERN_C_BEGIN 642 #undef __FUNCT__ 643 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 644 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 645 { 646 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 647 PetscErrorCode ierr; 648 649 PetscFunctionBegin; 650 if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 651 *arktype = ark->tableau->name; 652 PetscFunctionReturn(0); 653 } 654 #undef __FUNCT__ 655 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 656 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 657 { 658 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 659 PetscErrorCode ierr; 660 PetscBool match; 661 ARKTableauLink link; 662 663 PetscFunctionBegin; 664 if (ark->tableau) { 665 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 666 if (match) PetscFunctionReturn(0); 667 } 668 for (link = ARKTableauList; link; link=link->next) { 669 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 670 if (match) { 671 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 672 ark->tableau = &link->tab; 673 PetscFunctionReturn(0); 674 } 675 } 676 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 677 PetscFunctionReturn(0); 678 } 679 #undef __FUNCT__ 680 #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 681 PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 682 { 683 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 684 685 PetscFunctionBegin; 686 ark->imex = (PetscBool)!flg; 687 PetscFunctionReturn(0); 688 } 689 EXTERN_C_END 690 691 /* ------------------------------------------------------------ */ 692 /*MC 693 TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 694 695 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 696 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 697 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 698 699 Notes: 700 This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 701 702 Level: beginner 703 704 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXRegister() 705 706 M*/ 707 EXTERN_C_BEGIN 708 #undef __FUNCT__ 709 #define __FUNCT__ "TSCreate_ARKIMEX" 710 PetscErrorCode TSCreate_ARKIMEX(TS ts) 711 { 712 TS_ARKIMEX *th; 713 PetscErrorCode ierr; 714 715 PetscFunctionBegin; 716 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 717 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 718 #endif 719 720 ts->ops->reset = TSReset_ARKIMEX; 721 ts->ops->destroy = TSDestroy_ARKIMEX; 722 ts->ops->view = TSView_ARKIMEX; 723 ts->ops->setup = TSSetUp_ARKIMEX; 724 ts->ops->step = TSStep_ARKIMEX; 725 ts->ops->interpolate = TSInterpolate_ARKIMEX; 726 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 727 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 728 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 729 730 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 731 ts->data = (void*)th; 732 th->imex = PETSC_TRUE; 733 734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 735 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 736 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 737 PetscFunctionReturn(0); 738 } 739 EXTERN_C_END 740