OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 02:00:59 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #34 from totto82/bugfixes

Browse Source 
Changes in gravity potential to match Eclipse and Mrst
This commit is contained in:
Bård Skaflestad 2013-09-24 04:02:55 -07:00
commit cf9d417cad
3 changed files with 52 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
opm/autodiff/FullyImplicitBlackoilSolver.cpp
View File

@ -1006,10 +1006,16 @@ namespace {
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 ADB dp = ops_.ngrad * state.pressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
head = transi*dp;
//head = transi*(ops_.ngrad * state.pressure) + gflux;
UpwindSelector<double> upwind(grid_, ops_, head.value());

17
opm/autodiff/GeoProps.hpp
View File

@ -46,6 +46,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;
@ -59,6 +60,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) {
@ -82,15 +89,17 @@ namespace Opm
}
}
const Vector& poreVolume() const { return pvol_ ; }
const Vector& transmissibility() const { return trans_; }
const Vector& gravityPotential() const { return gpot_ ; }
const double* gravity() const { return gravity_; }
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.
};
}

31
opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
View File

@ -199,7 +199,35 @@ namespace Opm
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
static void outputWellStateMatlab(const Opm::WellState& well_state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["bhp"] = &well_state.bhp();
dm["wellrates"] = &well_state.wellRates();
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath);
}
fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
std::ofstream file(fname.str().c_str());
if (!file) {
OPM_THROW(std::runtime_error,"Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
#if 0
static void outputWaterCut(const Opm::Watercut& watercut,
@ -341,6 +369,8 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
}
SimulatorReport sreport;
@ -455,6 +485,7 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
#if 0
outputWaterCut(watercut, output_dir_);
if (wells_) {