00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
00037 #define OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
00038
00039 #include <cassert>
00040 #include <cstddef>
00041
00042 #include <algorithm>
00043 #include <vector>
00044 #include <iostream>
00045
00046 namespace Opm {
00047 namespace polymer_reorder {
00048 class ModelParameterStorage {
00049 public:
00050 ModelParameterStorage(int nc, int totconn)
00051 : drho_(0.0), rockdensity_(0.0), mob_(0),
00052 dmobds_(0), dmobwatdc_(0), mc_(0),
00053 dmcdc_(0), porevol_(0), porosity_(0), dg_(0), sw_(0), c_(0), cmax_(0),
00054 ds_(0), dsc_(0), dcads_(0), dcadsdc_(0), pc_(0), dpc_(0),
00055 trans_(0), data_()
00056 {
00057 size_t alloc_sz;
00058
00059 alloc_sz = 2 * nc;
00060 alloc_sz += 2 * nc;
00061 alloc_sz += nc;
00062 alloc_sz += nc;
00063 alloc_sz += nc;
00064 alloc_sz += 1 * nc;
00065 alloc_sz += 1 * nc;
00066 alloc_sz += 1 * totconn;
00067 alloc_sz += 1 * nc;
00068 alloc_sz += 1 * nc;
00069 alloc_sz += 1 * nc;
00070 alloc_sz += 1 * nc;
00071 alloc_sz += 1 * nc;
00072 alloc_sz += 1 * nc;
00073 alloc_sz += 1 * nc;
00074 alloc_sz += 1 * nc;
00075 alloc_sz += 1 * nc;
00076 alloc_sz += 1 * totconn;
00077 data_.resize(alloc_sz);
00078
00079 mob_ = &data_[0];
00080 dmobds_ = mob_ + (2 * nc );
00081 dmobwatdc_ = dmobds_ + (2 * nc );
00082 mc_ = dmobwatdc_ + (1 * nc );
00083 dmcdc_ = mc_ + (1 * nc );
00084 porevol_ = dmcdc_ + (1 * nc );
00085 porosity_ = porevol_ + (1 * nc );
00086 dg_ = porosity_ + (1 * nc );
00087 sw_ = dg_ + (1 * totconn);
00088 c_ = sw_ + (1 * nc );
00089 cmax_ = c_ + (1 * nc );
00090 ds_ = cmax_ + (1 * nc );
00091 dsc_ = ds_ + (1 * nc );
00092 dcads_ = dsc_ + (1 * nc );
00093 dcadsdc_ = dcads_ + (1 * nc );
00094 pc_ = dcadsdc_ + (1 * nc );
00095 dpc_ = pc_ + (1 * nc );
00096 trans_ = dpc_ + (1 * nc );
00097 }
00098
00099 double& drho () { return drho_ ; }
00100 double drho () const { return drho_ ; }
00101
00102 double& rockdensity() { return rockdensity_ ; }
00103 double rockdensity() const { return rockdensity_ ; }
00104
00105 double* mob (int cell) { return mob_ + (2*cell + 0); }
00106 const double* mob (int cell) const { return mob_ + (2*cell + 0); }
00107
00108 double* dmobds (int cell) { return dmobds_ + (2*cell + 0); }
00109 const double* dmobds (int cell) const { return dmobds_ + (2*cell + 0); }
00110
00111 double& dmobwatdc (int cell) { return dmobwatdc_[cell]; }
00112 double dmobwatdc (int cell) const { return dmobwatdc_[cell]; }
00113
00114 double& mc (int cell) { return mc_[cell]; }
00115 double mc (int cell) const { return mc_[cell]; }
00116
00117 double& dmcdc (int cell) { return dmcdc_[cell]; }
00118 double dmcdc (int cell) const { return dmcdc_[cell]; }
00119
00120 double* porevol() { return porevol_ ; }
00121 double porevol(int cell) const { return porevol_[cell] ; }
00122
00123 double* porosity() { return porosity_ ; }
00124 double porosity(int cell) const { return porosity_[cell] ; }
00125
00126
00127 double& dg(int i) { return dg_[i] ; }
00128 double dg(int i) const { return dg_[i] ; }
00129
00130 double& sw(int cell) { return sw_[cell] ; }
00131 double sw(int cell) const { return sw_[cell] ; }
00132
00133 double& c(int cell) { return c_[cell] ; }
00134 double c(int cell) const { return c_[cell] ; }
00135
00136 double& cmax(int cell) { return cmax_[cell] ; }
00137 double cmax(int cell) const { return cmax_[cell] ; }
00138
00139 double& ds(int cell) { return ds_[cell] ; }
00140 double ds(int cell) const { return ds_[cell] ; }
00141
00142 double& dsc(int cell) { return dsc_[cell] ; }
00143 double dsc(int cell) const { return dsc_[cell] ; }
00144
00145 double& dcads(int cell) { return dcads_[cell] ; }
00146 double dcads(int cell) const { return dcads_[cell] ; }
00147
00148 double& dcadsdc(int cell) { return dcadsdc_[cell] ; }
00149 double dcadsdc(int cell) const { return dcadsdc_[cell] ; }
00150
00151 double& pc(int cell) { return pc_[cell] ; }
00152 double pc(int cell) const { return pc_[cell] ; }
00153
00154 double& dpc(int cell) { return dpc_[cell] ; }
00155 double dpc(int cell) const { return dpc_[cell] ; }
00156
00157 double& trans(int f) { return trans_[f] ; }
00158 double trans(int f) const { return trans_[f] ; }
00159
00160 private:
00161 double drho_ ;
00162 double rockdensity_ ;
00163 double *mob_ ;
00164 double *dmobds_ ;
00165 double *dmobwatdc_ ;
00166 double *mc_ ;
00167 double *dmcdc_ ;
00168 double *porevol_ ;
00169 double *porosity_ ;
00170 double *dg_ ;
00171 double *sw_ ;
00172 double *c_ ;
00173 double *cmax_ ;
00174 double *ds_ ;
00175 double *dsc_ ;
00176 double *dcads_ ;
00177 double *dcadsdc_ ;
00178 double *pc_ ;
00179 double *dpc_ ;
00180 double *trans_ ;
00181
00182 std::vector<double> data_;
00183 };
00184 }
00185
00186
00187 template <class TwophaseFluidPolymer>
00188 class SinglePointUpwindTwoPhasePolymer {
00189 public:
00190 template <class Grid>
00191 SinglePointUpwindTwoPhasePolymer(const TwophaseFluidPolymer& fluid ,
00192 const Grid& g ,
00193 const std::vector<double>& porevol ,
00194 const double* grav = 0,
00195 const bool guess_previous = true)
00196 : fluid_ (fluid) ,
00197 gravity_ (grav) ,
00198 f2hf_ (2 * g.number_of_faces, -1) ,
00199 store_ (g.number_of_cells,
00200 g.cell_facepos[ g.number_of_cells ]),
00201 init_step_use_previous_sol_(guess_previous) ,
00202 sat_tol_ (1e-5)
00203 {
00204
00205 if (gravity_) {
00206 store_.drho() = fluid_.density(0) - fluid_.density(1);
00207 }
00208
00209 for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
00210 for (; i < g.cell_facepos[c + 1]; ++i) {
00211 const int f = g.cell_faces[i];
00212 const int p = 1 - (g.face_cells[2*f + 0] == c);
00213 f2hf_[2*f + p] = i;
00214 }
00215 }
00216
00217 std::copy(porevol.begin(), porevol.end(), store_.porevol());
00218 const double* poro = fluid.porosity();
00219 std::copy(poro, poro + g.number_of_cells, store_.porosity());
00220 store_.rockdensity() = fluid.rockdensity();
00221 }
00222
00223 void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces,
00224 const std::vector<int>& nb_faces)
00225 {
00226 std::vector<int> nbhf(hf_faces.size());
00227 for(unsigned int i=0; i<p_faces.size(); ++i){
00228 int nbf = nb_faces[i];
00229 if(f2hf_[2*nbf] == -1){
00230 nbhf[i] = f2hf_[2*nbf+1];
00231 }else{
00232 assert(f2hf_[2*nbf+1]==-1);
00233 nbhf[i] = f2hf_[2*nbf];
00234 }
00235 }
00236 for(unsigned int i=0; i<p_faces.size(); ++i){
00237
00238 int f = p_faces[i];
00239 int hf = hf_faces[i];
00240 bool changed=false;
00241
00242 if(f2hf_[2*f] == hf){
00243 assert(f2hf_[2*f+1]==-1);
00244 }else{
00245 assert(f2hf_[2*f]==-1);
00246 f2hf_[2*f]=nbhf[i];
00247 changed=true;
00248 }
00249 if(!changed){
00250 if(f2hf_[2*f+1]== hf){
00251 assert(f2hf_[2*f]==-1);
00252 }else{
00253 assert(f2hf_[2*f+1]==-1);
00254 f2hf_[2*f+1]=nbhf[i];
00255 changed=true;
00256 }
00257 }
00258 assert(changed);
00259 }
00260 }
00261
00262
00263
00264
00265
00266 enum { DofPerCell = 1 };
00267
00268 void
00269 initResidual(const int c, double* Fs, double* Fc) const {
00270 (void) c;
00271 *Fs = 0.0;
00272 *Fc = 0.0;
00273 }
00274
00275 template <class ReservoirState,
00276 class Grid >
00277 void
00278 fluxConnection(const ReservoirState& state ,
00279 const Grid& g ,
00280 const double dt ,
00281 const int cell ,
00282 const int f ,
00283 double* F ,
00284 double* dFd1 ,
00285 double* dFd2
00286
00287
00288
00289 ) const {
00290
00291 const int *n = g.face_cells + (2 * f);
00292 double dflux = state.faceflux()[f];
00293 double gflux = gravityFlux(f);
00294 double pcflux, dpcflux[2];
00295 capFlux(f, n, pcflux, dpcflux);
00296 gflux += pcflux;
00297
00298 int pix[2];
00299 double m[2], dmds[2], dmobwatdc;
00300 double mc, dmcdc;
00301 upwindMobility(dflux, gflux, n, pix, m, dmds, dmobwatdc, mc, dmcdc);
00302
00303 assert ((m[0] >= 0.0) && (m[1] >= 0.0));
00304
00305 double mt = m[0] + m[1];
00306 assert (mt >= 0.0);
00307
00308 double sgn = 2.0*(n[0] == cell) - 1.0;
00309 dflux *= sgn;
00310 gflux *= sgn;
00311
00312
00313 double f1 = m[0] / mt;
00314 const double v1 = dflux + m[1]*gflux;
00315
00316
00317 F[0] += dt * f1 * v1;
00318 F[1] += dt * mc * f1 * v1;
00319
00320
00321 double *dFsds[2];
00322 double *dFsdc[2];
00323 double *dFcds[2];
00324 double *dFcdc[2];
00325 if (n[0] == cell) {
00326 dFsds[0] = &dFd1[0]; dFsds[1] = &dFd2[0];
00327 dFsdc[0] = &dFd1[1]; dFsdc[1] = &dFd2[1];
00328 dFcds[0] = &dFd1[2]; dFcds[1] = &dFd2[2];
00329 dFcdc[0] = &dFd1[3]; dFcdc[1] = &dFd2[3];
00330
00331 dFd1[0] += sgn*dt * f1 * dpcflux[0] * m[1];
00332 dFd2[0] += sgn*dt * f1 * dpcflux[1] * m[1];
00333 dFd1[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
00334 dFd2[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
00335
00336
00337 } else {
00338 dFsds[0] = &dFd2[0]; dFsds[1] = &dFd1[0];
00339 dFsdc[0] = &dFd2[1]; dFsdc[1] = &dFd1[1];
00340 dFcds[0] = &dFd2[2]; dFcds[1] = &dFd1[2];
00341 dFcdc[0] = &dFd2[3]; dFcdc[1] = &dFd1[3];
00342
00343 dFd1[0] += sgn*dt * f1 * dpcflux[1] * m[1];
00344 dFd2[0] += sgn*dt * f1 * dpcflux[0] * m[1];
00345 dFd1[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
00346 dFd2[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
00347
00348
00349 }
00350
00351
00352 *dFsds[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmds[0];
00353
00354 *dFcds[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmds[0];
00355
00356
00357 *dFsds[ pix[1] ] -= dt * f1 / mt * v1 * dmds[1];
00358 *dFsds[ pix[1] ] += dt * f1 * gflux * dmds[1];
00359
00360 *dFcds[ pix[1] ] -= dt * f1 / mt * v1 * mc * dmds[1];
00361 *dFcds[ pix[1] ] += dt * f1 * gflux * mc * dmds[1];
00362
00363
00364 *dFsdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmobwatdc;
00365
00366 *dFcdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmobwatdc;
00367 *dFcdc[ pix[0] ] += dt * f1 * v1 * dmcdc;
00368 }
00369
00370 template <class Grid>
00371 void
00372 accumulation(const Grid& g,
00373 const int cell,
00374 double* F,
00375 double* dF
00376 ) const {
00377 (void) g;
00378
00379 const double pv = store_.porevol(cell);
00380 const double dps = fluid_.deadporespace();
00381 const double rhor = fluid_.rockdensity();
00382 const double poro = store_.porosity(cell);
00383
00384 F[0] += pv * store_.ds(cell);
00385 F[1] += pv * (1 - dps) * store_.dsc(cell) + rhor*(1 - poro)/poro*pv*store_.dcads(cell);
00386 dF[0] += pv;
00387 dF[1] += 0.;
00388 dF[2] += pv * (1 - dps) * store_.c(cell);
00389 dF[3] += pv * (1 - dps) * store_.sw(cell) + rhor*(1 - poro)/poro*pv*store_.dcadsdc(cell);
00390 }
00391
00392 template <class Grid ,
00393 class SourceTerms>
00394 void
00395 sourceTerms(const Grid& g ,
00396 const SourceTerms* src,
00397 const int i ,
00398 const double dt ,
00399 double* J ,
00400 double* F ) const {
00401
00402 (void) g;
00403
00404 double dflux = -src->flux[i];
00405
00406 if (dflux < 0) {
00407
00408 *F += dt * dflux * src->saturation[2*i + 0];
00409 } else {
00410
00411 const int cell = src->cell[i];
00412 const double* m = store_.mob (cell);
00413 const double* dm = store_.dmobds(cell);
00414
00415 const double mt = m[0] + m[1];
00416
00417 assert (! ((m[0] < 0) || (m[1] < 0)));
00418 assert (mt > 0);
00419
00420 const double f = m[0] / mt;
00421 const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
00422
00423 *F += dt * dflux * f;
00424 *J += dt * dflux * df;
00425 }
00426 }
00427 template <class Grid>
00428 void
00429 initGravityTrans(const Grid& g ,
00430 const std::vector<double> & htrans) {
00431
00432 assert (htrans.size() ==
00433 static_cast<std::vector<double>::size_type>(g.cell_facepos[ g.number_of_cells ]));
00434
00435 for (int f = 0; f < g.number_of_faces; ++f) {
00436 store_.trans(f) = 0.0;
00437 }
00438
00439 for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
00440 for (; i < g.cell_facepos[c + 1]; ++i) {
00441 int f = g.cell_faces[i];
00442
00443 assert (htrans[i] > 0.0);
00444
00445 store_.trans(f) += 1.0 / htrans[i];
00446 }
00447 }
00448
00449 for (int f = 0; f < g.number_of_faces; ++f) {
00450 store_.trans(f) = 1.0 / store_.trans(f);
00451 }
00452
00453 if (gravity_) {
00454 this->computeStaticGravity(g);
00455 }
00456 }
00457
00458
00459
00460
00461
00462 template <class ReservoirState,
00463 class Grid ,
00464 class JacobianSystem>
00465 void
00466 initStep(const ReservoirState& state,
00467 const Grid& g ,
00468 JacobianSystem& sys ) {
00469
00470 (void) state;
00471
00472 typename JacobianSystem::vector_type& x =
00473 sys.vector().writableSolution();
00474
00475 assert (x.size() == (::std::size_t) (2*g.number_of_cells));
00476
00477 if (init_step_use_previous_sol_) {
00478 std::fill(x.begin(), x.end(), 0.0);
00479 } else {
00480 std::fill(x.begin(), x.end(), 0.0);
00481 const std::vector<double>& s = state.saturation();
00482 const std::vector<double>& c = state.concentration();
00483 for (int cell = 0, ncell = g.number_of_cells; cell < ncell; ++cell) {
00484
00485 x[2*cell + 0] = 0.5 - s[2*cell + 0];
00486 x[2*cell + 1] = 1e-5 - c[cell];
00487 }
00488 }
00489 }
00490
00491 template <class ReservoirState,
00492 class Grid ,
00493 class JacobianSystem>
00494 bool
00495 initIteration(const ReservoirState& state,
00496 const Grid& g ,
00497 JacobianSystem& sys) {
00498
00499 double s[2];
00500 double mob[2];
00501 double dmobds[4];
00502 double dmobwatdc;
00503 double c, cmax;
00504 double mc, dmcdc;
00505 double pc, dpc;
00506
00507 const typename JacobianSystem::vector_type& x =
00508 sys.vector().solution();
00509 const ::std::vector<double>& sat = state.saturation();
00510 const ::std::vector<double>& cpoly = state.concentration();
00511 const ::std::vector<double>& cmaxpoly = state.maxconcentration();
00512
00513 bool in_range = true;
00514 for (int cell = 0; cell < g.number_of_cells; ++cell) {
00515
00516 store_.ds(cell) = x[2*cell + 0];
00517 s[0] = sat[cell*2 + 0] + x[2*cell + 0];
00518 c = cpoly[cell] + x[2*cell + 1];
00519 store_.sw(cell) = s[0];
00520 store_.c(cell) = c;
00521 cmax = std::max(c, cmaxpoly[cell]);
00522 store_.cmax(cell) = cmax;
00523 store_.dsc(cell) = s[0]*c - sat[cell*2 + 0]*cpoly[cell];
00524 double dcadsdc;
00525 double cads;
00526 fluid_.adsorption(cpoly[cell], cmaxpoly[cell], cads, dcadsdc);
00527 store_.dcads(cell) = -cads;
00528 fluid_.adsorption(c, cmax, cads, dcadsdc);
00529 store_.dcads(cell) += cads;
00530 store_.dcadsdc(cell) = dcadsdc;
00531 double s_min = fluid_.s_min(cell);
00532 double s_max = fluid_.s_max(cell);
00533
00534 if ( s[0] < (s_min - sat_tol_) || s[0] > (s_max + sat_tol_) ) {
00535
00536
00537
00538
00539
00540
00541 in_range = false;
00542 }
00543 s[0] = std::max(s_min, s[0]);
00544 s[0] = std::min(s_max, s[0]);
00545 s[1] = 1 - s[0];
00546
00547 fluid_.mobility(cell, s, c, cmax, mob, dmobds, dmobwatdc);
00548 fluid_.computeMc(c, mc, dmcdc);
00549 fluid_.pc(cell, s, pc, dpc);
00550
00551 store_.mob (cell)[0] = mob [0];
00552 store_.mob (cell)[1] = mob [1];
00553 store_.dmobds(cell)[0] = dmobds[0*2 + 0];
00554 store_.dmobds(cell)[1] = -dmobds[1*2 + 1];
00555 store_.dmobwatdc(cell) = dmobwatdc;
00556 store_.mc(cell) = mc;
00557 store_.dmcdc(cell) = dmcdc;
00558 store_.pc(cell) = pc;
00559 store_.dpc(cell) = dpc;
00560 }
00561 if (!in_range) {
00562 std::cout << "Warning: initIteration() - s was clamped in some cells.\n";
00563 }
00564 return in_range;
00565 }
00566
00567 template <class ReservoirState,
00568 class Grid ,
00569 class NewtonIterate >
00570 void
00571 finishIteration(const ReservoirState& state,
00572 const Grid& g ,
00573 NewtonIterate& it ) {
00574
00575 (void) state; (void) g; (void) it;
00576 typedef typename NewtonIterate::vector_type vector_t;
00577 }
00578
00579 template <class Grid ,
00580 class SolutionVector,
00581 class ReservoirState>
00582 void
00583 finishStep(const Grid& g ,
00584 const SolutionVector& x ,
00585 ReservoirState& state) {
00586
00587 double *s = &state.saturation()[0*2 + 0];
00588 double *c = &state.concentration()[0*1 + 0];
00589 double *cmax = &state.maxconcentration()[0*1 + 0];
00590
00591 for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2, c += 1, cmax +=1) {
00592 s[0] += x[2*cell + 0];
00593 c[0] += x[2*cell + 1];
00594 cmax[0] = std::max(c[0], cmax[0]);
00595 double s_min = fluid_.s_min(cell);
00596 double s_max = fluid_.s_max(cell);
00597 assert(s[0] >= s_min - sat_tol_);
00598 assert(s[0] <= s_max + sat_tol_);
00599 s[0] = std::max(s_min, s[0]);
00600 s[0] = std::min(s_max, s[0]);
00601 s[1] = 1.0 - s[0];
00602 }
00603 }
00604
00605 private:
00606 void
00607 upwindMobility(const double dflux,
00608 const double gflux,
00609 const int* n ,
00610 int* pix ,
00611 double* m ,
00612 double* dmds ,
00613 double& dmobwatdc ,
00614 double& mc,
00615 double& dmcdc) const {
00616 bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
00617 ( (! (dflux > 0)) && (! (gflux > 0)) );
00618
00619 if (equal_sign) {
00620
00621 if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
00622 else { pix[0] = 1; }
00623
00624 m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
00625 mc = store_.mc(n[ pix[0] ]);
00626
00627 if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
00628 else { pix[1] = 1; }
00629
00630 m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
00631
00632 } else {
00633
00634 if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
00635 else { pix[1] = 1; }
00636
00637 m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
00638
00639 if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
00640 else { pix[0] = 1; }
00641
00642 m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
00643 mc = store_.mc(n[ pix[0] ]);
00644 }
00645
00646 dmds[0] = store_.dmobds(n[ pix[0] ]) [ 0 ];
00647 dmds[1] = store_.dmobds(n[ pix[1] ]) [ 1 ];
00648 dmobwatdc = store_.dmobwatdc(n[ pix[0] ]);
00649 dmcdc = store_.dmcdc(n[ pix[0] ]);
00650 }
00651
00652 template <class Grid>
00653 void
00654 computeStaticGravity(const Grid& g) {
00655
00656 const int d = g.dimensions;
00657
00658 for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
00659 const double* cc = g.cell_centroids + (c * d);
00660
00661 for (; i < g.cell_facepos[c + 1]; ++i) {
00662 const int f = g.cell_faces[i];
00663 const double* fc = g.face_centroids + (f * d);
00664
00665 double dg = 0.0;
00666 for (int j = 0; j < d; ++j) {
00667 dg += gravity_[j] * (fc[j] - cc[j]);
00668 }
00669
00670 store_.dg(i) = store_.trans(f) * dg;
00671 }
00672 }
00673 }
00674
00675 double
00676 gravityFlux(const int f) const {
00677 double gflux;
00678
00679 if (gravity_) {
00680 int i1 = f2hf_[2*f + 0];
00681 int i2 = f2hf_[2*f + 1];
00682
00683 assert ((i1 >= 0) && (i2 >= 0));
00684
00685 gflux = store_.dg(i1) - store_.dg(i2);
00686 gflux *= store_.drho();
00687 } else {
00688 gflux = 0.0;
00689 }
00690
00691 return gflux;
00692 }
00693 void
00694 capFlux(const int f,const int* n,double& pcflux, double* dpcflux) const {
00695
00696 int i1 = n[0];
00697 int i2 = n[1];
00698 assert ((i1 >= 0) && (i2 >= 0));
00699
00700 pcflux = store_.trans(f)*(store_.pc(i2) - store_.pc(i1));
00701 dpcflux[0] = -store_.trans(f)*store_.dpc(i1);
00702 dpcflux[1] = store_.trans(f)*store_.dpc(i2);
00703 }
00704
00705 TwophaseFluidPolymer fluid_ ;
00706 const double* gravity_;
00707 std::vector<int> f2hf_ ;
00708 polymer_reorder::ModelParameterStorage store_ ;
00709 bool init_step_use_previous_sol_;
00710 double sat_tol_;
00711 };
00712 }
00713 #endif