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 TSRosWType TSRosWDefault = TSROSW2E; 15 static PetscBool TSRosWRegisterAllCalled; 16 static PetscBool TSRosWPackageInitialized; 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_RosW; 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "TSRosWRegisterAll" 51 /*@C 52 TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW 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, TSRosW, register, all 59 60 .seealso: TSRosWRegisterDestroy() 61 @*/ 62 PetscErrorCode TSRosWRegisterAll(void) 63 { 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 68 TSRosWRegisterAllCalled = 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 = TSRosWRegister(TSROSW2D,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 = PetscSqrtReal((PetscReal)2.0), 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 = TSRosWRegister(TSROSW2E,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 = TSRosWRegister(TSROSW3,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 = TSRosWRegister(TSROSW4,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 = TSRosWRegister(TSROSW5,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__ "TSRosWRegisterDestroy" 164 /*@C 165 TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 166 167 Not Collective 168 169 Level: advanced 170 171 .keywords: TSRosW, register, destroy 172 .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 173 @*/ 174 PetscErrorCode TSRosWRegisterDestroy(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 TSRosWRegisterAllCalled = PETSC_FALSE; 189 PetscFunctionReturn(0); 190 } 191 192 #undef __FUNCT__ 193 #define __FUNCT__ "TSRosWInitializePackage" 194 /*@C 195 TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 196 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 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, TSRosW, initialize, package 205 .seealso: PetscInitialize() 206 @*/ 207 PetscErrorCode TSRosWInitializePackage(const char path[]) 208 { 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 if (TSRosWPackageInitialized) PetscFunctionReturn(0); 213 TSRosWPackageInitialized = PETSC_TRUE; 214 ierr = TSRosWRegisterAll();CHKERRQ(ierr); 215 ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "TSRosWFinalizePackage" 221 /*@C 222 TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 223 called from PetscFinalize(). 224 225 Level: developer 226 227 .keywords: Petsc, destroy, package 228 .seealso: PetscFinalize() 229 @*/ 230 PetscErrorCode TSRosWFinalizePackage(void) 231 { 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 TSRosWPackageInitialized = PETSC_FALSE; 236 ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "TSRosWRegister" 242 /*@C 243 TSRosWRegister - 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: TSRosW 269 @*/ 270 PetscErrorCode TSRosWRegister(const TSRosWType 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_RosW" 309 static PetscErrorCode TSStep_RosW(TS ts) 310 { 311 TS_RosW *ark = (TS_RosW*)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_RosW" 373 static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X) 374 { 375 TS_RosW *ark = (TS_RosW*)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,"TSRosW %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_RosW" 403 static PetscErrorCode TSReset_RosW(TS ts) 404 { 405 TS_RosW *ark = (TS_RosW*)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_RosW" 424 static PetscErrorCode TSDestroy_RosW(TS ts) 425 { 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 430 ierr = PetscFree(ts->data);CHKERRQ(ierr); 431 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_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_RosW" 442 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts) 443 { 444 TS_RosW *ark = (TS_RosW*)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_RosW" 455 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 456 { 457 TS_RosW *ark = (TS_RosW*)ts->data; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 /* ark->Ydot has already been computed in SNESTSFormFunction_RosW (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_RosW" 468 static PetscErrorCode TSSetUp_RosW(TS ts) 469 { 470 TS_RosW *ark = (TS_RosW*)ts->data; 471 ARKTableau tab = ark->tableau; 472 PetscInt s = tab->s; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 if (!ark->tableau) { 477 ierr = TSRosWSetType(ts,TSRosWDefault);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_RosW" 492 static PetscErrorCode TSSetFromOptions_RosW(TS ts) 493 { 494 TS_RosW *ark = (TS_RosW*)ts->data; 495 PetscErrorCode ierr; 496 char arktype[256]; 497 498 PetscFunctionBegin; 499 ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 500 { 501 ARKTableauLink link; 502 PetscInt count,choice; 503 PetscBool flg; 504 const char **namelist; 505 ierr = PetscStrncpy(arktype,TSRosWDefault,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_rosw_type","Family of ARK IMEX method","TSRosWSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 510 ierr = TSRosWSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 511 ierr = PetscFree(namelist);CHKERRQ(ierr); 512 flg = (PetscBool)!ark->imex; 513 ierr = PetscOptionsBool("-ts_rosw_fully_implicit","Solve the problem fully implicitly","TSRosWSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 514 ark->imex = (PetscBool)!flg; 515 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 516 } 517 ierr = PetscOptionsTail();CHKERRQ(ierr); 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNCT__ 522 #define __FUNCT__ "PetscFormatRealArray" 523 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 524 { 525 PetscErrorCode ierr; 526 int i,left,count; 527 char *p; 528 529 PetscFunctionBegin; 530 for (i=0,p=buf,left=(int)len; i<n; i++) { 531 ierr = PetscSNPrintf(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 532 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 533 left -= count; 534 p += count; 535 *p++ = ' '; 536 } 537 p[i ? 0 : -1] = 0; 538 PetscFunctionReturn(0); 539 } 540 541 #undef __FUNCT__ 542 #define __FUNCT__ "TSView_RosW" 543 static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 544 { 545 TS_RosW *ark = (TS_RosW*)ts->data; 546 ARKTableau tab = ark->tableau; 547 PetscBool iascii; 548 PetscErrorCode ierr; 549 550 PetscFunctionBegin; 551 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 552 if (iascii) { 553 const TSRosWType arktype; 554 char buf[512]; 555 ierr = TSRosWGetType(ts,&arktype);CHKERRQ(ierr); 556 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 557 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 558 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 559 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 560 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 561 } 562 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "TSRosWSetType" 568 /*@C 569 TSRosWSetType - Set the type of ARK IMEX scheme 570 571 Logically collective 572 573 Input Parameter: 574 + ts - timestepping context 575 - arktype - type of ARK-IMEX scheme 576 577 Level: intermediate 578 579 .seealso: TSRosWGetType() 580 @*/ 581 PetscErrorCode TSRosWSetType(TS ts,const TSRosWType arktype) 582 { 583 PetscErrorCode ierr; 584 585 PetscFunctionBegin; 586 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 587 ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,arktype));CHKERRQ(ierr); 588 PetscFunctionReturn(0); 589 } 590 591 #undef __FUNCT__ 592 #define __FUNCT__ "TSRosWGetType" 593 /*@C 594 TSRosWGetType - Get the type of ARK IMEX scheme 595 596 Logically collective 597 598 Input Parameter: 599 . ts - timestepping context 600 601 Output Parameter: 602 . arktype - type of ARK-IMEX scheme 603 604 Level: intermediate 605 606 .seealso: TSRosWGetType() 607 @*/ 608 PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *arktype) 609 { 610 PetscErrorCode ierr; 611 612 PetscFunctionBegin; 613 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 614 ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,arktype));CHKERRQ(ierr); 615 PetscFunctionReturn(0); 616 } 617 618 #undef __FUNCT__ 619 #define __FUNCT__ "TSRosWSetFullyImplicit" 620 /*@C 621 TSRosWSetFullyImplicit - Solve both parts of the equation implicitly 622 623 Logically collective 624 625 Input Parameter: 626 + ts - timestepping context 627 - flg - PETSC_TRUE for fully implicit 628 629 Level: intermediate 630 631 .seealso: TSRosWGetType() 632 @*/ 633 PetscErrorCode TSRosWSetFullyImplicit(TS ts,PetscBool flg) 634 { 635 PetscErrorCode ierr; 636 637 PetscFunctionBegin; 638 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 639 ierr = PetscTryMethod(ts,"TSRosWSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 EXTERN_C_BEGIN 644 #undef __FUNCT__ 645 #define __FUNCT__ "TSRosWGetType_RosW" 646 PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *arktype) 647 { 648 TS_RosW *ark = (TS_RosW*)ts->data; 649 PetscErrorCode ierr; 650 651 PetscFunctionBegin; 652 if (!ark->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 653 *arktype = ark->tableau->name; 654 PetscFunctionReturn(0); 655 } 656 #undef __FUNCT__ 657 #define __FUNCT__ "TSRosWSetType_RosW" 658 PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType arktype) 659 { 660 TS_RosW *ark = (TS_RosW*)ts->data; 661 PetscErrorCode ierr; 662 PetscBool match; 663 ARKTableauLink link; 664 665 PetscFunctionBegin; 666 if (ark->tableau) { 667 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 668 if (match) PetscFunctionReturn(0); 669 } 670 for (link = ARKTableauList; link; link=link->next) { 671 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 672 if (match) { 673 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 674 ark->tableau = &link->tab; 675 PetscFunctionReturn(0); 676 } 677 } 678 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 679 PetscFunctionReturn(0); 680 } 681 #undef __FUNCT__ 682 #define __FUNCT__ "TSRosWSetFullyImplicit_RosW" 683 PetscErrorCode TSRosWSetFullyImplicit_RosW(TS ts,PetscBool flg) 684 { 685 TS_RosW *ark = (TS_RosW*)ts->data; 686 687 PetscFunctionBegin; 688 ark->imex = (PetscBool)!flg; 689 PetscFunctionReturn(0); 690 } 691 EXTERN_C_END 692 693 /* ------------------------------------------------------------ */ 694 /*MC 695 TSRosW - ODE solver using Rosenbrock-W schemes 696 697 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 698 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 699 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 700 701 Notes: 702 This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 703 704 Level: beginner 705 706 .seealso: TSCreate(), TS, TSSetType(), TSRosWRegister() 707 708 M*/ 709 EXTERN_C_BEGIN 710 #undef __FUNCT__ 711 #define __FUNCT__ "TSCreate_RosW" 712 PetscErrorCode TSCreate_RosW(TS ts) 713 { 714 TS_RosW *th; 715 PetscErrorCode ierr; 716 717 PetscFunctionBegin; 718 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 719 ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 720 #endif 721 722 ts->ops->reset = TSReset_RosW; 723 ts->ops->destroy = TSDestroy_RosW; 724 ts->ops->view = TSView_RosW; 725 ts->ops->setup = TSSetUp_RosW; 726 ts->ops->step = TSStep_RosW; 727 ts->ops->interpolate = TSInterpolate_RosW; 728 ts->ops->setfromoptions = TSSetFromOptions_RosW; 729 ts->ops->snesfunction = SNESTSFormFunction_RosW; 730 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 731 732 ierr = PetscNewLog(ts,TS_RosW,&th);CHKERRQ(ierr); 733 ts->data = (void*)th; 734 th->imex = PETSC_TRUE; 735 736 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 738 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetFullyImplicit_C","TSRosWSetFullyImplicit_RosW",TSRosWSetFullyImplicit_RosW);CHKERRQ(ierr); 739 PetscFunctionReturn(0); 740 } 741 EXTERN_C_END 742