Compute gravity potensial using the face average as in Eclipse and Mrst

This commit is contained in:
Tor Harald Sandve
2013-09-19 14:47:12 +02:00
parent 9109dad0de
commit f3d51c4a22
2 changed files with 20 additions and 2 deletions

View File

@@ -1003,10 +1003,19 @@ namespace Opm {
rq_[ actph ].mob = tr_mult * kr[ phase ] / mu;
const ADB rho = fluidDensity(phase, state.pressure, state.rs, cells_);
const ADB gflux = grav_ * rho;
ADB& head = rq_[ actph ].head;
head = transi*(ops_.ngrad * state.pressure) + gflux;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rho;
const V vz = geo_.z();
std::vector<int> blocksizes = rhoavg.blockPattern();
ADB z = ADB::constant(vz,blocksizes);
const ADB dp = ops_.ngrad * state.pressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * z));
head = transi*dp;
//head = transi*(ops_.ngrad * state.pressure) + gflux;
UpwindSelector<double> upwind(grid_, ops_, head.value());

View File

@@ -41,6 +41,7 @@ namespace Opm
: pvol_ (grid.number_of_cells)
, trans_(grid.number_of_faces)
, gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
, z_(grid.number_of_cells)
{
// Pore volume
const typename Vector::Index nc = grid.number_of_cells;
@@ -54,6 +55,12 @@ namespace Opm
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
tpfa_trans_compute (ug, htrans.data() , trans_.data());
// Compute z coordinates
for (int c = 0; c<nc; ++c){
z_[c] = grid.cell_centroids[c*3 + 2];
}
// Gravity potential
std::fill(gravity_, gravity_ + 3, 0.0);
if (grav != 0) {
@@ -80,12 +87,14 @@ namespace Opm
const Vector& poreVolume() const { return pvol_ ; }
const Vector& transmissibility() const { return trans_; }
const Vector& gravityPotential() const { return gpot_ ; }
const Vector& z() const {return z_;}
const double* gravity() const { return gravity_; }
private:
Vector pvol_ ;
Vector trans_;
Vector gpot_ ;
Vector z_;
double gravity_[3]; // Size 3 even if grid is 2-dim.
};
}