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