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