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 next_time_step; 320 PetscReal h,t; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 325 next_time_step = ts->time_step; 326 h = ts->time_step; 327 t = ts->ptime; 328 329 for (i=0; i<s; i++) { 330 if (At[i*s+i] == 0) { /* This stage is explicit */ 331 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 332 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 333 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 334 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 335 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 336 } else { 337 ark->stage_time = t + h*ct[i]; 338 ark->shift = 1./(h*At[i*s+i]); 339 /* Affine part */ 340 ierr = VecZeroEntries(W);CHKERRQ(ierr); 341 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 342 ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 343 /* Ydot = shift*(Y-Z) */ 344 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 345 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 346 ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr); 347 /* Initial guess taken from last stage */ 348 ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 349 ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 350 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 351 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 352 ts->nonlinear_its += its; ts->linear_its += lits; 353 } 354 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 355 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr); 356 if (ark->imex) { 357 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 358 } else { 359 ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr); 360 } 361 } 362 for (j=0; j<s; j++) w[j] = -h*bt[j]; 363 ierr = VecMAXPY(ts->vec_sol,s,w,YdotI);CHKERRQ(ierr); 364 for (j=0; j<s; j++) w[j] = h*b[j]; 365 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 366 367 ts->ptime += ts->time_step; 368 ts->time_step = next_time_step; 369 ts->steps++; 370 PetscFunctionReturn(0); 371 } 372 373 #undef __FUNCT__ 374 #define __FUNCT__ "TSInterpolate_RosW" 375 static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X) 376 { 377 TS_RosW *ark = (TS_RosW*)ts->data; 378 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 379 PetscReal tt,t = (itime - ts->ptime)/ts->time_step + 1; /* In the interval [0,1] */ 380 PetscScalar *bt,*b; 381 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 382 PetscErrorCode ierr; 383 384 PetscFunctionBegin; 385 if (!Bt || !B) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ark->tableau->name); 386 ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 387 for (i=0; i<s; i++) bt[i] = b[i] = 0; 388 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 389 for (i=0; i<s; i++) { 390 bt[i] += ts->time_step * Bt[i*pinterp+j] * tt * -1.0; 391 b[i] += ts->time_step * B[i*pinterp+j] * tt; 392 } 393 } 394 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"); 395 ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 396 ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 397 ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 398 ierr = PetscFree2(bt,b);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 /*------------------------------------------------------------*/ 403 #undef __FUNCT__ 404 #define __FUNCT__ "TSReset_RosW" 405 static PetscErrorCode TSReset_RosW(TS ts) 406 { 407 TS_RosW *ark = (TS_RosW*)ts->data; 408 PetscInt s; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 if (!ark->tableau) PetscFunctionReturn(0); 413 s = ark->tableau->s; 414 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 415 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 416 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 417 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 418 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 419 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 420 ierr = PetscFree(ark->work);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 #undef __FUNCT__ 425 #define __FUNCT__ "TSDestroy_RosW" 426 static PetscErrorCode TSDestroy_RosW(TS ts) 427 { 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 432 ierr = PetscFree(ts->data);CHKERRQ(ierr); 433 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 434 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 /* 439 This defines the nonlinear equation that is to be solved with SNES 440 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 441 */ 442 #undef __FUNCT__ 443 #define __FUNCT__ "SNESTSFormFunction_RosW" 444 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts) 445 { 446 TS_RosW *ark = (TS_RosW*)ts->data; 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 451 ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,ark->imex);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 455 #undef __FUNCT__ 456 #define __FUNCT__ "SNESTSFormJacobian_RosW" 457 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 458 { 459 TS_RosW *ark = (TS_RosW*)ts->data; 460 PetscErrorCode ierr; 461 462 PetscFunctionBegin; 463 /* ark->Ydot has already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 464 ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 465 PetscFunctionReturn(0); 466 } 467 468 #undef __FUNCT__ 469 #define __FUNCT__ "TSSetUp_RosW" 470 static PetscErrorCode TSSetUp_RosW(TS ts) 471 { 472 TS_RosW *ark = (TS_RosW*)ts->data; 473 ARKTableau tab = ark->tableau; 474 PetscInt s = tab->s; 475 PetscErrorCode ierr; 476 477 PetscFunctionBegin; 478 if (!ark->tableau) { 479 ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 480 } 481 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 482 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 483 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 484 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 485 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 486 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 487 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 488 PetscFunctionReturn(0); 489 } 490 /*------------------------------------------------------------*/ 491 492 #undef __FUNCT__ 493 #define __FUNCT__ "TSSetFromOptions_RosW" 494 static PetscErrorCode TSSetFromOptions_RosW(TS ts) 495 { 496 TS_RosW *ark = (TS_RosW*)ts->data; 497 PetscErrorCode ierr; 498 char arktype[256]; 499 500 PetscFunctionBegin; 501 ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 502 { 503 ARKTableauLink link; 504 PetscInt count,choice; 505 PetscBool flg; 506 const char **namelist; 507 ierr = PetscStrncpy(arktype,TSRosWDefault,sizeof arktype);CHKERRQ(ierr); 508 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 509 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 510 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 511 ierr = PetscOptionsEList("-ts_rosw_type","Family of ARK IMEX method","TSRosWSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 512 ierr = TSRosWSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 513 ierr = PetscFree(namelist);CHKERRQ(ierr); 514 flg = (PetscBool)!ark->imex; 515 ierr = PetscOptionsBool("-ts_rosw_fully_implicit","Solve the problem fully implicitly","TSRosWSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 516 ark->imex = (PetscBool)!flg; 517 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 518 } 519 ierr = PetscOptionsTail();CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 523 #undef __FUNCT__ 524 #define __FUNCT__ "PetscFormatRealArray" 525 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 526 { 527 PetscErrorCode ierr; 528 int i,left,count; 529 char *p; 530 531 PetscFunctionBegin; 532 for (i=0,p=buf,left=(int)len; i<n; i++) { 533 ierr = PetscSNPrintf(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 534 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 535 left -= count; 536 p += count; 537 *p++ = ' '; 538 } 539 p[i ? 0 : -1] = 0; 540 PetscFunctionReturn(0); 541 } 542 543 #undef __FUNCT__ 544 #define __FUNCT__ "TSView_RosW" 545 static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 546 { 547 TS_RosW *ark = (TS_RosW*)ts->data; 548 ARKTableau tab = ark->tableau; 549 PetscBool iascii; 550 PetscErrorCode ierr; 551 552 PetscFunctionBegin; 553 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 554 if (iascii) { 555 const TSRosWType arktype; 556 char buf[512]; 557 ierr = TSRosWGetType(ts,&arktype);CHKERRQ(ierr); 558 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 559 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 560 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 561 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 562 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 563 } 564 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 568 #undef __FUNCT__ 569 #define __FUNCT__ "TSRosWSetType" 570 /*@C 571 TSRosWSetType - Set the type of ARK IMEX scheme 572 573 Logically collective 574 575 Input Parameter: 576 + ts - timestepping context 577 - arktype - type of ARK-IMEX scheme 578 579 Level: intermediate 580 581 .seealso: TSRosWGetType() 582 @*/ 583 PetscErrorCode TSRosWSetType(TS ts,const TSRosWType arktype) 584 { 585 PetscErrorCode ierr; 586 587 PetscFunctionBegin; 588 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 589 ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,arktype));CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 593 #undef __FUNCT__ 594 #define __FUNCT__ "TSRosWGetType" 595 /*@C 596 TSRosWGetType - Get the type of ARK IMEX scheme 597 598 Logically collective 599 600 Input Parameter: 601 . ts - timestepping context 602 603 Output Parameter: 604 . arktype - type of ARK-IMEX scheme 605 606 Level: intermediate 607 608 .seealso: TSRosWGetType() 609 @*/ 610 PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *arktype) 611 { 612 PetscErrorCode ierr; 613 614 PetscFunctionBegin; 615 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 616 ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,arktype));CHKERRQ(ierr); 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "TSRosWSetFullyImplicit" 622 /*@C 623 TSRosWSetFullyImplicit - Solve both parts of the equation implicitly 624 625 Logically collective 626 627 Input Parameter: 628 + ts - timestepping context 629 - flg - PETSC_TRUE for fully implicit 630 631 Level: intermediate 632 633 .seealso: TSRosWGetType() 634 @*/ 635 PetscErrorCode TSRosWSetFullyImplicit(TS ts,PetscBool flg) 636 { 637 PetscErrorCode ierr; 638 639 PetscFunctionBegin; 640 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 641 ierr = PetscTryMethod(ts,"TSRosWSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 642 PetscFunctionReturn(0); 643 } 644 645 EXTERN_C_BEGIN 646 #undef __FUNCT__ 647 #define __FUNCT__ "TSRosWGetType_RosW" 648 PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *arktype) 649 { 650 TS_RosW *ark = (TS_RosW*)ts->data; 651 PetscErrorCode ierr; 652 653 PetscFunctionBegin; 654 if (!ark->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 655 *arktype = ark->tableau->name; 656 PetscFunctionReturn(0); 657 } 658 #undef __FUNCT__ 659 #define __FUNCT__ "TSRosWSetType_RosW" 660 PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType arktype) 661 { 662 TS_RosW *ark = (TS_RosW*)ts->data; 663 PetscErrorCode ierr; 664 PetscBool match; 665 ARKTableauLink link; 666 667 PetscFunctionBegin; 668 if (ark->tableau) { 669 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 670 if (match) PetscFunctionReturn(0); 671 } 672 for (link = ARKTableauList; link; link=link->next) { 673 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 674 if (match) { 675 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 676 ark->tableau = &link->tab; 677 PetscFunctionReturn(0); 678 } 679 } 680 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 681 PetscFunctionReturn(0); 682 } 683 #undef __FUNCT__ 684 #define __FUNCT__ "TSRosWSetFullyImplicit_RosW" 685 PetscErrorCode TSRosWSetFullyImplicit_RosW(TS ts,PetscBool flg) 686 { 687 TS_RosW *ark = (TS_RosW*)ts->data; 688 689 PetscFunctionBegin; 690 ark->imex = (PetscBool)!flg; 691 PetscFunctionReturn(0); 692 } 693 EXTERN_C_END 694 695 /* ------------------------------------------------------------ */ 696 /*MC 697 TSRosW - ODE solver using Rosenbrock-W schemes 698 699 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 700 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 701 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 702 703 Notes: 704 This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 705 706 Level: beginner 707 708 .seealso: TSCreate(), TS, TSSetType(), TSRosWRegister() 709 710 M*/ 711 EXTERN_C_BEGIN 712 #undef __FUNCT__ 713 #define __FUNCT__ "TSCreate_RosW" 714 PetscErrorCode TSCreate_RosW(TS ts) 715 { 716 TS_RosW *th; 717 PetscErrorCode ierr; 718 719 PetscFunctionBegin; 720 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 721 ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 722 #endif 723 724 ts->ops->reset = TSReset_RosW; 725 ts->ops->destroy = TSDestroy_RosW; 726 ts->ops->view = TSView_RosW; 727 ts->ops->setup = TSSetUp_RosW; 728 ts->ops->step = TSStep_RosW; 729 ts->ops->interpolate = TSInterpolate_RosW; 730 ts->ops->setfromoptions = TSSetFromOptions_RosW; 731 ts->ops->snesfunction = SNESTSFormFunction_RosW; 732 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 733 734 ierr = PetscNewLog(ts,TS_RosW,&th);CHKERRQ(ierr); 735 ts->data = (void*)th; 736 th->imex = PETSC_TRUE; 737 738 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 739 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetFullyImplicit_C","TSRosWSetFullyImplicit_RosW",TSRosWSetFullyImplicit_RosW);CHKERRQ(ierr); 741 PetscFunctionReturn(0); 742 } 743 EXTERN_C_END 744