OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
Bård Skaflestad 2012-09-05 11:46:54 +02:00
commit c5cae3268d
46 changed files with 2134 additions and 1053 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Makefile.am
View File

@ -57,6 +57,7 @@ opm/core/fluid/SatFuncStone2.cpp \
opm/core/fluid/SatFuncSimple.cpp \ opm/core/fluid/SatFuncSimple.cpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \
opm/core/fluid/blackoil/SinglePvtDead.cpp \ opm/core/fluid/blackoil/SinglePvtDead.cpp \
opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp \
opm/core/fluid/blackoil/SinglePvtInterface.cpp \ opm/core/fluid/blackoil/SinglePvtInterface.cpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \ opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \ opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \
@ -150,15 +151,18 @@ opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \
opm/core/fluid/RockBasic.hpp \ opm/core/fluid/RockBasic.hpp \
opm/core/fluid/RockCompressibility.hpp \ opm/core/fluid/RockCompressibility.hpp \
opm/core/fluid/RockFromDeck.hpp \ opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \
opm/core/fluid/SatFuncStone2.hpp \ opm/core/fluid/SatFuncStone2.hpp \
opm/core/fluid/SatFuncSimple.hpp \ opm/core/fluid/SatFuncSimple.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \
opm/core/fluid/SaturationPropsFromDeck_impl.hpp \
opm/core/fluid/SaturationPropsInterface.hpp \
opm/core/fluid/SimpleFluid2p.hpp \ opm/core/fluid/SimpleFluid2p.hpp \
opm/core/fluid/blackoil/BlackoilPhases.hpp \ opm/core/fluid/blackoil/BlackoilPhases.hpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \ opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \
opm/core/fluid/blackoil/SinglePvtDead.hpp \ opm/core/fluid/blackoil/SinglePvtDead.hpp \
opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp \
opm/core/fluid/blackoil/SinglePvtInterface.hpp \ opm/core/fluid/blackoil/SinglePvtInterface.hpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \ opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \ opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \
@ -233,6 +237,7 @@ opm/core/utility/ColumnExtract.hpp \
opm/core/utility/ErrorMacros.hpp \ opm/core/utility/ErrorMacros.hpp \
opm/core/utility/Factory.hpp \ opm/core/utility/Factory.hpp \
opm/core/utility/MonotCubicInterpolator.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \
opm/core/utility/NonuniformTableLinear.hpp \
opm/core/utility/RootFinders.hpp \ opm/core/utility/RootFinders.hpp \
opm/core/utility/SparseTable.hpp \ opm/core/utility/SparseTable.hpp \
opm/core/utility/SparseVector.hpp \ opm/core/utility/SparseVector.hpp \

2
examples/sim_2p_comp_reorder.cpp
View File

@ -94,7 +94,7 @@ main(int argc, char** argv)
// Grid init // Grid init
grid.reset(new GridManager(*deck)); grid.reset(new GridManager(*deck));
// Rock and fluid init // Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid())); props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param));
// check_well_controls = param.getDefault("check_well_controls", false); // check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility. // Rock compressibility.

2
examples/sim_wateroil.cpp
View File

@ -175,7 +175,7 @@ main(int argc, char** argv)
// Grid init // Grid init
grid.reset(new Opm::GridManager(deck)); grid.reset(new Opm::GridManager(deck));
// Rock and fluid init // Rock and fluid init
props.reset(new Opm::BlackoilPropertiesFromDeck(deck, *grid->c_grid())); props.reset(new BlackoilPropertiesFromDeck(deck, *grid->c_grid(), param));
// Wells init. // Wells init.
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
check_well_controls = param.getDefault("check_well_controls", false); check_well_controls = param.getDefault("check_well_controls", false);

72
opm/core/fluid/BlackoilPropertiesBasic.cpp
View File

@ -26,20 +26,20 @@ namespace Opm
{ {
BlackoilPropertiesBasic::BlackoilPropertiesBasic(const parameter::ParameterGroup& param, BlackoilPropertiesBasic::BlackoilPropertiesBasic(const parameter::ParameterGroup& param,
const int dim, const int dim,
const int num_cells) const int num_cells)
{ {
double poro = param.getDefault("porosity", 1.0); double poro = param.getDefault("porosity", 1.0);
using namespace Opm::unit; using namespace Opm::unit;
using namespace Opm::prefix; using namespace Opm::prefix;
double perm = param.getDefault("permeability", 100.0)*milli*darcy; double perm = param.getDefault("permeability", 100.0)*milli*darcy;
rock_.init(dim, num_cells, poro, perm); rock_.init(dim, num_cells, poro, perm);
pvt_.init(param); pvt_.init(param);
satprops_.init(param); satprops_.init(param);
if (pvt_.numPhases() != satprops_.numPhases()) { if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
} }
} }
BlackoilPropertiesBasic::~BlackoilPropertiesBasic() BlackoilPropertiesBasic::~BlackoilPropertiesBasic()
@ -90,11 +90,11 @@ namespace Opm
/// \param[out] dmudp If non-null: array of nP viscosity derivative values, /// \param[out] dmudp If non-null: array of nP viscosity derivative values,
/// array must be valid before calling. /// array must be valid before calling.
void BlackoilPropertiesBasic::viscosity(const int n, void BlackoilPropertiesBasic::viscosity(const int n,
const double* p, const double* p,
const double* z, const double* z,
const int* /*cells*/, const int* /*cells*/,
double* mu, double* mu,
double* dmudp) const double* dmudp) const
{ {
if (dmudp) { if (dmudp) {
THROW("BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented."); THROW("BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented.");
@ -114,16 +114,16 @@ namespace Opm
/// array must be valid before calling. The matrices are output /// array must be valid before calling. The matrices are output
/// in Fortran order. /// in Fortran order.
void BlackoilPropertiesBasic::matrix(const int n, void BlackoilPropertiesBasic::matrix(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*z*/, const double* /*z*/,
const int* /*cells*/, const int* /*cells*/,
double* A, double* A,
double* dAdp) const double* dAdp) const
{ {
const int np = numPhases(); const int np = numPhases();
ASSERT(np <= 2); ASSERT(np <= 2);
double B[2]; // Must be enough since component classes do not handle more than 2. double B[2]; // Must be enough since component classes do not handle more than 2.
pvt_.B(1, 0, 0, B); pvt_.B(1, 0, 0, B);
// Compute A matrix // Compute A matrix
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
@ -152,8 +152,8 @@ namespace Opm
/// of a call to the method matrix(). /// of a call to the method matrix().
/// \param[out] rho Array of nP density values, array must be valid before calling. /// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesBasic::density(const int n, void BlackoilPropertiesBasic::density(const int n,
const double* A, const double* A,
double* rho) const double* rho) const
{ {
const int np = numPhases(); const int np = numPhases();
const double* sdens = pvt_.surfaceDensities(); const double* sdens = pvt_.surfaceDensities();
@ -186,10 +186,10 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j}, /// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void BlackoilPropertiesBasic::relperm(const int n, void BlackoilPropertiesBasic::relperm(const int n,
const double* s, const double* s,
const int* /*cells*/, const int* /*cells*/,
double* kr, double* kr,
double* dkrds) const double* dkrds) const
{ {
satprops_.relperm(n, s, kr, dkrds); satprops_.relperm(n, s, kr, dkrds);
} }
@ -205,10 +205,10 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j}, /// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void BlackoilPropertiesBasic::capPress(const int n, void BlackoilPropertiesBasic::capPress(const int n,
const double* s, const double* s,
const int* /*cells*/, const int* /*cells*/,
double* pc, double* pc,
double* dpcds) const double* dpcds) const
{ {
satprops_.capPress(n, s, pc, dpcds); satprops_.capPress(n, s, pc, dpcds);
} }
@ -226,7 +226,7 @@ namespace Opm
double* smin, double* smin,
double* smax) const double* smax) const
{ {
satprops_.satRange(n, smin, smax); satprops_.satRange(n, smin, smax);
} }

20
opm/core/fluid/BlackoilPropertiesBasic.hpp
View File

@ -35,16 +35,16 @@ namespace Opm
{ {
public: public:
/// Construct from parameters. /// Construct from parameters.
/// The following parameters are accepted (defaults): /// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2. /// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP /// mu1 [mu2, mu3] (1.0) Viscosity in cP
/// porosity (1.0) Porosity /// porosity (1.0) Porosity
/// permeability (100.0) Permeability in mD /// permeability (100.0) Permeability in mD
BlackoilPropertiesBasic(const parameter::ParameterGroup& param, BlackoilPropertiesBasic(const parameter::ParameterGroup& param,
const int dim, const int dim,
const int num_cells); const int num_cells);
/// Destructor. /// Destructor.
virtual ~BlackoilPropertiesBasic(); virtual ~BlackoilPropertiesBasic();
@ -151,7 +151,7 @@ namespace Opm
double* dpcds) const; double* dpcds) const;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be /// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]]. /// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points. /// \param[in] n Number of data points.

64
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -19,6 +19,7 @@
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp> #include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm namespace Opm
{ {
@ -26,12 +27,16 @@ namespace Opm
const UnstructuredGrid& grid) const UnstructuredGrid& grid)
{ {
rock_.init(deck, grid); rock_.init(deck, grid);
pvt_.init(deck); pvt_.init(deck, 200);
satprops_.init(deck, grid); SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
if (pvt_.numPhases() != satprops_.numPhases()) { = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" satprops_.reset(ptr);
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); ptr->init(deck, grid, 200);
}
if (pvt_.numPhases() != satprops_->numPhases()) {
THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
} }
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
@ -39,13 +44,42 @@ namespace Opm
const parameter::ParameterGroup& param) const parameter::ParameterGroup& param)
{ {
rock_.init(deck, grid); rock_.init(deck, grid);
const int samples = param.getDefault("dead_tab_size", 1025); const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(deck, samples); pvt_.init(deck, pvt_samples);
satprops_.init(deck, grid, param);
if (pvt_.numPhases() != satprops_.numPhases()) { // Unfortunate lack of pointer smartness here...
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" const int sat_samples = param.getDefault("sat_tab_size", 200);
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
bool use_stone2 = (threephase_model == "stone2");
if (sat_samples > 1) {
if (use_stone2) {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
}
} else {
if (use_stone2) {
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
}
}
if (pvt_.numPhases() != satprops_->numPhases()) {
THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
} }
} }
@ -250,7 +284,7 @@ namespace Opm
double* kr, double* kr,
double* dkrds) const double* dkrds) const
{ {
satprops_.relperm(n, s, cells, kr, dkrds); satprops_->relperm(n, s, cells, kr, dkrds);
} }
@ -269,7 +303,7 @@ namespace Opm
double* pc, double* pc,
double* dpcds) const double* dpcds) const
{ {
satprops_.capPress(n, s, cells, pc, dpcds); satprops_->capPress(n, s, cells, pc, dpcds);
} }
@ -285,7 +319,7 @@ namespace Opm
double* smin, double* smin,
double* smax) const double* smax) const
{ {
satprops_.satRange(n, cells, smin, smax); satprops_->satRange(n, cells, smin, smax);
} }

16
opm/core/fluid/BlackoilPropertiesFromDeck.hpp
View File

@ -27,6 +27,7 @@
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <boost/scoped_ptr.hpp>
struct UnstructuredGrid; struct UnstructuredGrid;
@ -52,8 +53,11 @@ namespace Opm
/// mapping from cell indices (typically from a processed grid) /// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck. /// to logical cartesian indices consistent with the deck.
/// \param[in] param Parameters. Accepted parameters include: /// \param[in] param Parameters. Accepted parameters include:
/// dead_tab_size (1025) number of uniform sample points for dead-oil pvt tables. /// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables.
/// tab_size_kr (200) number of uniform sample points for saturation tables. /// sat_tab_size (200) number of uniform sample points for saturation tables.
/// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2").
/// For both size parameters, a 0 or negative value indicates that no spline fitting is to
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck, BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid, const UnstructuredGrid& grid,
const parameter::ParameterGroup& param); const parameter::ParameterGroup& param);
@ -163,9 +167,9 @@ namespace Opm
double* dpcds) const; double* dpcds) const;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be /// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]]. /// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices. /// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
@ -178,7 +182,7 @@ namespace Opm
private: private:
RockFromDeck rock_; RockFromDeck rock_;
BlackoilPvtProperties pvt_; BlackoilPvtProperties pvt_;
SaturationPropsFromDeck satprops_; boost::scoped_ptr<SaturationPropsInterface> satprops_;
mutable std::vector<double> B_; mutable std::vector<double> B_;
mutable std::vector<double> dB_; mutable std::vector<double> dB_;
mutable std::vector<double> R_; mutable std::vector<double> R_;

6
opm/core/fluid/BlackoilPropertiesInterface.hpp
View File

@ -138,9 +138,9 @@ namespace Opm
double* dpcds) const = 0; double* dpcds) const = 0;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be /// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]]. /// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices. /// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.

70
opm/core/fluid/IncompPropertiesBasic.cpp
View File

@ -28,22 +28,22 @@ namespace Opm
{ {
IncompPropertiesBasic::IncompPropertiesBasic(const parameter::ParameterGroup& param, IncompPropertiesBasic::IncompPropertiesBasic(const parameter::ParameterGroup& param,
const int dim, const int dim,
const int num_cells) const int num_cells)
{ {
double poro = param.getDefault("porosity", 1.0); double poro = param.getDefault("porosity", 1.0);
using namespace Opm::unit; using namespace Opm::unit;
using namespace Opm::prefix; using namespace Opm::prefix;
double perm = param.getDefault("permeability", 100.0)*milli*darcy; double perm = param.getDefault("permeability", 100.0)*milli*darcy;
rock_.init(dim, num_cells, poro, perm); rock_.init(dim, num_cells, poro, perm);
pvt_.init(param); pvt_.init(param);
satprops_.init(param); satprops_.init(param);
if (pvt_.numPhases() != satprops_.numPhases()) { if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
} }
viscosity_.resize(pvt_.numPhases()); viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]); pvt_.mu(1, 0, 0, &viscosity_[0]);
} }
IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases, IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases,
@ -56,14 +56,14 @@ namespace Opm
const int num_cells) const int num_cells)
{ {
rock_.init(dim, num_cells, por, perm); rock_.init(dim, num_cells, por, perm);
pvt_.init(num_phases, rho, mu); pvt_.init(num_phases, rho, mu);
satprops_.init(num_phases, relpermfunc); satprops_.init(num_phases, relpermfunc);
if (pvt_.numPhases() != satprops_.numPhases()) { if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
} }
viscosity_.resize(pvt_.numPhases()); viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]); pvt_.mu(1, 0, 0, &viscosity_[0]);
} }
IncompPropertiesBasic::~IncompPropertiesBasic() IncompPropertiesBasic::~IncompPropertiesBasic()
@ -109,7 +109,7 @@ namespace Opm
/// \return Array of P viscosity values. /// \return Array of P viscosity values.
const double* IncompPropertiesBasic::viscosity() const const double* IncompPropertiesBasic::viscosity() const
{ {
return &viscosity_[0]; return &viscosity_[0];
} }
/// \return Array of P density values. /// \return Array of P density values.
@ -117,7 +117,7 @@ namespace Opm
{ {
// No difference between reservoir and surface densities // No difference between reservoir and surface densities
// modelled by this class. // modelled by this class.
return pvt_.surfaceDensities(); return pvt_.surfaceDensities();
} }
/// \return Array of P density values. /// \return Array of P density values.
@ -125,7 +125,7 @@ namespace Opm
{ {
// No difference between reservoir and surface densities // No difference between reservoir and surface densities
// modelled by this class. // modelled by this class.
return pvt_.surfaceDensities(); return pvt_.surfaceDensities();
} }
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
@ -138,10 +138,10 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j}, /// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesBasic::relperm(const int n, void IncompPropertiesBasic::relperm(const int n,
const double* s, const double* s,
const int* /*cells*/, const int* /*cells*/,
double* kr, double* kr,
double* dkrds) const double* dkrds) const
{ {
satprops_.relperm(n, s, kr, dkrds); satprops_.relperm(n, s, kr, dkrds);
} }
@ -157,10 +157,10 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j}, /// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesBasic::capPress(const int n, void IncompPropertiesBasic::capPress(const int n,
const double* s, const double* s,
const int* /*cells*/, const int* /*cells*/,
double* pc, double* pc,
double* dpcds) const double* dpcds) const
{ {
satprops_.capPress(n, s, pc, dpcds); satprops_.capPress(n, s, pc, dpcds);
} }
@ -174,11 +174,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void IncompPropertiesBasic::satRange(const int n, void IncompPropertiesBasic::satRange(const int n,
const int* /*cells*/, const int* /*cells*/,
double* smin, double* smin,
double* smax) const double* smax) const
{ {
satprops_.satRange(n, smin, smax); satprops_.satRange(n, smin, smax);
} }
} // namespace Opm } // namespace Opm

34
opm/core/fluid/IncompPropertiesBasic.hpp
View File

@ -42,29 +42,29 @@ namespace Opm
{ {
public: public:
/// Construct from parameters. /// Construct from parameters.
/// The following parameters are accepted (defaults): /// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2. /// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP /// mu1 [mu2, mu3] (1.0) Viscosity in cP
/// porosity (1.0) Porosity /// porosity (1.0) Porosity
/// permeability (100.0) Permeability in mD /// permeability (100.0) Permeability in mD
IncompPropertiesBasic(const parameter::ParameterGroup& param, IncompPropertiesBasic(const parameter::ParameterGroup& param,
const int dim, const int dim,
const int num_cells); const int num_cells);
/// Construct from arguments a basic two phase fluid. /// Construct from arguments a basic two phase fluid.
IncompPropertiesBasic(const int num_phases, IncompPropertiesBasic(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc, const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho, const std::vector<double>& rho,
const std::vector<double>& mu, const std::vector<double>& mu,
const double porosity, const double porosity,
const double permeability, const double permeability,
const int dim, const int dim,
const int num_cells); const int num_cells);
/// Destructor. /// Destructor.
virtual ~IncompPropertiesBasic(); virtual ~IncompPropertiesBasic();
// ---- Rock interface ---- // ---- Rock interface ----
@ -132,9 +132,9 @@ namespace Opm
double* dpcds) const; double* dpcds) const;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be /// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]]. /// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices. /// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
@ -145,9 +145,9 @@ namespace Opm
double* smax) const; double* smax) const;
private: private:
RockBasic rock_; RockBasic rock_;
PvtPropertiesBasic pvt_; PvtPropertiesBasic pvt_;
SaturationPropsBasic satprops_; SaturationPropsBasic satprops_;
std::vector<double> viscosity_; std::vector<double> viscosity_;
}; };

44
opm/core/fluid/IncompPropertiesFromDeck.cpp
View File

@ -27,15 +27,15 @@ namespace Opm
{ {
IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck, IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid) const UnstructuredGrid& grid)
{ {
rock_.init(deck, grid); rock_.init(deck, grid);
pvt_.init(deck); pvt_.init(deck);
satprops_.init(deck, grid); satprops_.init(deck, grid, 200);
if (pvt_.numPhases() != satprops_.numPhases()) { if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
} }
} }
IncompPropertiesFromDeck::~IncompPropertiesFromDeck() IncompPropertiesFromDeck::~IncompPropertiesFromDeck()
@ -81,19 +81,19 @@ namespace Opm
/// \return Array of P viscosity values. /// \return Array of P viscosity values.
const double* IncompPropertiesFromDeck::viscosity() const const double* IncompPropertiesFromDeck::viscosity() const
{ {
return pvt_.viscosity(); return pvt_.viscosity();
} }
/// \return Array of P density values. /// \return Array of P density values.
const double* IncompPropertiesFromDeck::density() const const double* IncompPropertiesFromDeck::density() const
{ {
return pvt_.reservoirDensities(); return pvt_.reservoirDensities();
} }
/// \return Array of P density values. /// \return Array of P density values.
const double* IncompPropertiesFromDeck::surfaceDensity() const const double* IncompPropertiesFromDeck::surfaceDensity() const
{ {
return pvt_.surfaceDensities(); return pvt_.surfaceDensities();
} }
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
@ -106,10 +106,10 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j}, /// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesFromDeck::relperm(const int n, void IncompPropertiesFromDeck::relperm(const int n,
const double* s, const double* s,
const int* cells, const int* cells,
double* kr, double* kr,
double* dkrds) const double* dkrds) const
{ {
satprops_.relperm(n, s, cells, kr, dkrds); satprops_.relperm(n, s, cells, kr, dkrds);
} }
@ -125,10 +125,10 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j}, /// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesFromDeck::capPress(const int n, void IncompPropertiesFromDeck::capPress(const int n,
const double* s, const double* s,
const int* cells, const int* cells,
double* pc, double* pc,
double* dpcds) const double* dpcds) const
{ {
satprops_.capPress(n, s, cells, pc, dpcds); satprops_.capPress(n, s, cells, pc, dpcds);
} }
@ -142,11 +142,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void IncompPropertiesFromDeck::satRange(const int n, void IncompPropertiesFromDeck::satRange(const int n,
const int* cells, const int* cells,
double* smin, double* smin,
double* smax) const double* smax) const
{ {
satprops_.satRange(n, cells, smin, smax); satprops_.satRange(n, cells, smin, smax);
} }
} // namespace Opm } // namespace Opm

14
opm/core/fluid/IncompPropertiesFromDeck.hpp
View File

@ -51,9 +51,9 @@ namespace Opm
/// mapping from cell indices (typically from a processed grid) /// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck. /// to logical cartesian indices consistent with the deck.
IncompPropertiesFromDeck(const EclipseGridParser& deck, IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid); const UnstructuredGrid& grid);
/// Destructor. /// Destructor.
virtual ~IncompPropertiesFromDeck(); virtual ~IncompPropertiesFromDeck();
// ---- Rock interface ---- // ---- Rock interface ----
@ -121,9 +121,9 @@ namespace Opm
double* dpcds) const; double* dpcds) const;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be /// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]]. /// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices. /// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
@ -134,8 +134,8 @@ namespace Opm
double* smax) const; double* smax) const;
private: private:
RockFromDeck rock_; RockFromDeck rock_;
PvtPropertiesIncompFromDeck pvt_; PvtPropertiesIncompFromDeck pvt_;
SaturationPropsFromDeck satprops_; SaturationPropsFromDeck<SatFuncStone2Uniform> satprops_;
}; };

6
opm/core/fluid/IncompPropertiesInterface.hpp
View File

@ -109,9 +109,9 @@ namespace Opm
double* pc, double* pc,
double* dpcds) const = 0; double* dpcds) const = 0;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// In cell cells[i], saturation of phase p is allowed to be /// In cell cells[i], saturation of phase p is allowed to be
/// in the interval [smin[i*P + p], smax[i*P + p]]. /// in the interval [smin[i*P + p], smax[i*P + p]].
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices. /// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.

112
opm/core/fluid/PvtPropertiesBasic.cpp
View File

@ -34,41 +34,41 @@ namespace Opm
void PvtPropertiesBasic::init(const parameter::ParameterGroup& param) void PvtPropertiesBasic::init(const parameter::ParameterGroup& param)
{ {
int num_phases = param.getDefault("num_phases", 2); int num_phases = param.getDefault("num_phases", 2);
if (num_phases > 3 || num_phases < 1) { if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
} }
density_.resize(num_phases); density_.resize(num_phases);
viscosity_.resize(num_phases); viscosity_.resize(num_phases);
// We currently do not allow the user to set B. // We currently do not allow the user to set B.
formation_volume_factor_.clear(); formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0); formation_volume_factor_.resize(num_phases, 1.0);
// Setting mu and rho from parameters // Setting mu and rho from parameters
using namespace Opm::prefix; using namespace Opm::prefix;
using namespace Opm::unit; using namespace Opm::unit;
const double kgpm3 = kilogram/cubic(meter); const double kgpm3 = kilogram/cubic(meter);
const double cP = centi*Poise; const double cP = centi*Poise;
std::string rname[3] = { "rho1", "rho2", "rho3" }; std::string rname[3] = { "rho1", "rho2", "rho3" };
double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 }; double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 };
std::string vname[3] = { "mu1", "mu2", "mu3" }; std::string vname[3] = { "mu1", "mu2", "mu3" };
double vdefault[3] = { 1.0, 1.0, 1.0 }; double vdefault[3] = { 1.0, 1.0, 1.0 };
for (int phase = 0; phase < num_phases; ++phase) { for (int phase = 0; phase < num_phases; ++phase) {
density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]); density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]);
viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]); viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]);
} }
} }
void PvtPropertiesBasic::init(const int num_phases, void PvtPropertiesBasic::init(const int num_phases,
const std::vector<double>& rho, const std::vector<double>& rho,
const std::vector<double>& visc) const std::vector<double>& visc)
{ {
if (num_phases > 3 || num_phases < 1) { if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
} }
// We currently do not allow the user to set B. // We currently do not allow the user to set B.
formation_volume_factor_.clear(); formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0); formation_volume_factor_.resize(num_phases, 1.0);
density_ = rho; density_ = rho;
viscosity_ = visc; viscosity_ = visc;
} }
@ -87,69 +87,69 @@ namespace Opm
void PvtPropertiesBasic::mu(const int n, void PvtPropertiesBasic::mu(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*z*/, const double* /*z*/,
double* output_mu) const double* output_mu) const
{ {
const int np = numPhases(); const int np = numPhases();
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
output_mu[np*i + phase] = viscosity_[phase]; output_mu[np*i + phase] = viscosity_[phase];
} }
} }
} }
void PvtPropertiesBasic::B(const int n, void PvtPropertiesBasic::B(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*z*/, const double* /*z*/,
double* output_B) const double* output_B) const
{ {
const int np = numPhases(); const int np = numPhases();
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
output_B[np*i + phase] = formation_volume_factor_[phase]; output_B[np*i + phase] = formation_volume_factor_[phase];
} }
} }
} }
void PvtPropertiesBasic::dBdp(const int n, void PvtPropertiesBasic::dBdp(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*z*/, const double* /*z*/,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
{ {
const int np = numPhases(); const int np = numPhases();
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
output_B[np*i + phase] = formation_volume_factor_[phase]; output_B[np*i + phase] = formation_volume_factor_[phase];
output_dBdp[np*i + phase] = 0.0; output_dBdp[np*i + phase] = 0.0;
} }
} }
} }
void PvtPropertiesBasic::R(const int n, void PvtPropertiesBasic::R(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*z*/, const double* /*z*/,
double* output_R) const double* output_R) const
{ {
const int np = numPhases(); const int np = numPhases();
std::fill(output_R, output_R + n*np, 0.0); std::fill(output_R, output_R + n*np, 0.0);
} }
void PvtPropertiesBasic::dRdp(const int n, void PvtPropertiesBasic::dRdp(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*z*/, const double* /*z*/,
double* output_R, double* output_R,
double* output_dRdp) const double* output_dRdp) const
{ {
const int np = numPhases(); const int np = numPhases();
std::fill(output_R, output_R + n*np, 0.0); std::fill(output_R, output_R + n*np, 0.0);
std::fill(output_dRdp, output_dRdp + n*np, 0.0); std::fill(output_dRdp, output_dRdp + n*np, 0.0);
} }
} // namespace Opm } // namespace Opm

18
opm/core/fluid/PvtPropertiesBasic.hpp
View File

@ -38,11 +38,11 @@ namespace Opm
PvtPropertiesBasic(); PvtPropertiesBasic();
/// Initialize from parameters. /// Initialize from parameters.
/// The following parameters are accepted (defaults): /// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1, 2 or 3. /// num_phases (2) Must be 1, 2 or 3.
/// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3
/// mu1 [mu2, mu3] (1.0) Viscosity in cP /// mu1 [mu2, mu3] (1.0) Viscosity in cP
void init(const parameter::ParameterGroup& param); void init(const parameter::ParameterGroup& param);
/// Initialize from arguments. /// Initialize from arguments.
/// Basic multi phase fluid pvt properties. /// Basic multi phase fluid pvt properties.
@ -55,7 +55,7 @@ namespace Opm
/// Densities of stock components at surface conditions. /// Densities of stock components at surface conditions.
/// \return Array of size numPhases(). /// \return Array of size numPhases().
const double* surfaceDensities() const; const double* surfaceDensities() const;
/// Viscosity as a function of p and z. /// Viscosity as a function of p and z.
void mu(const int n, void mu(const int n,
@ -90,9 +90,9 @@ namespace Opm
double* output_dRdp) const; double* output_dRdp) const;
private: private:
std::vector<double> density_; std::vector<double> density_;
std::vector<double> viscosity_; std::vector<double> viscosity_;
std::vector<double> formation_volume_factor_; std::vector<double> formation_volume_factor_;
}; };
} }

70
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
View File

@ -38,54 +38,54 @@ namespace Opm
{ {
typedef std::vector<std::vector<std::vector<double> > > table_t; typedef std::vector<std::vector<std::vector<double> > > table_t;
// If we need multiple regions, this class and the SinglePvt* classes must change. // If we need multiple regions, this class and the SinglePvt* classes must change.
int region_number = 0; int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(deck); PhaseUsage phase_usage = phaseUsageFromDeck(deck);
if (phase_usage.phase_used[PhaseUsage::Vapour] || if (phase_usage.phase_used[PhaseUsage::Vapour] ||
!phase_usage.phase_used[PhaseUsage::Aqua] || !phase_usage.phase_used[PhaseUsage::Aqua] ||
!phase_usage.phase_used[PhaseUsage::Liquid]) { !phase_usage.phase_used[PhaseUsage::Liquid]) {
THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n"); THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n");
} }
// Surface densities. Accounting for different orders in eclipse and our code. // Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) { if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number]; const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water]; surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil]; surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
} else { } else {
THROW("Input is missing DENSITY\n"); THROW("Input is missing DENSITY\n");
} }
// Make reservoir densities the same as surface densities initially. // Make reservoir densities the same as surface densities initially.
// We will modify them with formation volume factors if found. // We will modify them with formation volume factors if found.
reservoir_density_ = surface_density_; reservoir_density_ = surface_density_;
// Water viscosity. // Water viscosity.
if (deck.hasField("PVTW")) { if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number]; const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) { if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
MESSAGE("Compressibility effects in PVTW are ignored."); MESSAGE("Compressibility effects in PVTW are ignored.");
} }
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1]; reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3]; viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else { } else {
// Eclipse 100 default. // Eclipse 100 default.
// viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise; // viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise;
THROW("Input is missing PVTW\n"); THROW("Input is missing PVTW\n");
} }
// Oil viscosity. // Oil viscosity.
if (deck.hasField("PVCDO")) { if (deck.hasField("PVCDO")) {
const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number]; const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number];
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) { if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
MESSAGE("Compressibility effects in PVCDO are ignored."); MESSAGE("Compressibility effects in PVCDO are ignored.");
} }
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1]; reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3]; viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else { } else {
THROW("Input is missing PVCDO\n"); THROW("Input is missing PVCDO\n");
} }
} }
const double* PvtPropertiesIncompFromDeck::surfaceDensities() const const double* PvtPropertiesIncompFromDeck::surfaceDensities() const

12
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
View File

@ -39,14 +39,14 @@ namespace Opm
PvtPropertiesIncompFromDeck(); PvtPropertiesIncompFromDeck();
/// Initialize from deck. /// Initialize from deck.
void init(const EclipseGridParser& deck); void init(const EclipseGridParser& deck);
/// Number of active phases. /// Number of active phases.
int numPhases() const; int numPhases() const;
/// Densities of stock components at surface conditions. /// Densities of stock components at surface conditions.
/// \return Array of size numPhases(). /// \return Array of size numPhases().
const double* surfaceDensities() const; const double* surfaceDensities() const;
/// Densities of stock components at reservoir conditions. /// Densities of stock components at reservoir conditions.
/// Note: a reasonable question to ask is why there can be /// Note: a reasonable question to ask is why there can be
@ -58,15 +58,15 @@ namespace Opm
/// reporting and using data given in terms of surface values, /// reporting and using data given in terms of surface values,
/// we need to handle this difference. /// we need to handle this difference.
/// \return Array of size numPhases(). /// \return Array of size numPhases().
const double* reservoirDensities() const; const double* reservoirDensities() const;
/// Viscosities. /// Viscosities.
const double* viscosity() const; const double* viscosity() const;
private: private:
std::tr1::array<double, 2> surface_density_; std::tr1::array<double, 2> surface_density_;
std::tr1::array<double, 2> reservoir_density_; std::tr1::array<double, 2> reservoir_density_;
std::tr1::array<double, 2> viscosity_; std::tr1::array<double, 2> viscosity_;
}; };
} }

30
opm/core/fluid/RockBasic.cpp
View File

@ -25,29 +25,29 @@ namespace Opm
/// Default constructor. /// Default constructor.
RockBasic::RockBasic() RockBasic::RockBasic()
: dimensions_(-1) : dimensions_(-1)
{ {
} }
/// Initialize with homogenous porosity and permeability. /// Initialize with homogenous porosity and permeability.
void RockBasic::init(const int dimensions, void RockBasic::init(const int dimensions,
const int num_cells, const int num_cells,
const double poro, const double poro,
const double perm) const double perm)
{ {
dimensions_ = dimensions; dimensions_ = dimensions;
porosity_.clear(); porosity_.clear();
porosity_.resize(num_cells, poro); porosity_.resize(num_cells, poro);
permeability_.clear(); permeability_.clear();
const int dsq = dimensions*dimensions; const int dsq = dimensions*dimensions;
permeability_.resize(num_cells*dsq, 0.0); permeability_.resize(num_cells*dsq, 0.0);
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < num_cells; ++i) { for (int i = 0; i < num_cells; ++i) {
for (int d = 0; d < dimensions; ++d) { for (int d = 0; d < dimensions; ++d) {
permeability_[dsq*i + dimensions*d + d] = perm; permeability_[dsq*i + dimensions*d + d] = perm;
} }
} }
} }

8
opm/core/fluid/RockBasic.hpp
View File

@ -35,9 +35,9 @@ namespace Opm
/// Initialize with homogenous porosity and permeability. /// Initialize with homogenous porosity and permeability.
void init(const int dimensions, void init(const int dimensions,
const int num_cells, const int num_cells,
const double poro, const double poro,
const double perm); const double perm);
/// \return D, the number of spatial dimensions. /// \return D, the number of spatial dimensions.
int numDimensions() const int numDimensions() const
@ -66,7 +66,7 @@ namespace Opm
} }
private: private:
int dimensions_; int dimensions_;
std::vector<double> porosity_; std::vector<double> porosity_;
std::vector<double> permeability_; std::vector<double> permeability_;
}; };

6
opm/core/fluid/RockCompressibility.cpp
View File

@ -69,8 +69,8 @@ namespace Opm
const double cpnorm = rock_comp_*(pressure - pref_); const double cpnorm = rock_comp_*(pressure - pref_);
return (1.0 + cpnorm + 0.5*cpnorm*cpnorm); return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
} else { } else {
// return Opm::linearInterpolation(p_, poromult_, pressure); // return Opm::linearInterpolation(p_, poromult_, pressure);
return Opm::linearInterpolationExtrap(p_, poromult_, pressure); return Opm::linearInterpolationExtrap(p_, poromult_, pressure);
} }
} }
@ -81,7 +81,7 @@ namespace Opm
} else { } else {
//const double poromult = Opm::linearInterpolation(p_, poromult_, pressure); //const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
//const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure); //const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure); const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivativeExtrap(p_, poromult_, pressure); const double dporomultdp = Opm::linearInterpolationDerivativeExtrap(p_, poromult_, pressure);
return dporomultdp/poromult; return dporomultdp/poromult;

237
opm/core/fluid/SatFuncSimple.cpp
View File

@ -1,3 +1,22 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SatFuncSimple.hpp> #include <opm/core/fluid/SatFuncSimple.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp> #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
@ -6,21 +25,14 @@
#include <opm/core/utility/buildUniformMonotoneTable.hpp> #include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <iostream> #include <iostream>
namespace Opm namespace Opm
{ {
void SatFuncSimple::init(const EclipseGridParser& deck, void SatFuncSimpleUniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg)
{
init(deck, table_num, phase_usg, 200);
}
void SatFuncSimple::init(const EclipseGridParser& deck,
const int table_num, const int table_num,
const PhaseUsage phase_usg, const PhaseUsage phase_usg,
const int samples) const int samples)
@ -65,7 +77,7 @@ namespace Opm
} }
void SatFuncSimple::evalKr(const double* s, double* kr) const void SatFuncSimpleUniform::evalKr(const double* s, double* kr) const
{ {
if (phase_usage.num_phases == 3) { if (phase_usage.num_phases == 3) {
// A simplified relative permeability model. // A simplified relative permeability model.
@ -106,7 +118,7 @@ namespace Opm
} }
void SatFuncSimple::evalKrDeriv(const double* s, double* kr, double* dkrds) const void SatFuncSimpleUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{ {
const int np = phase_usage.num_phases; const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0); std::fill(dkrds, dkrds + np*np, 0.0);
@ -172,7 +184,7 @@ namespace Opm
} }
void SatFuncSimple::evalPc(const double* s, double* pc) const void SatFuncSimpleUniform::evalPc(const double* s, double* pc) const
{ {
pc[phase_usage.phase_pos[Liquid]] = 0.0; pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) { if (phase_usage.phase_used[Aqua]) {
@ -185,7 +197,206 @@ namespace Opm
} }
} }
void SatFuncSimple::evalPcDeriv(const double* s, double* pc, double* dpcds) const void SatFuncSimpleUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncSimpleNonuniform ======
void SatFuncSimpleNonuniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int /*samples*/)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
krw_ = NonuniformTableLinear<double>(sw, krw);
krow_ = NonuniformTableLinear<double>(sw, krow);
pcow_ = NonuniformTableLinear<double>(sw, pcow);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
krg_ = NonuniformTableLinear<double>(sg, krg);
krog_ = NonuniformTableLinear<double>(sg, krog);
pcog_ = NonuniformTableLinear<double>(sg, pcog);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncSimpleNonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncSimpleNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[Liquid + Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + Vapour*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncSimpleNonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncSimpleNonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{ {
// The problem of determining three-phase capillary pressures // The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase // is very hard experimentally, usually one extends two-phase

36
opm/core/fluid/SatFuncSimple.hpp
View File

@ -18,17 +18,21 @@
*/ */
#ifndef SATFUNCSIMPLE_HPP #ifndef SATFUNCSIMPLE_HPP
#define SATFUNCSIMPLE_HPP #define SATFUNCSIMPLE_HPP
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp> #include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp> #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <vector> #include <vector>
namespace Opm namespace Opm
{ {
class SatFuncSimple: public BlackoilPhases class SatFuncSimpleUniform : public BlackoilPhases
{ {
public: public:
void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); void init(const EclipseGridParser& deck,
void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg, const int table_num,
const PhaseUsage phase_usg,
const int samples); const int samples);
void evalKr(const double* s, double* kr) const; void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
@ -46,5 +50,31 @@ namespace Opm
UniformTableLinear<double> pcog_; UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc) double krocw_; // = krow_(s_wc)
}; };
class SatFuncSimpleNonuniform : public BlackoilPhases
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
} // namespace Opm } // namespace Opm
#endif // SATFUNCSIMPLE_HPP #endif // SATFUNCSIMPLE_HPP

232
opm/core/fluid/SatFuncStone2.cpp
View File

@ -1,3 +1,22 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SatFuncStone2.hpp> #include <opm/core/fluid/SatFuncStone2.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp> #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
@ -6,20 +25,14 @@
#include <opm/core/utility/buildUniformMonotoneTable.hpp> #include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <iostream> #include <iostream>
namespace Opm namespace Opm
{ {
void SatFuncStone2::init(const EclipseGridParser& deck, void SatFuncStone2Uniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg)
{
init(deck, table_num, phase_usg, 200);
}
void SatFuncStone2::init(const EclipseGridParser& deck,
const int table_num, const int table_num,
const PhaseUsage phase_usg, const PhaseUsage phase_usg,
const int samples) const int samples)
@ -64,7 +77,7 @@ namespace Opm
} }
void SatFuncStone2::evalKr(const double* s, double* kr) const void SatFuncStone2Uniform::evalKr(const double* s, double* kr) const
{ {
if (phase_usage.num_phases == 3) { if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model. // Stone-II relative permeability model.
@ -105,7 +118,7 @@ namespace Opm
} }
void SatFuncStone2::evalKrDeriv(const double* s, double* kr, double* dkrds) const void SatFuncStone2Uniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{ {
const int np = phase_usage.num_phases; const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0); std::fill(dkrds, dkrds + np*np, 0.0);
@ -167,7 +180,7 @@ namespace Opm
} }
void SatFuncStone2::evalPc(const double* s, double* pc) const void SatFuncStone2Uniform::evalPc(const double* s, double* pc) const
{ {
pc[phase_usage.phase_pos[Liquid]] = 0.0; pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) { if (phase_usage.phase_used[Aqua]) {
@ -180,7 +193,7 @@ namespace Opm
} }
} }
void SatFuncStone2::evalPcDeriv(const double* s, double* pc, double* dpcds) const void SatFuncStone2Uniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{ {
// The problem of determining three-phase capillary pressures // The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase // is very hard experimentally, usually one extends two-phase
@ -206,4 +219,199 @@ namespace Opm
// ====== Methods for SatFuncSimpleNonuniform ======
void SatFuncStone2Nonuniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int /*samples*/)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
krw_ = NonuniformTableLinear<double>(sw, krw);
krow_ = NonuniformTableLinear<double>(sw, krow);
pcow_ = NonuniformTableLinear<double>(sw, pcow);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
krg_ = NonuniformTableLinear<double>(sg, krg);
krog_ = NonuniformTableLinear<double>(sg, krog);
pcog_ = NonuniformTableLinear<double>(sg, pcog);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncStone2Nonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
double krog = krog_(sg); // = 1 - so - sw
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncStone2Nonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncStone2Nonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncStone2Nonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
} // namespace Opm } // namespace Opm

36
opm/core/fluid/SatFuncStone2.hpp
View File

@ -18,17 +18,21 @@
*/ */
#ifndef SATFUNCSTONE2_HPP #ifndef SATFUNCSTONE2_HPP
#define SATFUNCSTONE2_HPP #define SATFUNCSTONE2_HPP
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp> #include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp> #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <vector> #include <vector>
namespace Opm namespace Opm
{ {
class SatFuncStone2: public BlackoilPhases class SatFuncStone2Uniform : public BlackoilPhases
{ {
public: public:
void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); void init(const EclipseGridParser& deck,
void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg, const int table_num,
const PhaseUsage phase_usg,
const int samples); const int samples);
void evalKr(const double* s, double* kr) const; void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
@ -46,5 +50,31 @@ namespace Opm
UniformTableLinear<double> pcog_; UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc) double krocw_; // = krow_(s_wc)
}; };
class SatFuncStone2Nonuniform : public BlackoilPhases
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
} // namespace Opm } // namespace Opm
#endif // SATFUNCSTONE2_HPP #endif // SATFUNCSTONE2_HPP

204
opm/core/fluid/SaturationPropsBasic.cpp
View File

@ -29,64 +29,64 @@ namespace Opm
namespace { namespace {
struct KrFunConstant struct KrFunConstant
{ {
double kr(double) double kr(double)
{ {
return 1.0; return 1.0;
} }
double dkrds(double) double dkrds(double)
{ {
return 0.0; return 0.0;
} }
}; };
struct KrFunLinear struct KrFunLinear
{ {
double kr(double s) double kr(double s)
{ {
return s; return s;
} }
double dkrds(double) double dkrds(double)
{ {
return 1.0; return 1.0;
} }
}; };
struct KrFunQuadratic struct KrFunQuadratic
{ {
double kr(double s) double kr(double s)
{ {
return s*s; return s*s;
} }
double dkrds(double s) double dkrds(double s)
{ {
return 2.0*s; return 2.0*s;
} }
}; };
template <class Fun> template <class Fun>
static inline void evalAllKrDeriv(const int n, const int np, static inline void evalAllKrDeriv(const int n, const int np,
const double* s, double* kr, double* dkrds, Fun fun) const double* s, double* kr, double* dkrds, Fun fun)
{ {
if (dkrds == 0) { if (dkrds == 0) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n*np; ++i) { for (int i = 0; i < n*np; ++i) {
kr[i] = fun.kr(s[i]); kr[i] = fun.kr(s[i]);
} }
return; return;
} }
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0);
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
kr[i*np + phase] = fun.kr(s[i*np + phase]); kr[i*np + phase] = fun.kr(s[i*np + phase]);
// Only diagonal elements in derivative. // Only diagonal elements in derivative.
dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]); dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]);
} }
} }
} }
} // anon namespace } // anon namespace
@ -109,25 +109,25 @@ namespace Opm
/// Initialize from parameters. /// Initialize from parameters.
void SaturationPropsBasic::init(const parameter::ParameterGroup& param) void SaturationPropsBasic::init(const parameter::ParameterGroup& param)
{ {
int num_phases = param.getDefault("num_phases", 2); int num_phases = param.getDefault("num_phases", 2);
if (num_phases > 2 || num_phases < 1) { if (num_phases > 2 || num_phases < 1) {
THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases); THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases);
} }
num_phases_ = num_phases; num_phases_ = num_phases;
//std::string rpf = param.getDefault("relperm_func", std::string("Unset")); //std::string rpf = param.getDefault("relperm_func", std::string("Unset"));
std::string rpf = param.getDefault("relperm_func", std::string("Linear")); std::string rpf = param.getDefault("relperm_func", std::string("Linear"));
if (rpf == "Constant") { if (rpf == "Constant") {
relperm_func_ = Constant; relperm_func_ = Constant;
if(num_phases!=1){ if(num_phases!=1){
THROW("Constant relperm with more than one phase???"); THROW("Constant relperm with more than one phase???");
} }
} else if (rpf == "Linear") { } else if (rpf == "Linear") {
relperm_func_ = Linear; relperm_func_ = Linear;
} else if (rpf == "Quadratic") { } else if (rpf == "Quadratic") {
relperm_func_ = Quadratic; relperm_func_ = Quadratic;
} else { } else {
THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf); THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf);
} }
} }
@ -136,7 +136,7 @@ namespace Opm
/// \return P, the number of phases. /// \return P, the number of phases.
int SaturationPropsBasic::numPhases() const int SaturationPropsBasic::numPhases() const
{ {
return num_phases_; return num_phases_;
} }
@ -152,29 +152,29 @@ namespace Opm
/// m_{ij} = \frac{dkr_i}{ds^j}, /// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsBasic::relperm(const int n, void SaturationPropsBasic::relperm(const int n,
const double* s, const double* s,
double* kr, double* kr,
double* dkrds) const double* dkrds) const
{ {
switch (relperm_func_) { switch (relperm_func_) {
case Constant: case Constant:
{ {
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant()); evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant());
break; break;
} }
case Linear: case Linear:
{ {
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear()); evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear());
break; break;
} }
case Quadratic: case Quadratic:
{ {
evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic()); evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic());
break; break;
} }
default: default:
THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_); THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_);
} }
} }
@ -190,13 +190,13 @@ namespace Opm
/// m_{ij} = \frac{dpc_i}{ds^j}, /// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...) /// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsBasic::capPress(const int n, void SaturationPropsBasic::capPress(const int n,
const double* /*s*/, const double* /*s*/,
double* pc, double* pc,
double* dpcds) const double* dpcds) const
{ {
std::fill(pc, pc + num_phases_*n, 0.0); std::fill(pc, pc + num_phases_*n, 0.0);
if (dpcds) { if (dpcds) {
std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0); std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0);
} }
} }
@ -207,11 +207,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void SaturationPropsBasic::satRange(const int n, void SaturationPropsBasic::satRange(const int n,
double* smin, double* smin,
double* smax) const double* smax) const
{ {
std::fill(smin, smin + num_phases_*n, 0.0); std::fill(smin, smin + num_phases_*n, 0.0);
std::fill(smax, smax + num_phases_*n, 1.0); std::fill(smax, smax + num_phases_*n, 1.0);
} }

20
opm/core/fluid/SaturationPropsBasic.hpp
View File

@ -40,12 +40,12 @@ namespace Opm
SaturationPropsBasic(); SaturationPropsBasic();
/// Initialize from parameters. /// Initialize from parameters.
/// The following parameters are accepted (defaults): /// The following parameters are accepted (defaults):
/// num_phases (2) Must be 1 or 2. /// num_phases (2) Must be 1 or 2.
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
void init(const parameter::ParameterGroup& param); void init(const parameter::ParameterGroup& param);
enum RelPermFunc { Constant, Linear, Quadratic }; enum RelPermFunc { Constant, Linear, Quadratic };
/// Initialize from arguments a basic Saturation property. /// Initialize from arguments a basic Saturation property.
void init(const int num_phases, void init(const int num_phases,
@ -86,18 +86,18 @@ namespace Opm
double* pc, double* pc,
double* dpcds) const; double* dpcds) const;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void satRange(const int n, void satRange(const int n,
double* smin, double* smin,
double* smax) const; double* smax) const;
private: private:
int num_phases_; int num_phases_;
RelPermFunc relperm_func_; RelPermFunc relperm_func_;
}; };

232
opm/core/fluid/SaturationPropsFromDeck.cpp
View File

@ -19,7 +19,6 @@
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp> #include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <iostream> #include <iostream>
@ -27,236 +26,9 @@
namespace Opm namespace Opm
{ {
/// Default constructor. // This file should be removed in the future.
SaturationPropsFromDeck::SaturationPropsFromDeck() // Holding off until refactoring of SaturationPropsFromDeck class is done.
{
}
/// Initialize from deck.
void SaturationPropsFromDeck::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
THROW("SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = grid.number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = grid.global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
cell_to_func_[cell] = satnum[deck_pos] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
num_tables = swof_table.size();
if (num_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
int num_sgof_tables = sgof_table.size();
if (num_sgof_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
THROW("Inconsistent number of tables in SWOF and SGOF.");
}
}
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_);
}
}
/// Initialize from deck, grid and parameters
void SaturationPropsFromDeck::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
THROW("SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = grid.number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = grid.global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
cell_to_func_[cell] = satnum[deck_pos] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
num_tables = swof_table.size();
if (num_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
int num_sgof_tables = sgof_table.size();
if (num_sgof_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
THROW("Inconsistent number of tables in SWOF and SGOF.");
}
}
// Initialize tables.
const int tab_size = param.getDefault("tab_size_kr", 200);
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, tab_size);
}
}
/// \return P, the number of phases.
int SaturationPropsFromDeck::numPhases() const
{
return phase_usage_.num_phases;
}
/// Relative permeability.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsFromDeck::relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
}
}
/// Capillary pressure.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsFromDeck::capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
if (dpcds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
}
}
}
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void SaturationPropsFromDeck::satRange(const int n,
const int* cells,
double* smin,
double* smax) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
for (int i = 0; i < n; ++i) {
for (int p = 0; p < np; ++p) {
smin[np*i + p] = funcForCell(cells[i]).smin_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
}
}
}
// Map the cell number to the correct function set.
const SaturationPropsFromDeck::satfunc_t&
SaturationPropsFromDeck::funcForCell(const int cell) const
{
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
} // namespace Opm } // namespace Opm

41
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

@ -19,9 +19,10 @@
#ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
#define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
#include <opm/core/fluid/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp> #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SatFuncStone2.hpp> #include <opm/core/fluid/SatFuncStone2.hpp>
#include <opm/core/fluid/SatFuncSimple.hpp> #include <opm/core/fluid/SatFuncSimple.hpp>
@ -32,23 +33,31 @@ struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
class SaturationPropsFromDeck : public BlackoilPhases
/// Interface to saturation functions from deck.
/// Possible values for template argument (for now):
/// SatFuncSetStone2Nonuniform,
/// SatFuncSetStone2Uniform.
/// SatFuncSetSimpleNonuniform,
/// SatFuncSetSimpleUniform.
template <class SatFuncSet>
class SaturationPropsFromDeck : public SaturationPropsInterface
{ {
public: public:
/// Default constructor. /// Default constructor.
SaturationPropsFromDeck(); SaturationPropsFromDeck();
/// Initialize from deck and grid. /// Initialize from deck and grid.
/// \param deck Deck input parser /// \param[in] deck Deck input parser
/// \param grid Grid to which property object applies, needed for the /// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid) /// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck. /// to logical cartesian indices consistent with the deck.
void init(const EclipseGridParser& deck, /// \param[in] samples Number of uniform sample points for saturation tables.
const UnstructuredGrid& grid); /// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(const EclipseGridParser& deck, void init(const EclipseGridParser& deck,
const UnstructuredGrid& grid, const UnstructuredGrid& grid,
const parameter::ParameterGroup& param); const int samples);
/// \return P, the number of phases. /// \return P, the number of phases.
int numPhases() const; int numPhases() const;
@ -83,22 +92,23 @@ namespace Opm
double* pc, double* pc,
double* dpcds) const; double* dpcds) const;
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling. /// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling. /// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void satRange(const int n, void satRange(const int n,
const int* cells, const int* cells,
double* smin, double* smin,
double* smax) const; double* smax) const;
private: private:
PhaseUsage phase_usage_; PhaseUsage phase_usage_;
typedef SatFuncSimple satfunc_t; std::vector<SatFuncSet> satfuncset_;
std::vector<satfunc_t> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1 std::vector<int> cell_to_func_; // = SATNUM - 1
const satfunc_t& funcForCell(const int cell) const; typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
}; };
@ -106,6 +116,7 @@ namespace Opm
} // namespace Opm } // namespace Opm
#include <opm/core/fluid/SaturationPropsFromDeck_impl.hpp>
#endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED

221
opm/core/fluid/SaturationPropsFromDeck_impl.hpp Normal file
View File

@ -0,0 +1,221 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
#define OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/grid.h>
namespace Opm
{
// ----------- Methods of SaturationPropsFromDeck ---------
/// Default constructor.
template <class SatFuncSet>
SaturationPropsFromDeck<SatFuncSet>::SaturationPropsFromDeck()
{
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const int samples)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
THROW("SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = grid.number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = grid.global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
cell_to_func_[cell] = satnum[deck_pos] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
num_tables = swof_table.size();
if (num_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
int num_sgof_tables = sgof_table.size();
if (num_sgof_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
THROW("Inconsistent number of tables in SWOF and SGOF.");
}
}
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
}
}
/// \return P, the number of phases.
template <class SatFuncSet>
int SaturationPropsFromDeck<SatFuncSet>::numPhases() const
{
return phase_usage_.num_phases;
}
/// Relative permeability.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
}
}
/// Capillary pressure.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
if (dpcds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
}
}
}
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::satRange(const int n,
const int* cells,
double* smin,
double* smax) const
{
ASSERT (cells != 0);
const int np = phase_usage_.num_phases;
for (int i = 0; i < n; ++i) {
for (int p = 0; p < np; ++p) {
smin[np*i + p] = funcForCell(cells[i]).smin_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
}
}
}
// Map the cell number to the correct function set.
template <class SatFuncSet>
const typename SaturationPropsFromDeck<SatFuncSet>::Funcs&
SaturationPropsFromDeck<SatFuncSet>::funcForCell(const int cell) const
{
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
} // namespace Opm
#endif // OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED

85
opm/core/fluid/SaturationPropsInterface.hpp Normal file
View File

@ -0,0 +1,85 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
#define OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
namespace Opm
{
class SaturationPropsInterface : public BlackoilPhases
{
public:
/// Virtual destructor.
virtual ~SaturationPropsInterface() {};
/// \return P, the number of phases.
virtual int numPhases() const = 0;
/// Relative permeability.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
virtual void relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const = 0;
/// Capillary pressure.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
/// array must be valid before calling.
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
virtual void capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const = 0;
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
virtual void satRange(const int n,
const int* cells,
double* smin,
double* smax) const = 0;
};
} // namespace Opm
#endif // OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED

31
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp
View File

@ -21,6 +21,7 @@
#include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp> #include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
#include <opm/core/fluid/blackoil/SinglePvtDead.hpp> #include <opm/core/fluid/blackoil/SinglePvtDead.hpp>
#include <opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp>
#include <opm/core/fluid/blackoil/SinglePvtLiveOil.hpp> #include <opm/core/fluid/blackoil/SinglePvtLiveOil.hpp>
#include <opm/core/fluid/blackoil/SinglePvtLiveGas.hpp> #include <opm/core/fluid/blackoil/SinglePvtLiveGas.hpp>
#include <opm/core/fluid/blackoil/SinglePvtConstCompr.hpp> #include <opm/core/fluid/blackoil/SinglePvtConstCompr.hpp>
@ -43,14 +44,14 @@ namespace Opm
{ {
typedef std::vector<std::vector<std::vector<double> > > table_t; typedef std::vector<std::vector<std::vector<double> > > table_t;
// If we need multiple regions, this class and the SinglePvt* classes must change. // If we need multiple regions, this class and the SinglePvt* classes must change.
region_number_ = 0; region_number_ = 0;
phase_usage_ = phaseUsageFromDeck(deck); phase_usage_ = phaseUsageFromDeck(deck);
// Surface densities. Accounting for different orders in eclipse and our code. // Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) { if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number_]; const std::vector<double>& d = deck.getDENSITY().densities_[region_number_];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
if (phase_usage_.phase_used[Aqua]) { if (phase_usage_.phase_used[Aqua]) {
densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water]; densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water];
} }
@ -60,9 +61,9 @@ namespace Opm
if (phase_usage_.phase_used[Liquid]) { if (phase_usage_.phase_used[Liquid]) {
densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil]; densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil];
} }
} else { } else {
THROW("Input is missing DENSITY\n"); THROW("Input is missing DENSITY\n");
} }
// Set the properties. // Set the properties.
props_.resize(phase_usage_.num_phases); props_.resize(phase_usage_.num_phases);
@ -78,7 +79,11 @@ namespace Opm
// Oil PVT // Oil PVT
if (phase_usage_.phase_used[Liquid]) { if (phase_usage_.phase_used[Liquid]) {
if (deck.hasField("PVDO")) { if (deck.hasField("PVDO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_, samples)); if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_, samples));
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_));
}
} else if (deck.hasField("PVTO")) { } else if (deck.hasField("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_)); props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_));
} else if (deck.hasField("PVCDO")) { } else if (deck.hasField("PVCDO")) {
@ -87,10 +92,14 @@ namespace Opm
THROW("Input is missing PVDO or PVTO\n"); THROW("Input is missing PVDO or PVTO\n");
} }
} }
// Gas PVT // Gas PVT
if (phase_usage_.phase_used[Vapour]) { if (phase_usage_.phase_used[Vapour]) {
if (deck.hasField("PVDG")) { if (deck.hasField("PVDG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_, samples)); if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_, samples));
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_));
}
} else if (deck.hasField("PVTG")) { } else if (deck.hasField("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_)); props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_));
} else { } else {

14
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp
View File

@ -47,7 +47,13 @@ namespace Opm
BlackoilPvtProperties(); BlackoilPvtProperties();
/// Initialize from deck. /// Initialize from deck.
void init(const EclipseGridParser& deck, const int samples = 16); /// \param deck An input deck.
/// \param samples If greater than zero, indicates the number of
/// uniform samples to be taken from monotone spline
/// curves interpolating the fluid data.
/// Otherwise, interpolate linearly in the original
/// data without fitting a spline.
void init(const EclipseGridParser& deck, const int samples);
/// Number of active phases. /// Number of active phases.
int numPhases() const; int numPhases() const;
@ -64,7 +70,7 @@ namespace Opm
/// Densities of stock components at surface conditions. /// Densities of stock components at surface conditions.
/// \return Array of size numPhases(). /// \return Array of size numPhases().
const double* surfaceDensities() const; const double* surfaceDensities() const;
/// Viscosity as a function of p and z. /// Viscosity as a function of p and z.
void mu(const int n, void mu(const int n,
@ -105,11 +111,11 @@ namespace Opm
PhaseUsage phase_usage_; PhaseUsage phase_usage_;
int region_number_; int region_number_;
std::vector<std::tr1::shared_ptr<SinglePvtInterface> > props_; std::vector<std::tr1::shared_ptr<SinglePvtInterface> > props_;
double densities_[MaxNumPhases]; double densities_[MaxNumPhases];
mutable std::vector<double> data1_; mutable std::vector<double> data1_;
mutable std::vector<double> data2_; mutable std::vector<double> data2_;
}; };

14
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp
View File

@ -40,12 +40,12 @@ namespace Opm
public: public:
typedef std::vector<std::vector<double> > table_t; typedef std::vector<std::vector<double> > table_t;
SinglePvtConstCompr(const table_t& pvtw) SinglePvtConstCompr(const table_t& pvtw)
{ {
const int region_number = 0; const int region_number = 0;
if (pvtw.size() != 1) { if (pvtw.size() != 1) {
THROW("More than one PVD-region"); THROW("More than one PVD-region");
} }
ref_press_ = pvtw[region_number][0]; ref_press_ = pvtw[region_number][0];
ref_B_ = pvtw[region_number][1]; ref_B_ = pvtw[region_number][1];
comp_ = pvtw[region_number][2]; comp_ = pvtw[region_number][2];
@ -53,7 +53,7 @@ namespace Opm
visc_comp_ = pvtw[region_number][4]; visc_comp_ = pvtw[region_number][4];
} }
SinglePvtConstCompr(double visc) SinglePvtConstCompr(double visc)
: ref_press_(0.0), : ref_press_(0.0),
ref_B_(1.0), ref_B_(1.0),
comp_(0.0), comp_(0.0),
@ -62,7 +62,7 @@ namespace Opm
{ {
} }
virtual ~SinglePvtConstCompr() virtual ~SinglePvtConstCompr()
{ {
} }

26
opm/core/fluid/blackoil/SinglePvtDead.cpp
View File

@ -1,5 +1,5 @@
/* /*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <opm/core/fluid/blackoil/SinglePvtDead.hpp> #include <opm/core/fluid/blackoil/SinglePvtDead.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <algorithm> #include <algorithm>
// Extra includes for debug dumping of tables. // Extra includes for debug dumping of tables.
@ -33,25 +33,25 @@ namespace Opm
// Member functions // Member functions
//------------------------------------------------------------------------- //-------------------------------------------------------------------------
/// Constructor /// Constructor
SinglePvtDead::SinglePvtDead(const table_t& pvd_table, const int samples) SinglePvtDead::SinglePvtDead(const table_t& pvd_table)
{ {
const int region_number = 0; const int region_number = 0;
if (pvd_table.size() != 1) { if (pvd_table.size() != 1) {
THROW("More than one PVT-region"); THROW("More than one PVT-region");
} }
// Copy data // Copy data
const int sz = pvd_table[region_number][0].size(); const int sz = pvd_table[region_number][0].size();
std::vector<double> press(sz); std::vector<double> press(sz);
std::vector<double> B_inv(sz); std::vector<double> B_inv(sz);
std::vector<double> visc(sz); std::vector<double> visc(sz);
for (int i = 0; i < sz; ++i) { for (int i = 0; i < sz; ++i) {
press[i] = pvd_table[region_number][0][i]; press[i] = pvd_table[region_number][0][i];
B_inv[i] = 1.0 / pvd_table[region_number][1][i]; B_inv[i] = 1.0 / pvd_table[region_number][1][i];
visc[i] = pvd_table[region_number][2][i]; visc[i] = pvd_table[region_number][2][i];
} }
buildUniformMonotoneTable(press, B_inv, samples, one_over_B_); one_over_B_ = NonuniformTableLinear<double>(press, B_inv);
buildUniformMonotoneTable(press, visc, samples, viscosity_); viscosity_ = NonuniformTableLinear<double>(press, visc);
// Dumping the created tables. // Dumping the created tables.
// static int count = 0; // static int count = 0;

18
opm/core/fluid/blackoil/SinglePvtDead.hpp
View File

@ -1,5 +1,5 @@
/* /*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -22,7 +22,7 @@
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp> #include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp> #include <opm/core/utility/NonuniformTableLinear.hpp>
#include <vector> #include <vector>
namespace Opm namespace Opm
@ -36,9 +36,9 @@ namespace Opm
class SinglePvtDead : public SinglePvtInterface class SinglePvtDead : public SinglePvtInterface
{ {
public: public:
typedef std::vector<std::vector<std::vector<double> > > table_t; typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDead(const table_t& pvd_table, const int samples = 16); SinglePvtDead(const table_t& pvd_table);
virtual ~SinglePvtDead(); virtual ~SinglePvtDead();
/// Viscosity as a function of p and z. /// Viscosity as a function of p and z.
virtual void mu(const int n, virtual void mu(const int n,
@ -72,12 +72,12 @@ namespace Opm
double* output_R, double* output_R,
double* output_dRdp) const; double* output_dRdp) const;
private: private:
// PVT properties of dry gas or dead oil // PVT properties of dry gas or dead oil
UniformTableLinear<double> one_over_B_; NonuniformTableLinear<double> one_over_B_;
UniformTableLinear<double> viscosity_; NonuniformTableLinear<double> viscosity_;
}; };
} }
#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED
#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED

127
opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp Normal file
View File

@ -0,0 +1,127 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.
// #include <boost/lexical_cast.hpp>
// #include <string>
// #include <fstream>
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtDeadSpline::SinglePvtDeadSpline(const table_t& pvd_table, const int samples)
{
const int region_number = 0;
if (pvd_table.size() != 1) {
THROW("More than one PVT-region");
}
// Copy data
const int sz = pvd_table[region_number][0].size();
std::vector<double> press(sz);
std::vector<double> B_inv(sz);
std::vector<double> visc(sz);
for (int i = 0; i < sz; ++i) {
press[i] = pvd_table[region_number][0][i];
B_inv[i] = 1.0 / pvd_table[region_number][1][i];
visc[i] = pvd_table[region_number][2][i];
}
buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
buildUniformMonotoneTable(press, visc, samples, viscosity_);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
// Destructor
SinglePvtDeadSpline::~SinglePvtDeadSpline()
{
}
void SinglePvtDeadSpline::mu(const int n,
const double* p,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
}
}
void SinglePvtDeadSpline::B(const int n,
const double* p,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = 1.0/one_over_B_(p[i]);
}
}
void SinglePvtDeadSpline::dBdp(const int n,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, p, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*one_over_B_.derivative(p[i]);
}
}
void SinglePvtDeadSpline::R(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
void SinglePvtDeadSpline::dRdp(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
}

84
opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp Normal file
View File

@ -0,0 +1,84 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <vector>
namespace Opm
{
/// Class for immiscible dead oil and dry gas.
/// For all the virtual methods, the following apply: p and z
/// are expected to be of size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class SinglePvtDeadSpline : public SinglePvtInterface
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDeadSpline(const table_t& pvd_table, const int samples);
virtual ~SinglePvtDeadSpline();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const double* p,
const double* z,
double* output_mu) const;
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const double* p,
const double* z,
double* output_B,
double* output_dBdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
UniformTableLinear<double> one_over_B_;
UniformTableLinear<double> viscosity_;
};
}
#endif // OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED

2
opm/core/fluid/blackoil/SinglePvtInterface.hpp
View File

@ -32,7 +32,7 @@ namespace Opm
public: public:
SinglePvtInterface(); SinglePvtInterface();
virtual ~SinglePvtInterface(); virtual ~SinglePvtInterface();
/// \param[in] num_phases The number of active phases. /// \param[in] num_phases The number of active phases.
/// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases /// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases

236
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp
View File

@ -47,37 +47,37 @@ namespace Opm
/// Constructor /// Constructor
SinglePvtLiveGas::SinglePvtLiveGas(const table_t& pvtg) SinglePvtLiveGas::SinglePvtLiveGas(const table_t& pvtg)
{ {
// GAS, PVTG // GAS, PVTG
const int region_number = 0; const int region_number = 0;
if (pvtg.size() != 1) { if (pvtg.size() != 1) {
THROW("More than one PVD-region"); THROW("More than one PVD-region");
} }
saturated_gas_table_.resize(4); saturated_gas_table_.resize(4);
const int sz = pvtg[region_number].size(); const int sz = pvtg[region_number].size();
for (int k=0; k<4; ++k) { for (int k=0; k<4; ++k) {
saturated_gas_table_[k].resize(sz); saturated_gas_table_[k].resize(sz);
} }
for (int i=0; i<sz; ++i) { for (int i=0; i<sz; ++i) {
saturated_gas_table_[0][i] = pvtg[region_number][i][0]; // p saturated_gas_table_[0][i] = pvtg[region_number][i][0]; // p
saturated_gas_table_[1][i] = pvtg[region_number][i][2]; // Bg saturated_gas_table_[1][i] = pvtg[region_number][i][2]; // Bg
saturated_gas_table_[2][i] = pvtg[region_number][i][3]; // mu_g saturated_gas_table_[2][i] = pvtg[region_number][i][3]; // mu_g
saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv
} }
undersat_gas_tables_.resize(sz); undersat_gas_tables_.resize(sz);
for (int i=0; i<sz; ++i) { for (int i=0; i<sz; ++i) {
undersat_gas_tables_[i].resize(3); undersat_gas_tables_[i].resize(3);
int tsize = (pvtg[region_number][i].size() - 1)/3; int tsize = (pvtg[region_number][i].size() - 1)/3;
undersat_gas_tables_[i][0].resize(tsize); undersat_gas_tables_[i][0].resize(tsize);
undersat_gas_tables_[i][1].resize(tsize); undersat_gas_tables_[i][1].resize(tsize);
undersat_gas_tables_[i][2].resize(tsize); undersat_gas_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) { for (int j=0, k=0; j<tsize; ++j) {
undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv
undersat_gas_tables_[i][1][j] = pvtg[region_number][i][++k]; // Bg undersat_gas_tables_[i][1][j] = pvtg[region_number][i][++k]; // Bg
undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]; // mu_g undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]; // mu_g
} }
} }
} }
// Destructor // Destructor
@ -184,16 +184,16 @@ namespace Opm
// To handle no-gas case. // To handle no-gas case.
return 0.0; return 0.0;
} }
double satR = linearInterpolationExtrap(saturated_gas_table_[0], double satR = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press); saturated_gas_table_[3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) { if (satR < maxR ) {
// Saturated case // Saturated case
return satR; return satR;
} else { } else {
// Undersaturated case // Undersaturated case
return maxR; return maxR;
} }
} }
void SinglePvtLiveGas::evalRDeriv(const double press, const double* surfvol, void SinglePvtLiveGas::evalRDeriv(const double press, const double* surfvol,
@ -205,20 +205,20 @@ namespace Opm
dRdpval = 0.0; dRdpval = 0.0;
return; return;
} }
double satR = linearInterpolationExtrap(saturated_gas_table_[0], double satR = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press); saturated_gas_table_[3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) { if (satR < maxR ) {
// Saturated case // Saturated case
Rval = satR; Rval = satR;
dRdpval = linearInterpolDerivative(saturated_gas_table_[0], dRdpval = linearInterpolDerivative(saturated_gas_table_[0],
saturated_gas_table_[3], saturated_gas_table_[3],
press); press);
} else { } else {
// Undersaturated case // Undersaturated case
Rval = maxR; Rval = maxR;
dRdpval = 0.0; dRdpval = 0.0;
} }
} }
double SinglePvtLiveGas::miscible_gas(const double press, double SinglePvtLiveGas::miscible_gas(const double press,
@ -226,81 +226,81 @@ namespace Opm
const int item, const int item,
const bool deriv) const const bool deriv) const
{ {
int section; int section;
double Rval = linearInterpolationExtrap(saturated_gas_table_[0], double Rval = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press, saturated_gas_table_[3], press,
section); section);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (deriv) { if (deriv) {
if (Rval < maxR ) { // Saturated case if (Rval < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_gas_table_[0], return linearInterpolDerivative(saturated_gas_table_[0],
saturated_gas_table_[item], saturated_gas_table_[item],
press); press);
} else { // Undersaturated case } else { // Undersaturated case
int is = section; int is = section;
if (undersat_gas_tables_[is][0].size() < 2) { if (undersat_gas_tables_[is][0].size() < 2) {
double val = (saturated_gas_table_[item][is+1] double val = (saturated_gas_table_[item][is+1]
- saturated_gas_table_[item][is]) / - saturated_gas_table_[item][is]) /
(saturated_gas_table_[0][is+1] - (saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]); saturated_gas_table_[0][is]);
return val; return val;
} }
double val1 = double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0], linearInterpolationExtrap(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item], undersat_gas_tables_[is][item],
maxR); maxR);
double val2 = double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0], linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item], undersat_gas_tables_[is+1][item],
maxR); maxR);
double val = (val2 - val1)/ double val = (val2 - val1)/
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]); (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
return val; return val;
} }
} else { } else {
if (Rval < maxR ) { // Saturated case if (Rval < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_gas_table_[0], return linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[item], saturated_gas_table_[item],
press); press);
} else { // Undersaturated case } else { // Undersaturated case
int is = section; int is = section;
// Extrapolate from first table section // Extrapolate from first table section
if (is == 0 && press < saturated_gas_table_[0][0]) { if (is == 0 && press < saturated_gas_table_[0][0]) {
return linearInterpolationExtrap(undersat_gas_tables_[0][0], return linearInterpolationExtrap(undersat_gas_tables_[0][0],
undersat_gas_tables_[0][item], undersat_gas_tables_[0][item],
maxR); maxR);
} }
// Extrapolate from last table section // Extrapolate from last table section
int ltp = saturated_gas_table_[0].size() - 1; int ltp = saturated_gas_table_[0].size() - 1;
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) { if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
return linearInterpolationExtrap(undersat_gas_tables_[ltp][0], return linearInterpolationExtrap(undersat_gas_tables_[ltp][0],
undersat_gas_tables_[ltp][item], undersat_gas_tables_[ltp][item],
maxR); maxR);
} }
// Interpolate between table sections // Interpolate between table sections
double w = (press - saturated_gas_table_[0][is]) / double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] - (saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]); saturated_gas_table_[0][is]);
if (undersat_gas_tables_[is][0].size() < 2) { if (undersat_gas_tables_[is][0].size() < 2) {
double val = saturated_gas_table_[item][is] + double val = saturated_gas_table_[item][is] +
w*(saturated_gas_table_[item][is+1] - w*(saturated_gas_table_[item][is+1] -
saturated_gas_table_[item][is]); saturated_gas_table_[item][is]);
return val; return val;
} }
double val1 = double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0], linearInterpolationExtrap(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item], undersat_gas_tables_[is][item],
maxR); maxR);
double val2 = double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0], linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item], undersat_gas_tables_[is+1][item],
maxR); maxR);
double val = val1 + w*(val2 - val1); double val = val1 + w*(val2 - val1);
return val; return val;
} }
} }
} }

18
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp
View File

@ -33,10 +33,10 @@ namespace Opm
class SinglePvtLiveGas : public SinglePvtInterface class SinglePvtLiveGas : public SinglePvtInterface
{ {
public: public:
typedef std::vector<std::vector<std::vector<double> > > table_t; typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtLiveGas(const table_t& pvto); SinglePvtLiveGas(const table_t& pvto);
virtual ~SinglePvtLiveGas(); virtual ~SinglePvtLiveGas();
/// Viscosity as a function of p and z. /// Viscosity as a function of p and z.
virtual void mu(const int n, virtual void mu(const int n,
@ -76,14 +76,14 @@ namespace Opm
double evalR(double press, const double* surfvol) const; double evalR(double press, const double* surfvol) const;
void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const; void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
// item: 1=>B 2=>mu; // item: 1=>B 2=>mu;
double miscible_gas(const double press, double miscible_gas(const double press,
const double* surfvol, const double* surfvol,
const int item, const int item,
const bool deriv = false) const; const bool deriv = false) const;
// PVT properties of wet gas (with vaporised oil) // PVT properties of wet gas (with vaporised oil)
std::vector<std::vector<double> > saturated_gas_table_; std::vector<std::vector<double> > saturated_gas_table_;
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_; std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
}; };

356
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp
View File

@ -38,122 +38,122 @@ namespace Opm
/// Constructor /// Constructor
SinglePvtLiveOil::SinglePvtLiveOil(const table_t& pvto) SinglePvtLiveOil::SinglePvtLiveOil(const table_t& pvto)
{ {
// OIL, PVTO // OIL, PVTO
const int region_number = 0; const int region_number = 0;
if (pvto.size() != 1) { if (pvto.size() != 1) {
THROW("More than one PVD-region"); THROW("More than one PVD-region");
} }
saturated_oil_table_.resize(4); saturated_oil_table_.resize(4);
const int sz = pvto[region_number].size(); const int sz = pvto[region_number].size();
for (int k=0; k<4; ++k) { for (int k=0; k<4; ++k) {
saturated_oil_table_[k].resize(sz); saturated_oil_table_[k].resize(sz);
} }
for (int i=0; i<sz; ++i) { for (int i=0; i<sz; ++i) {
saturated_oil_table_[0][i] = pvto[region_number][i][1]; // p saturated_oil_table_[0][i] = pvto[region_number][i][1]; // p
saturated_oil_table_[1][i] = 1.0/pvto[region_number][i][2]; // 1/Bo saturated_oil_table_[1][i] = 1.0/pvto[region_number][i][2]; // 1/Bo
saturated_oil_table_[2][i] = pvto[region_number][i][3]; // mu_o saturated_oil_table_[2][i] = pvto[region_number][i][3]; // mu_o
saturated_oil_table_[3][i] = pvto[region_number][i][0]; // Rs saturated_oil_table_[3][i] = pvto[region_number][i][0]; // Rs
} }
undersat_oil_tables_.resize(sz); undersat_oil_tables_.resize(sz);
for (int i=0; i<sz; ++i) { for (int i=0; i<sz; ++i) {
undersat_oil_tables_[i].resize(3); undersat_oil_tables_[i].resize(3);
int tsize = (pvto[region_number][i].size() - 1)/3; int tsize = (pvto[region_number][i].size() - 1)/3;
undersat_oil_tables_[i][0].resize(tsize); undersat_oil_tables_[i][0].resize(tsize);
undersat_oil_tables_[i][1].resize(tsize); undersat_oil_tables_[i][1].resize(tsize);
undersat_oil_tables_[i][2].resize(tsize); undersat_oil_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) { for (int j=0, k=0; j<tsize; ++j) {
undersat_oil_tables_[i][0][j] = pvto[region_number][i][++k]; // p undersat_oil_tables_[i][0][j] = pvto[region_number][i][++k]; // p
undersat_oil_tables_[i][1][j] = 1.0/pvto[region_number][i][++k]; // 1/Bo undersat_oil_tables_[i][1][j] = 1.0/pvto[region_number][i][++k]; // 1/Bo
undersat_oil_tables_[i][2][j] = pvto[region_number][i][++k]; // mu_o undersat_oil_tables_[i][2][j] = pvto[region_number][i][++k]; // mu_o
} }
} }
// Fill in additional entries in undersaturated tables by interpolating/extrapolating 1/Bo and mu_o ... // Fill in additional entries in undersaturated tables by interpolating/extrapolating 1/Bo and mu_o ...
int iPrev = -1; int iPrev = -1;
int iNext = 1; int iNext = 1;
while (undersat_oil_tables_[iNext][0].size() < 2) { while (undersat_oil_tables_[iNext][0].size() < 2) {
++iNext; ++iNext;
} }
ASSERT(iNext < sz); ASSERT(iNext < sz);
for (int i=0; i<sz; ++i) { for (int i=0; i<sz; ++i) {
if (undersat_oil_tables_[i][0].size() > 1) { if (undersat_oil_tables_[i][0].size() > 1) {
iPrev = i; iPrev = i;
continue; continue;
} }
bool flagPrev = (iPrev >= 0); bool flagPrev = (iPrev >= 0);
bool flagNext = true; bool flagNext = true;
if (iNext < i) { if (iNext < i) {
iPrev = iNext; iPrev = iNext;
flagPrev = true; flagPrev = true;
iNext = i+1; iNext = i+1;
while (undersat_oil_tables_[iNext][0].size() < 2) { while (undersat_oil_tables_[iNext][0].size() < 2) {
++iNext; ++iNext;
} }
} }
double slopePrevBinv = 0.0; double slopePrevBinv = 0.0;
double slopePrevVisc = 0.0; double slopePrevVisc = 0.0;
double slopeNextBinv = 0.0; double slopeNextBinv = 0.0;
double slopeNextVisc = 0.0; double slopeNextVisc = 0.0;
while (flagPrev || flagNext) { while (flagPrev || flagNext) {
double pressure0 = undersat_oil_tables_[i][0].back(); double pressure0 = undersat_oil_tables_[i][0].back();
double pressure = 1.0e47; double pressure = 1.0e47;
if (flagPrev) { if (flagPrev) {
std::vector<double>::iterator itPrev = upper_bound(undersat_oil_tables_[iPrev][0].begin(), std::vector<double>::iterator itPrev = upper_bound(undersat_oil_tables_[iPrev][0].begin(),
undersat_oil_tables_[iPrev][0].end(),pressure0+1.); undersat_oil_tables_[iPrev][0].end(),pressure0+1.);
if (itPrev == undersat_oil_tables_[iPrev][0].end()) { if (itPrev == undersat_oil_tables_[iPrev][0].end()) {
--itPrev; // Extrapolation ... --itPrev; // Extrapolation ...
} else if (itPrev == undersat_oil_tables_[iPrev][0].begin()) { } else if (itPrev == undersat_oil_tables_[iPrev][0].begin()) {
++itPrev; ++itPrev;
} }
if (itPrev == undersat_oil_tables_[iPrev][0].end()-1) { if (itPrev == undersat_oil_tables_[iPrev][0].end()-1) {
flagPrev = false; // Last data set for "prev" ... flagPrev = false; // Last data set for "prev" ...
} }
double dPPrev = *itPrev - *(itPrev-1); double dPPrev = *itPrev - *(itPrev-1);
pressure = *itPrev; pressure = *itPrev;
int index = int(itPrev - undersat_oil_tables_[iPrev][0].begin()); int index = int(itPrev - undersat_oil_tables_[iPrev][0].begin());
slopePrevBinv = (undersat_oil_tables_[iPrev][1][index] - undersat_oil_tables_[iPrev][1][index-1])/dPPrev; slopePrevBinv = (undersat_oil_tables_[iPrev][1][index] - undersat_oil_tables_[iPrev][1][index-1])/dPPrev;
slopePrevVisc = (undersat_oil_tables_[iPrev][2][index] - undersat_oil_tables_[iPrev][2][index-1])/dPPrev; slopePrevVisc = (undersat_oil_tables_[iPrev][2][index] - undersat_oil_tables_[iPrev][2][index-1])/dPPrev;
} }
if (flagNext) { if (flagNext) {
std::vector<double>::iterator itNext = upper_bound(undersat_oil_tables_[iNext][0].begin(), std::vector<double>::iterator itNext = upper_bound(undersat_oil_tables_[iNext][0].begin(),
undersat_oil_tables_[iNext][0].end(),pressure0+1.); undersat_oil_tables_[iNext][0].end(),pressure0+1.);
if (itNext == undersat_oil_tables_[iNext][0].end()) { if (itNext == undersat_oil_tables_[iNext][0].end()) {
--itNext; // Extrapolation ... --itNext; // Extrapolation ...
} else if (itNext == undersat_oil_tables_[iNext][0].begin()) { } else if (itNext == undersat_oil_tables_[iNext][0].begin()) {
++itNext; ++itNext;
} }
if (itNext == undersat_oil_tables_[iNext][0].end()-1) { if (itNext == undersat_oil_tables_[iNext][0].end()-1) {
flagNext = false; // Last data set for "next" ... flagNext = false; // Last data set for "next" ...
} }
double dPNext = *itNext - *(itNext-1); double dPNext = *itNext - *(itNext-1);
if (flagPrev) { if (flagPrev) {
pressure = std::min(pressure,*itNext); pressure = std::min(pressure,*itNext);
} else { } else {
pressure = *itNext; pressure = *itNext;
} }
int index = int(itNext - undersat_oil_tables_[iNext][0].begin()); int index = int(itNext - undersat_oil_tables_[iNext][0].begin());
slopeNextBinv = (undersat_oil_tables_[iNext][1][index] - undersat_oil_tables_[iNext][1][index-1])/dPNext; slopeNextBinv = (undersat_oil_tables_[iNext][1][index] - undersat_oil_tables_[iNext][1][index-1])/dPNext;
slopeNextVisc = (undersat_oil_tables_[iNext][2][index] - undersat_oil_tables_[iNext][2][index-1])/dPNext; slopeNextVisc = (undersat_oil_tables_[iNext][2][index] - undersat_oil_tables_[iNext][2][index-1])/dPNext;
} }
double dP = pressure - pressure0; double dP = pressure - pressure0;
if (iPrev >= 0) { if (iPrev >= 0) {
double w = (saturated_oil_table_[3][i] - saturated_oil_table_[3][iPrev]) / double w = (saturated_oil_table_[3][i] - saturated_oil_table_[3][iPrev]) /
(saturated_oil_table_[3][iNext] - saturated_oil_table_[3][iPrev]); (saturated_oil_table_[3][iNext] - saturated_oil_table_[3][iPrev]);
undersat_oil_tables_[i][0].push_back(pressure0+dP); undersat_oil_tables_[i][0].push_back(pressure0+dP);
undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back() + undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back() +
dP*(slopePrevBinv+w*(slopeNextBinv-slopePrevBinv))); dP*(slopePrevBinv+w*(slopeNextBinv-slopePrevBinv)));
undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back() + undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back() +
dP*(slopePrevVisc+w*(slopeNextVisc-slopePrevVisc))); dP*(slopePrevVisc+w*(slopeNextVisc-slopePrevVisc)));
} else { } else {
undersat_oil_tables_[i][0].push_back(pressure0+dP); undersat_oil_tables_[i][0].push_back(pressure0+dP);
undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back()+dP*slopeNextBinv); undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back()+dP*slopeNextBinv);
undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back()+dP*slopeNextVisc); undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back()+dP*slopeNextVisc);
} }
} }
} }
} }
/// Destructor. /// Destructor.
@ -238,15 +238,15 @@ namespace Opm
double SinglePvtLiveOil::evalB(double press, const double* surfvol) const double SinglePvtLiveOil::evalB(double press, const double* surfvol) const
{ {
// if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case. // if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
return 1.0/miscible_oil(press, surfvol, 1, false); return 1.0/miscible_oil(press, surfvol, 1, false);
} }
void SinglePvtLiveOil::evalBDeriv(const double press, const double* surfvol, void SinglePvtLiveOil::evalBDeriv(const double press, const double* surfvol,
double& Bval, double& dBdpval) const double& Bval, double& dBdpval) const
{ {
Bval = evalB(press, surfvol); Bval = evalB(press, surfvol);
dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true); dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true);
} }
double SinglePvtLiveOil::evalR(double press, const double* surfvol) const double SinglePvtLiveOil::evalR(double press, const double* surfvol) const
@ -254,14 +254,14 @@ namespace Opm
if (surfvol[phase_pos_[Vapour]] == 0.0) { if (surfvol[phase_pos_[Vapour]] == 0.0) {
return 0.0; return 0.0;
} }
double Rval = linearInterpolationExtrap(saturated_oil_table_[0], double Rval = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], press); saturated_oil_table_[3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { // Saturated case if (Rval < maxR ) { // Saturated case
return Rval; return Rval;
} else { } else {
return maxR; // Undersaturated case return maxR; // Undersaturated case
} }
} }
void SinglePvtLiveOil::evalRDeriv(const double press, const double* surfvol, void SinglePvtLiveOil::evalRDeriv(const double press, const double* surfvol,
@ -272,19 +272,19 @@ namespace Opm
dRdpval = 0.0; dRdpval = 0.0;
return; return;
} }
Rval = linearInterpolationExtrap(saturated_oil_table_[0], Rval = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], press); saturated_oil_table_[3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { if (Rval < maxR ) {
// Saturated case // Saturated case
dRdpval = linearInterpolDerivative(saturated_oil_table_[0], dRdpval = linearInterpolDerivative(saturated_oil_table_[0],
saturated_oil_table_[3], saturated_oil_table_[3],
press); press);
} else { } else {
// Undersaturated case // Undersaturated case
Rval = maxR; Rval = maxR;
dRdpval = 0.0; dRdpval = 0.0;
} }
} }
@ -293,57 +293,57 @@ namespace Opm
const int item, const int item,
const bool deriv) const const bool deriv) const
{ {
int section; int section;
double Rval = linearInterpolationExtrap(saturated_oil_table_[0], double Rval = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], saturated_oil_table_[3],
press, section); press, section);
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (deriv) { if (deriv) {
if (Rval < maxR ) { // Saturated case if (Rval < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_oil_table_[0], return linearInterpolDerivative(saturated_oil_table_[0],
saturated_oil_table_[item], saturated_oil_table_[item],
press); press);
} else { // Undersaturated case } else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
ASSERT(undersat_oil_tables_[is][0].size() >= 2);
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], maxR); int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) / double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]); (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
ASSERT(undersat_oil_tables_[is][0].size() >= 2); ASSERT(undersat_oil_tables_[is][0].size() >= 2);
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2); ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 = double val1 =
linearInterpolationExtrap(undersat_oil_tables_[is][0], linearInterpolDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item], undersat_oil_tables_[is][item],
press); press);
double val2 = double val2 =
linearInterpolationExtrap(undersat_oil_tables_[is+1][0], linearInterpolDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item], undersat_oil_tables_[is+1][item],
press); press);
double val = val1 + w*(val2 - val1); double val = val1 + w*(val2 - val1);
return val; return val;
} }
} } else {
if (Rval < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
ASSERT(undersat_oil_tables_[is][0].size() >= 2);
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationExtrap(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationExtrap(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
} }
} // namespace Opm } // namespace Opm

18
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp
View File

@ -34,10 +34,10 @@ namespace Opm
class SinglePvtLiveOil : public SinglePvtInterface class SinglePvtLiveOil : public SinglePvtInterface
{ {
public: public:
typedef std::vector<std::vector<std::vector<double> > > table_t; typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtLiveOil(const table_t& pvto); SinglePvtLiveOil(const table_t& pvto);
virtual ~SinglePvtLiveOil(); virtual ~SinglePvtLiveOil();
/// Viscosity as a function of p and z. /// Viscosity as a function of p and z.
virtual void mu(const int n, virtual void mu(const int n,
@ -77,15 +77,15 @@ namespace Opm
double evalR(double press, const double* surfvol) const; double evalR(double press, const double* surfvol) const;
void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const; void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
// item: 1=>1/B 2=>mu; // item: 1=>1/B 2=>mu;
double miscible_oil(const double press, double miscible_oil(const double press,
const double* surfvol, const double* surfvol,
const int item, const int item,
const bool deriv = false) const; const bool deriv = false) const;
// PVT properties of live oil (with dissolved gas) // PVT properties of live oil (with dissolved gas)
std::vector<std::vector<double> > saturated_oil_table_; std::vector<std::vector<double> > saturated_oil_table_;
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_; std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
}; };
} }

244
opm/core/utility/NonuniformTableLinear.hpp Normal file
View File

@ -0,0 +1,244 @@
/*
Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2011, 2012 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
#define OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
#include <cmath>
#include <exception>
#include <vector>
#include <utility>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
namespace Opm
{
/// @brief Exception used for domain errors.
struct OutsideDomainException : public std::exception {};
/// @brief This class uses linear interpolation to compute the value
/// (and its derivative) of a function f sampled at possibly
/// nonuniform points.
/// @tparam T the range type of the function (should be an algebraic ring type)
template<typename T>
class NonuniformTableLinear
{
public:
/// @brief Default constructor.
NonuniformTableLinear();
/// @brief Construct from vectors of x and y values.
/// @param x_values vector of domain values
/// @param y_values vector of corresponding range values.
NonuniformTableLinear(const std::vector<double>& x_values,
const std::vector<T>& y_values);
/// @brief Get the domain.
/// @return the domain as a pair of doubles.
std::pair<double, double> domain();
/// @brief Rescale the domain.
/// @param new_domain the new domain as a pair of doubles.
void rescaleDomain(std::pair<double, double> new_domain);
/// @brief Evaluate the value at x.
/// @param x a domain value
/// @return f(x)
double operator()(const double x) const;
/// @brief Evaluate the derivative at x.
/// @param x a domain value
/// @return f'(x)
double derivative(const double x) const;
/// @brief Evaluate the inverse at y. Requires T to be a double.
/// @param y a range value
/// @return f^{-1}(y)
double inverse(const double y) const;
/// @brief Equality operator.
/// @param other another NonuniformTableLinear.
/// @return true if they are represented exactly alike.
bool operator==(const NonuniformTableLinear& other) const;
/// @brief Policies for how to behave when trying to evaluate outside the domain.
enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2};
/// @brief Sets the behavioural policy for evaluation to the left of the domain.
/// @param rp the policy
void setLeftPolicy(RangePolicy rp);
/// @brief Sets the behavioural policy for evaluation to the right of the domain.
/// @param rp the policy
void setRightPolicy(RangePolicy rp);
protected:
std::vector<double> x_values_;
std::vector<T> y_values_;
mutable std::vector<T> x_values_reversed_;
mutable std::vector<T> y_values_reversed_;
RangePolicy left_;
RangePolicy right_;
};
// A utility function
/// @brief Detect if a sequence is nondecreasing.
/// @tparam FI a forward iterator whose value type has operator< defined.
/// @param beg start of sequence
/// @param end one-beyond-end of sequence
/// @return false if there exists two consecutive values (v1, v2) in the sequence
/// for which v2 < v1, else returns true.
template <typename FI>
bool isNondecreasing(const FI beg, const FI end)
{
if (beg == end) return true;
FI it = beg;
++it;
FI prev = beg;
for (; it != end; ++it, ++prev) {
if (*it < *prev) {
return false;
}
}
return true;
}
// Member implementations.
template<typename T>
inline
NonuniformTableLinear<T>
::NonuniformTableLinear()
: left_(ClosestValue), right_(ClosestValue)
{
}
template<typename T>
inline
NonuniformTableLinear<T>
::NonuniformTableLinear(const std::vector<double>& x_values,
const std::vector<T>& y_values)
: x_values_(x_values), y_values_(y_values),
left_(ClosestValue), right_(ClosestValue)
{
ASSERT(isNondecreasing(x_values.begin(), x_values.end()));
}
template<typename T>
inline std::pair<double, double>
NonuniformTableLinear<T>
::domain()
{
return std::make_pair(x_values_[0], x_values_.back());
}
template<typename T>
inline void
NonuniformTableLinear<T>
::rescaleDomain(std::pair<double, double> new_domain)
{
const double a = x_values_[0];
const double b = x_values_.back();
const double c = new_domain.first;
const double d = new_domain.second;
// x in [a, b] -> x in [c, d]
for (int i = 0; i < int(x_values_.size()); ++i) {
x_values_[i] = (x_values_[i] - a)*(d - c)/(b - a) + c;
}
}
template<typename T>
inline double
NonuniformTableLinear<T>
::operator()(const double x) const
{
return Opm::linearInterpolation(x_values_, y_values_, x);
}
template<typename T>
inline double
NonuniformTableLinear<T>
::derivative(const double x) const
{
return Opm::linearInterpolationDerivative(x_values_, y_values_, x);
}
template<typename T>
inline double
NonuniformTableLinear<T>
::inverse(const double y) const
{
if (y_values_.front() < y_values_.back()) {
return Opm::linearInterpolation(y_values_, x_values_, y);
} else {
if (y_values_reversed_.empty()) {
y_values_reversed_ = y_values_;
std::reverse(y_values_reversed_.begin(), y_values_reversed_.end());
ASSERT(isNondecreasing(y_values_reversed_.begin(), y_values_reversed_.end()));
x_values_reversed_ = x_values_;
std::reverse(x_values_reversed_.begin(), x_values_reversed_.end());
}
return Opm::linearInterpolation(y_values_reversed_, x_values_reversed_, y);
}
}
template<typename T>
inline bool
NonuniformTableLinear<T>
::operator==(const NonuniformTableLinear<T>& other) const
{
return x_values_ == other.x_values_
&& y_values_ == other.y_values_
&& left_ == other.left_
&& right_ == other.right_;
}
template<typename T>
inline void
NonuniformTableLinear<T>
::setLeftPolicy(RangePolicy rp)
{
if (rp != ClosestValue) {
THROW("Only ClosestValue RangePolicy implemented.");
}
left_ = rp;
}
template<typename T>
inline void
NonuniformTableLinear<T>
::setRightPolicy(RangePolicy rp)
{
if (rp != ClosestValue) {
THROW("Only ClosestValue RangePolicy implemented.");
}
right_ = rp;
}
} // namespace Opm
#endif // OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED

2
tests/bo_resprop_test.cpp
View File

@ -43,7 +43,7 @@ int main(int argc, char** argv)
UnstructuredGrid grid; UnstructuredGrid grid;
grid.number_of_cells = 1; grid.number_of_cells = 1;
grid.global_cell = NULL; grid.global_cell = NULL;
Opm::BlackoilPropertiesFromDeck props(deck, grid); Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
const int n = 1; const int n = 1;
double p[n] = { 150e5 }; double p[n] = { 150e5 };