Make default initial saturation equal to connate water. Use new total mobilities.

This commit is contained in:
Atgeirr Flø Rasmussen 2012-03-06 16:08:35 +01:00
parent f40aab431b
commit 5b21e24e79

View File

@ -40,6 +40,7 @@
#include <opm/polymer/TransportModelPolymer.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/polymerUtilities.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
@ -103,6 +104,19 @@ public:
}
}
virtual void satRange(const int n,
const int* /*cells*/,
double* smin,
double* smax) const
{
const int np = 2;
for (int i = 0; i < n; ++i) {
smin[np*i + 0] = sw_[0];
smax[np*i + 0] = sw_.back();
smin[np*i + 1] = 1.0 - sw_[0];
smax[np*i + 1] = 1.0 - sw_.back();
}
}
private:
double krw(double s) const
@ -136,17 +150,33 @@ private:
class ReservoirState
{
public:
ReservoirState(const UnstructuredGrid* g, const int num_phases = 2, const double init_sat = 0.0)
ReservoirState(const UnstructuredGrid* g, const double init_sat = 0.0)
: press_ (g->number_of_cells, 0.0),
fpress_(g->number_of_faces, 0.0),
flux_ (g->number_of_faces, 0.0),
sat_ (num_phases * g->number_of_cells, 0.0),
sat_ (2 * g->number_of_cells, 0.0),
concentration_(g->number_of_cells, 0.0),
cmax_(g->number_of_cells, 0.0)
{
for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[num_phases*cell] = init_sat;
sat_[num_phases*cell + num_phases - 1] = 1.0 - init_sat;
sat_[2*cell] = init_sat;
sat_[2*cell + 1] = 1.0 - init_sat;
}
}
void setToMinimumWaterSat(const Opm::IncompPropertiesInterface& props)
{
const int n = props.numCells();
std::vector<int> cells(n);
for (int i = 0; i < n; ++i) {
cells[i] = i;
}
std::vector<double> smin(2*n);
std::vector<double> smax(2*n);
props.satRange(n, &cells[0], &smin[0], &smax[0]);
for (int cell = 0; cell < n; ++cell) {
sat_[2*cell] = smin[2*cell];
sat_[2*cell + 1] = 1.0 - smin[2*cell];
}
}
@ -369,7 +399,10 @@ main(int argc, char** argv)
std::vector<double> totmob;
std::vector<double> omega; // Empty dummy unless/until we include gravity here.
double init_sat = param.getDefault("init_sat", 0.0);
ReservoirState state(grid->c_grid(), props->numPhases(), init_sat);
ReservoirState state(grid->c_grid(), init_sat);
if (!param.has("init_sat")) {
state.setToMinimumWaterSat(*props);
}
// We need a separate reorder_sat, because the reorder
// code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(num_cells);
@ -426,9 +459,11 @@ main(int argc, char** argv)
}
if (use_gravity) {
computeTotalMobilityOmega(*props, allcells, state.saturation(), totmob, omega);
computeTotalMobilityOmega(*props, polydata, allcells, state.saturation(), state.concentration(),
totmob, omega);
} else {
computeTotalMobility(*props, allcells, state.saturation(), totmob);
computeTotalMobility(*props, polydata, allcells, state.saturation(), state.concentration(),
totmob);
}
pressure_timer.start();
psolver.solve(totmob, omega, src, state.pressure(), state.faceflux());