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