diff --git a/Makefile.am b/Makefile.am index a477f733..2eeea276 100644 --- a/Makefile.am +++ b/Makefile.am @@ -57,6 +57,7 @@ opm/core/fluid/SatFuncStone2.cpp \ opm/core/fluid/SatFuncSimple.cpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \ opm/core/fluid/blackoil/SinglePvtDead.cpp \ +opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp \ opm/core/fluid/blackoil/SinglePvtInterface.cpp \ opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \ opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \ @@ -150,15 +151,18 @@ opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \ opm/core/fluid/RockBasic.hpp \ opm/core/fluid/RockCompressibility.hpp \ opm/core/fluid/RockFromDeck.hpp \ -opm/core/fluid/SaturationPropsBasic.hpp \ -opm/core/fluid/SaturationPropsFromDeck.hpp \ opm/core/fluid/SatFuncStone2.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/blackoil/BlackoilPhases.hpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \ opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \ opm/core/fluid/blackoil/SinglePvtDead.hpp \ +opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp \ opm/core/fluid/blackoil/SinglePvtInterface.hpp \ opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \ opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \ @@ -233,6 +237,7 @@ opm/core/utility/ColumnExtract.hpp \ opm/core/utility/ErrorMacros.hpp \ opm/core/utility/Factory.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \ +opm/core/utility/NonuniformTableLinear.hpp \ opm/core/utility/RootFinders.hpp \ opm/core/utility/SparseTable.hpp \ opm/core/utility/SparseVector.hpp \ diff --git a/examples/sim_2p_comp_reorder.cpp b/examples/sim_2p_comp_reorder.cpp index 9f2b516b..cdc4a5a7 100644 --- a/examples/sim_2p_comp_reorder.cpp +++ b/examples/sim_2p_comp_reorder.cpp @@ -94,7 +94,7 @@ main(int argc, char** argv) // Grid init grid.reset(new GridManager(*deck)); // 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); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Rock compressibility. diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index f7e0172f..73374f00 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -175,7 +175,7 @@ main(int argc, char** argv) // Grid init grid.reset(new Opm::GridManager(deck)); // 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.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); check_well_controls = param.getDefault("check_well_controls", false); diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp index be8211c1..2496d98e 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.cpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.cpp @@ -26,20 +26,20 @@ namespace Opm { BlackoilPropertiesBasic::BlackoilPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells) + const int dim, + const int num_cells) { - double poro = param.getDefault("porosity", 1.0); - using namespace Opm::unit; - using namespace Opm::prefix; - double perm = param.getDefault("permeability", 100.0)*milli*darcy; + double poro = param.getDefault("porosity", 1.0); + using namespace Opm::unit; + using namespace Opm::prefix; + double perm = param.getDefault("permeability", 100.0)*milli*darcy; rock_.init(dim, num_cells, poro, perm); - pvt_.init(param); + pvt_.init(param); satprops_.init(param); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } } BlackoilPropertiesBasic::~BlackoilPropertiesBasic() @@ -90,11 +90,11 @@ namespace Opm /// \param[out] dmudp If non-null: array of nP viscosity derivative values, /// array must be valid before calling. void BlackoilPropertiesBasic::viscosity(const int n, - const double* p, - const double* z, - const int* /*cells*/, - double* mu, - double* dmudp) const + const double* p, + const double* z, + const int* /*cells*/, + double* mu, + double* dmudp) const { if (dmudp) { 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 /// in Fortran order. void BlackoilPropertiesBasic::matrix(const int n, - const double* /*p*/, - const double* /*z*/, - const int* /*cells*/, - double* A, - double* dAdp) const + const double* /*p*/, + const double* /*z*/, + const int* /*cells*/, + double* A, + double* dAdp) const { - const int np = numPhases(); - ASSERT(np <= 2); - double B[2]; // Must be enough since component classes do not handle more than 2. - pvt_.B(1, 0, 0, B); + const int np = numPhases(); + ASSERT(np <= 2); + double B[2]; // Must be enough since component classes do not handle more than 2. + pvt_.B(1, 0, 0, B); // Compute A matrix // #pragma omp parallel for for (int i = 0; i < n; ++i) { @@ -152,8 +152,8 @@ namespace Opm /// of a call to the method matrix(). /// \param[out] rho Array of nP density values, array must be valid before calling. void BlackoilPropertiesBasic::density(const int n, - const double* A, - double* rho) const + const double* A, + double* rho) const { const int np = numPhases(); const double* sdens = pvt_.surfaceDensities(); @@ -186,10 +186,10 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void BlackoilPropertiesBasic::relperm(const int n, - const double* s, - const int* /*cells*/, - double* kr, - double* dkrds) const + const double* s, + const int* /*cells*/, + double* kr, + double* dkrds) const { satprops_.relperm(n, s, kr, dkrds); } @@ -205,10 +205,10 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void BlackoilPropertiesBasic::capPress(const int n, - const double* s, - const int* /*cells*/, - double* pc, - double* dpcds) const + const double* s, + const int* /*cells*/, + double* pc, + double* dpcds) const { satprops_.capPress(n, s, pc, dpcds); } @@ -226,7 +226,7 @@ namespace Opm double* smin, double* smax) const { - satprops_.satRange(n, smin, smax); + satprops_.satRange(n, smin, smax); } diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp index d879036e..44a76013 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.hpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -35,16 +35,16 @@ namespace Opm { public: /// Construct from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1 or 2. - /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". - /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 - /// mu1 [mu2, mu3] (1.0) Viscosity in cP - /// porosity (1.0) Porosity - /// permeability (100.0) Permeability in mD + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 + /// mu1 [mu2, mu3] (1.0) Viscosity in cP + /// porosity (1.0) Porosity + /// permeability (100.0) Permeability in mD BlackoilPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells); + const int dim, + const int num_cells); /// Destructor. virtual ~BlackoilPropertiesBasic(); @@ -151,7 +151,7 @@ namespace Opm 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 the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 5c834fe8..e2457f39 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -19,6 +19,7 @@ #include #include + namespace Opm { @@ -26,12 +27,16 @@ namespace Opm const UnstructuredGrid& grid) { rock_.init(deck, grid); - pvt_.init(deck); - satprops_.init(deck, grid); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } + pvt_.init(deck, 200); + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + 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, @@ -39,13 +44,42 @@ namespace Opm const parameter::ParameterGroup& param) { rock_.init(deck, grid); - const int samples = param.getDefault("dead_tab_size", 1025); - pvt_.init(deck, samples); - satprops_.init(deck, grid, param); + const int pvt_samples = param.getDefault("pvt_tab_size", 200); + pvt_.init(deck, pvt_samples); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + // Unfortunate lack of pointer smartness here... + const int sat_samples = param.getDefault("sat_tab_size", 200); + std::string threephase_model = param.getDefault("threephase_model", "simple"); + bool use_stone2 = (threephase_model == "stone2"); + if (sat_samples > 1) { + if (use_stone2) { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } + } else { + if (use_stone2) { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + 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* 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* 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* smax) const { - satprops_.satRange(n, cells, smin, smax); + satprops_->satRange(n, cells, smin, smax); } diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index 13f5799d..eb8b947b 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -27,6 +27,7 @@ #include #include #include +#include struct UnstructuredGrid; @@ -40,20 +41,23 @@ namespace Opm public: /// Initialize from deck and grid. /// \param[in] deck Deck input parser - /// \param[in] 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) /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid); /// Initialize from deck, grid and parameters. /// \param[in] deck Deck input parser - /// \param[in] 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) /// to logical cartesian indices consistent with the deck. /// \param[in] param Parameters. Accepted parameters include: - /// dead_tab_size (1025) number of uniform sample points for dead-oil pvt tables. - /// tab_size_kr (200) number of uniform sample points for saturation tables. + /// pvt_tab_size (200) number of uniform sample points for dead-oil pvt 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, const UnstructuredGrid& grid, const parameter::ParameterGroup& param); @@ -163,9 +167,9 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \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. @@ -178,7 +182,7 @@ namespace Opm private: RockFromDeck rock_; BlackoilPvtProperties pvt_; - SaturationPropsFromDeck satprops_; + boost::scoped_ptr satprops_; mutable std::vector B_; mutable std::vector dB_; mutable std::vector R_; diff --git a/opm/core/fluid/BlackoilPropertiesInterface.hpp b/opm/core/fluid/BlackoilPropertiesInterface.hpp index dd8a857f..501c8a40 100644 --- a/opm/core/fluid/BlackoilPropertiesInterface.hpp +++ b/opm/core/fluid/BlackoilPropertiesInterface.hpp @@ -138,9 +138,9 @@ namespace Opm double* dpcds) const = 0; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \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. diff --git a/opm/core/fluid/IncompPropertiesBasic.cpp b/opm/core/fluid/IncompPropertiesBasic.cpp index 0ce3c88a..ab68017c 100644 --- a/opm/core/fluid/IncompPropertiesBasic.cpp +++ b/opm/core/fluid/IncompPropertiesBasic.cpp @@ -28,22 +28,22 @@ namespace Opm { IncompPropertiesBasic::IncompPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells) + const int dim, + const int num_cells) { - double poro = param.getDefault("porosity", 1.0); - using namespace Opm::unit; - using namespace Opm::prefix; - double perm = param.getDefault("permeability", 100.0)*milli*darcy; + double poro = param.getDefault("porosity", 1.0); + using namespace Opm::unit; + using namespace Opm::prefix; + double perm = param.getDefault("permeability", 100.0)*milli*darcy; rock_.init(dim, num_cells, poro, perm); - pvt_.init(param); + pvt_.init(param); satprops_.init(param); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } - viscosity_.resize(pvt_.numPhases()); - pvt_.mu(1, 0, 0, &viscosity_[0]); + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } + viscosity_.resize(pvt_.numPhases()); + pvt_.mu(1, 0, 0, &viscosity_[0]); } IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases, @@ -56,14 +56,14 @@ namespace Opm const int num_cells) { rock_.init(dim, num_cells, por, perm); - pvt_.init(num_phases, rho, mu); + pvt_.init(num_phases, rho, mu); satprops_.init(num_phases, relpermfunc); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } - viscosity_.resize(pvt_.numPhases()); - pvt_.mu(1, 0, 0, &viscosity_[0]); + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } + viscosity_.resize(pvt_.numPhases()); + pvt_.mu(1, 0, 0, &viscosity_[0]); } IncompPropertiesBasic::~IncompPropertiesBasic() @@ -109,7 +109,7 @@ namespace Opm /// \return Array of P viscosity values. const double* IncompPropertiesBasic::viscosity() const { - return &viscosity_[0]; + return &viscosity_[0]; } /// \return Array of P density values. @@ -117,7 +117,7 @@ namespace Opm { // No difference between reservoir and surface densities // modelled by this class. - return pvt_.surfaceDensities(); + return pvt_.surfaceDensities(); } /// \return Array of P density values. @@ -125,7 +125,7 @@ namespace Opm { // No difference between reservoir and surface densities // modelled by this class. - return pvt_.surfaceDensities(); + return pvt_.surfaceDensities(); } /// \param[in] n Number of data points. @@ -138,10 +138,10 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesBasic::relperm(const int n, - const double* s, - const int* /*cells*/, - double* kr, - double* dkrds) const + const double* s, + const int* /*cells*/, + double* kr, + double* dkrds) const { satprops_.relperm(n, s, kr, dkrds); } @@ -157,10 +157,10 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesBasic::capPress(const int n, - const double* s, - const int* /*cells*/, - double* pc, - double* dpcds) const + const double* s, + const int* /*cells*/, + double* pc, + double* dpcds) const { 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] smax Array of nP maximum s values, array must be valid before calling. void IncompPropertiesBasic::satRange(const int n, - const int* /*cells*/, - double* smin, - double* smax) const + const int* /*cells*/, + double* smin, + double* smax) const { - satprops_.satRange(n, smin, smax); + satprops_.satRange(n, smin, smax); } } // namespace Opm diff --git a/opm/core/fluid/IncompPropertiesBasic.hpp b/opm/core/fluid/IncompPropertiesBasic.hpp index b90b0876..72c39b1b 100644 --- a/opm/core/fluid/IncompPropertiesBasic.hpp +++ b/opm/core/fluid/IncompPropertiesBasic.hpp @@ -42,29 +42,29 @@ namespace Opm { public: /// Construct from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1 or 2. - /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". - /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 - /// mu1 [mu2, mu3] (1.0) Viscosity in cP - /// porosity (1.0) Porosity - /// permeability (100.0) Permeability in mD + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 + /// mu1 [mu2, mu3] (1.0) Viscosity in cP + /// porosity (1.0) Porosity + /// permeability (100.0) Permeability in mD IncompPropertiesBasic(const parameter::ParameterGroup& param, - const int dim, - const int num_cells); + const int dim, + const int num_cells); /// Construct from arguments a basic two phase fluid. IncompPropertiesBasic(const int num_phases, const SaturationPropsBasic::RelPermFunc& relpermfunc, const std::vector& rho, - const std::vector& mu, + const std::vector& mu, const double porosity, const double permeability, const int dim, - const int num_cells); + const int num_cells); - /// Destructor. + /// Destructor. virtual ~IncompPropertiesBasic(); // ---- Rock interface ---- @@ -132,9 +132,9 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \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. @@ -145,9 +145,9 @@ namespace Opm double* smax) const; private: RockBasic rock_; - PvtPropertiesBasic pvt_; + PvtPropertiesBasic pvt_; SaturationPropsBasic satprops_; - std::vector viscosity_; + std::vector viscosity_; }; diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp index 570154bf..9aa69098 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.cpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp @@ -27,15 +27,15 @@ namespace Opm { IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid) + const UnstructuredGrid& grid) { rock_.init(deck, grid); - pvt_.init(deck); - satprops_.init(deck, grid); - if (pvt_.numPhases() != satprops_.numPhases()) { - THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); - } + pvt_.init(deck); + satprops_.init(deck, grid, 200); + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } } IncompPropertiesFromDeck::~IncompPropertiesFromDeck() @@ -81,19 +81,19 @@ namespace Opm /// \return Array of P viscosity values. const double* IncompPropertiesFromDeck::viscosity() const { - return pvt_.viscosity(); + return pvt_.viscosity(); } /// \return Array of P density values. const double* IncompPropertiesFromDeck::density() const { - return pvt_.reservoirDensities(); + return pvt_.reservoirDensities(); } /// \return Array of P density values. const double* IncompPropertiesFromDeck::surfaceDensity() const { - return pvt_.surfaceDensities(); + return pvt_.surfaceDensities(); } /// \param[in] n Number of data points. @@ -106,10 +106,10 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesFromDeck::relperm(const int n, - const double* s, - const int* cells, - double* kr, - double* dkrds) const + const double* s, + const int* cells, + double* kr, + double* dkrds) const { satprops_.relperm(n, s, cells, kr, dkrds); } @@ -125,10 +125,10 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m_01 ...) void IncompPropertiesFromDeck::capPress(const int n, - const double* s, - const int* cells, - double* pc, - double* dpcds) const + const double* s, + const int* cells, + double* pc, + double* dpcds) const { 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] smax Array of nP maximum s values, array must be valid before calling. void IncompPropertiesFromDeck::satRange(const int n, - const int* cells, - double* smin, - double* smax) const + const int* cells, + double* smin, + double* smax) const { - satprops_.satRange(n, cells, smin, smax); + satprops_.satRange(n, cells, smin, smax); } } // namespace Opm diff --git a/opm/core/fluid/IncompPropertiesFromDeck.hpp b/opm/core/fluid/IncompPropertiesFromDeck.hpp index 68623e5c..6680eaa4 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.hpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.hpp @@ -47,13 +47,13 @@ namespace Opm public: /// Initialize from deck and grid. /// \param deck Deck input parser - /// \param grid Grid to which property object applies, needed for the + /// \param grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. IncompPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid); - /// Destructor. + /// Destructor. virtual ~IncompPropertiesFromDeck(); // ---- Rock interface ---- @@ -121,9 +121,9 @@ namespace Opm double* dpcds) const; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \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. @@ -134,8 +134,8 @@ namespace Opm double* smax) const; private: RockFromDeck rock_; - PvtPropertiesIncompFromDeck pvt_; - SaturationPropsFromDeck satprops_; + PvtPropertiesIncompFromDeck pvt_; + SaturationPropsFromDeck satprops_; }; diff --git a/opm/core/fluid/IncompPropertiesInterface.hpp b/opm/core/fluid/IncompPropertiesInterface.hpp index a5025e8a..4f8bfb51 100644 --- a/opm/core/fluid/IncompPropertiesInterface.hpp +++ b/opm/core/fluid/IncompPropertiesInterface.hpp @@ -109,9 +109,9 @@ namespace Opm double* pc, double* dpcds) const = 0; - /// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \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. diff --git a/opm/core/fluid/PvtPropertiesBasic.cpp b/opm/core/fluid/PvtPropertiesBasic.cpp index 3d7d7bb2..5220502a 100644 --- a/opm/core/fluid/PvtPropertiesBasic.cpp +++ b/opm/core/fluid/PvtPropertiesBasic.cpp @@ -34,41 +34,41 @@ namespace Opm void PvtPropertiesBasic::init(const parameter::ParameterGroup& param) { - int num_phases = param.getDefault("num_phases", 2); - if (num_phases > 3 || num_phases < 1) { - THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); - } - density_.resize(num_phases); - viscosity_.resize(num_phases); - // We currently do not allow the user to set B. - formation_volume_factor_.clear(); - formation_volume_factor_.resize(num_phases, 1.0); + int num_phases = param.getDefault("num_phases", 2); + if (num_phases > 3 || num_phases < 1) { + THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); + } + density_.resize(num_phases); + viscosity_.resize(num_phases); + // We currently do not allow the user to set B. + formation_volume_factor_.clear(); + formation_volume_factor_.resize(num_phases, 1.0); - // Setting mu and rho from parameters - using namespace Opm::prefix; - using namespace Opm::unit; - const double kgpm3 = kilogram/cubic(meter); - const double cP = centi*Poise; - std::string rname[3] = { "rho1", "rho2", "rho3" }; - double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 }; - std::string vname[3] = { "mu1", "mu2", "mu3" }; - double vdefault[3] = { 1.0, 1.0, 1.0 }; - for (int phase = 0; phase < num_phases; ++phase) { - density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]); - viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]); - } + // Setting mu and rho from parameters + using namespace Opm::prefix; + using namespace Opm::unit; + const double kgpm3 = kilogram/cubic(meter); + const double cP = centi*Poise; + std::string rname[3] = { "rho1", "rho2", "rho3" }; + double rdefault[3] = { 1.0e3, 1.0e3, 1.0e3 }; + std::string vname[3] = { "mu1", "mu2", "mu3" }; + double vdefault[3] = { 1.0, 1.0, 1.0 }; + for (int phase = 0; phase < num_phases; ++phase) { + density_[phase] = kgpm3*param.getDefault(rname[phase], rdefault[phase]); + viscosity_[phase] = cP*param.getDefault(vname[phase], vdefault[phase]); + } } void PvtPropertiesBasic::init(const int num_phases, const std::vector& rho, const std::vector& visc) { - if (num_phases > 3 || num_phases < 1) { - THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); - } - // We currently do not allow the user to set B. - formation_volume_factor_.clear(); - formation_volume_factor_.resize(num_phases, 1.0); + if (num_phases > 3 || num_phases < 1) { + THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases); + } + // We currently do not allow the user to set B. + formation_volume_factor_.clear(); + formation_volume_factor_.resize(num_phases, 1.0); density_ = rho; viscosity_ = visc; } @@ -87,69 +87,69 @@ namespace Opm void PvtPropertiesBasic::mu(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_mu) const + const double* /*p*/, + const double* /*z*/, + double* output_mu) const { - const int np = numPhases(); + const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[np*i + phase] = viscosity_[phase]; } - } + } } void PvtPropertiesBasic::B(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_B) const + const double* /*p*/, + const double* /*z*/, + double* output_B) const { - const int np = numPhases(); + const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[np*i + phase] = formation_volume_factor_[phase]; } - } + } } void PvtPropertiesBasic::dBdp(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_B, - double* output_dBdp) const + const double* /*p*/, + const double* /*z*/, + double* output_B, + double* output_dBdp) const { - const int np = numPhases(); + const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[np*i + phase] = formation_volume_factor_[phase]; output_dBdp[np*i + phase] = 0.0; } - } + } } void PvtPropertiesBasic::R(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R) const + const double* /*p*/, + const double* /*z*/, + double* output_R) const { - const int np = numPhases(); - std::fill(output_R, output_R + n*np, 0.0); + const int np = numPhases(); + std::fill(output_R, output_R + n*np, 0.0); } void PvtPropertiesBasic::dRdp(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const + const double* /*p*/, + const double* /*z*/, + double* output_R, + double* output_dRdp) const { - const int np = numPhases(); - std::fill(output_R, output_R + n*np, 0.0); - std::fill(output_dRdp, output_dRdp + n*np, 0.0); + const int np = numPhases(); + std::fill(output_R, output_R + n*np, 0.0); + std::fill(output_dRdp, output_dRdp + n*np, 0.0); } } // namespace Opm diff --git a/opm/core/fluid/PvtPropertiesBasic.hpp b/opm/core/fluid/PvtPropertiesBasic.hpp index d1a3e9f3..3e30e807 100644 --- a/opm/core/fluid/PvtPropertiesBasic.hpp +++ b/opm/core/fluid/PvtPropertiesBasic.hpp @@ -38,11 +38,11 @@ namespace Opm PvtPropertiesBasic(); /// Initialize from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1, 2 or 3. - /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 - /// mu1 [mu2, mu3] (1.0) Viscosity in cP - void init(const parameter::ParameterGroup& param); + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1, 2 or 3. + /// rho1 [rho2, rho3] (1.0e3) Density in kg/m^3 + /// mu1 [mu2, mu3] (1.0) Viscosity in cP + void init(const parameter::ParameterGroup& param); /// Initialize from arguments. /// Basic multi phase fluid pvt properties. @@ -55,7 +55,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities() const; /// Viscosity as a function of p and z. void mu(const int n, @@ -90,9 +90,9 @@ namespace Opm double* output_dRdp) const; private: - std::vector density_; - std::vector viscosity_; - std::vector formation_volume_factor_; + std::vector density_; + std::vector viscosity_; + std::vector formation_volume_factor_; }; } diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp index 62cc1b23..34a4ba00 100644 --- a/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp +++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.cpp @@ -38,54 +38,54 @@ namespace Opm { typedef std::vector > > table_t; // 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); - if (phase_usage.phase_used[PhaseUsage::Vapour] || - !phase_usage.phase_used[PhaseUsage::Aqua] || - !phase_usage.phase_used[PhaseUsage::Liquid]) { - THROW("PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n"); - } + if (phase_usage.phase_used[PhaseUsage::Vapour] || + !phase_usage.phase_used[PhaseUsage::Aqua] || + !phase_usage.phase_used[PhaseUsage::Liquid]) { + 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. - if (deck.hasField("DENSITY")) { - const std::vector& d = deck.getDENSITY().densities_[region_number]; - 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::Liquid]] = d[ECL_oil]; - } else { - THROW("Input is missing DENSITY\n"); - } + // Surface densities. Accounting for different orders in eclipse and our code. + if (deck.hasField("DENSITY")) { + const std::vector& d = deck.getDENSITY().densities_[region_number]; + 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::Liquid]] = d[ECL_oil]; + } else { + THROW("Input is missing DENSITY\n"); + } // Make reservoir densities the same as surface densities initially. // We will modify them with formation volume factors if found. reservoir_density_ = surface_density_; // Water viscosity. - if (deck.hasField("PVTW")) { - const std::vector& pvtw = deck.getPVTW().pvtw_[region_number]; - if (pvtw[2] != 0.0 || pvtw[4] != 0.0) { - MESSAGE("Compressibility effects in PVTW are ignored."); - } + if (deck.hasField("PVTW")) { + const std::vector& pvtw = deck.getPVTW().pvtw_[region_number]; + if (pvtw[2] != 0.0 || pvtw[4] != 0.0) { + MESSAGE("Compressibility effects in PVTW are ignored."); + } reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1]; - viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3]; - } else { - // Eclipse 100 default. - // viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise; - THROW("Input is missing PVTW\n"); - } + viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3]; + } else { + // Eclipse 100 default. + // viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise; + THROW("Input is missing PVTW\n"); + } // Oil viscosity. - if (deck.hasField("PVCDO")) { - const std::vector& pvcdo = deck.getPVCDO().pvcdo_[region_number]; - if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) { - MESSAGE("Compressibility effects in PVCDO are ignored."); - } + if (deck.hasField("PVCDO")) { + const std::vector& pvcdo = deck.getPVCDO().pvcdo_[region_number]; + if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) { + MESSAGE("Compressibility effects in PVCDO are ignored."); + } reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1]; - viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3]; - } else { - THROW("Input is missing PVCDO\n"); - } + viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3]; + } else { + THROW("Input is missing PVCDO\n"); + } } const double* PvtPropertiesIncompFromDeck::surfaceDensities() const diff --git a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp index acbabaa4..01bbeed1 100644 --- a/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp +++ b/opm/core/fluid/PvtPropertiesIncompFromDeck.hpp @@ -39,14 +39,14 @@ namespace Opm PvtPropertiesIncompFromDeck(); /// Initialize from deck. - void init(const EclipseGridParser& deck); + void init(const EclipseGridParser& deck); /// Number of active phases. int numPhases() const; /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities() const; /// Densities of stock components at reservoir conditions. /// 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, /// we need to handle this difference. /// \return Array of size numPhases(). - const double* reservoirDensities() const; + const double* reservoirDensities() const; /// Viscosities. const double* viscosity() const; private: - std::tr1::array surface_density_; - std::tr1::array reservoir_density_; - std::tr1::array viscosity_; + std::tr1::array surface_density_; + std::tr1::array reservoir_density_; + std::tr1::array viscosity_; }; } diff --git a/opm/core/fluid/RockBasic.cpp b/opm/core/fluid/RockBasic.cpp index a1423b76..b66cccce 100644 --- a/opm/core/fluid/RockBasic.cpp +++ b/opm/core/fluid/RockBasic.cpp @@ -25,29 +25,29 @@ namespace Opm /// Default constructor. RockBasic::RockBasic() - : dimensions_(-1) + : dimensions_(-1) { } /// Initialize with homogenous porosity and permeability. void RockBasic::init(const int dimensions, - const int num_cells, - const double poro, - const double perm) + const int num_cells, + const double poro, + const double perm) { - dimensions_ = dimensions; - porosity_.clear(); - porosity_.resize(num_cells, poro); - permeability_.clear(); - const int dsq = dimensions*dimensions; - permeability_.resize(num_cells*dsq, 0.0); + dimensions_ = dimensions; + porosity_.clear(); + porosity_.resize(num_cells, poro); + permeability_.clear(); + const int dsq = dimensions*dimensions; + permeability_.resize(num_cells*dsq, 0.0); // #pragma omp parallel for - for (int i = 0; i < num_cells; ++i) { - for (int d = 0; d < dimensions; ++d) { - permeability_[dsq*i + dimensions*d + d] = perm; - } - } + for (int i = 0; i < num_cells; ++i) { + for (int d = 0; d < dimensions; ++d) { + permeability_[dsq*i + dimensions*d + d] = perm; + } + } } diff --git a/opm/core/fluid/RockBasic.hpp b/opm/core/fluid/RockBasic.hpp index 96563780..2adc3ffd 100644 --- a/opm/core/fluid/RockBasic.hpp +++ b/opm/core/fluid/RockBasic.hpp @@ -35,9 +35,9 @@ namespace Opm /// Initialize with homogenous porosity and permeability. void init(const int dimensions, - const int num_cells, - const double poro, - const double perm); + const int num_cells, + const double poro, + const double perm); /// \return D, the number of spatial dimensions. int numDimensions() const @@ -66,7 +66,7 @@ namespace Opm } private: - int dimensions_; + int dimensions_; std::vector porosity_; std::vector permeability_; }; diff --git a/opm/core/fluid/RockCompressibility.cpp b/opm/core/fluid/RockCompressibility.cpp index 10486101..af65d852 100644 --- a/opm/core/fluid/RockCompressibility.cpp +++ b/opm/core/fluid/RockCompressibility.cpp @@ -69,8 +69,8 @@ namespace Opm const double cpnorm = rock_comp_*(pressure - pref_); return (1.0 + cpnorm + 0.5*cpnorm*cpnorm); } else { - // return Opm::linearInterpolation(p_, poromult_, pressure); - return Opm::linearInterpolationExtrap(p_, poromult_, pressure); + // return Opm::linearInterpolation(p_, poromult_, pressure); + return Opm::linearInterpolationExtrap(p_, poromult_, pressure); } } @@ -81,7 +81,7 @@ namespace Opm } else { //const double poromult = Opm::linearInterpolation(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); return dporomultdp/poromult; diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 7ba8903f..ec7db9a6 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -51,7 +51,7 @@ namespace Opm /// Initialize from deck and cell mapping. /// \param deck Deck input parser - /// \param grid grid to which property object applies, needed for the + /// \param grid grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. void RockFromDeck::init(const EclipseGridParser& deck, diff --git a/opm/core/fluid/RockFromDeck.hpp b/opm/core/fluid/RockFromDeck.hpp index 65027e42..63e71a6e 100644 --- a/opm/core/fluid/RockFromDeck.hpp +++ b/opm/core/fluid/RockFromDeck.hpp @@ -37,7 +37,7 @@ namespace Opm /// Initialize from deck and grid. /// \param deck Deck input parser - /// \param grid Grid to which property object applies, needed for the + /// \param grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. void init(const EclipseGridParser& deck, diff --git a/opm/core/fluid/SatFuncSimple.cpp b/opm/core/fluid/SatFuncSimple.cpp index 740af43c..7e226a4a 100644 --- a/opm/core/fluid/SatFuncSimple.cpp +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -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 . +*/ + #include #include #include @@ -6,21 +25,14 @@ #include #include #include + namespace Opm { - void SatFuncSimple::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, + void SatFuncSimpleUniform::init(const EclipseGridParser& deck, const int table_num, const PhaseUsage phase_usg, 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) { // 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; 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; 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& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& pcow = swof_table[table_num][3]; + krw_ = NonuniformTableLinear(sw, krw); + krow_ = NonuniformTableLinear(sw, krow); + pcow_ = NonuniformTableLinear(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& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& pcog = sgof_table[table_num][3]; + krg_ = NonuniformTableLinear(sg, krg); + krog_ = NonuniformTableLinear(sg, krog); + pcog_ = NonuniformTableLinear(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 // is very hard experimentally, usually one extends two-phase diff --git a/opm/core/fluid/SatFuncSimple.hpp b/opm/core/fluid/SatFuncSimple.hpp index a3388848..338a359e 100644 --- a/opm/core/fluid/SatFuncSimple.hpp +++ b/opm/core/fluid/SatFuncSimple.hpp @@ -18,17 +18,21 @@ */ #ifndef SATFUNCSIMPLE_HPP #define SATFUNCSIMPLE_HPP + #include #include +#include #include #include + namespace Opm { - class SatFuncSimple: public BlackoilPhases + class SatFuncSimpleUniform : public BlackoilPhases { public: - void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); - void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg, + 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; @@ -46,5 +50,31 @@ namespace Opm UniformTableLinear pcog_; 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 krw_; + NonuniformTableLinear krow_; + NonuniformTableLinear pcow_; + NonuniformTableLinear krg_; + NonuniformTableLinear krog_; + NonuniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + } // namespace Opm #endif // SATFUNCSIMPLE_HPP diff --git a/opm/core/fluid/SatFuncStone2.cpp b/opm/core/fluid/SatFuncStone2.cpp index a03bfb66..062f8534 100644 --- a/opm/core/fluid/SatFuncStone2.cpp +++ b/opm/core/fluid/SatFuncStone2.cpp @@ -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 . +*/ + #include #include #include @@ -6,20 +25,14 @@ #include #include #include + namespace Opm { - void SatFuncStone2::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, + void SatFuncStone2Uniform::init(const EclipseGridParser& deck, const int table_num, const PhaseUsage phase_usg, 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) { // 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; 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; 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 // 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& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& pcow = swof_table[table_num][3]; + krw_ = NonuniformTableLinear(sw, krw); + krow_ = NonuniformTableLinear(sw, krow); + pcow_ = NonuniformTableLinear(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& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& pcog = sgof_table[table_num][3]; + krg_ = NonuniformTableLinear(sg, krg); + krog_ = NonuniformTableLinear(sg, krog); + pcog_ = NonuniformTableLinear(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 diff --git a/opm/core/fluid/SatFuncStone2.hpp b/opm/core/fluid/SatFuncStone2.hpp index ae4a6d38..532bf7ce 100644 --- a/opm/core/fluid/SatFuncStone2.hpp +++ b/opm/core/fluid/SatFuncStone2.hpp @@ -18,17 +18,21 @@ */ #ifndef SATFUNCSTONE2_HPP #define SATFUNCSTONE2_HPP + #include #include +#include #include #include + namespace Opm { - class SatFuncStone2: public BlackoilPhases + class SatFuncStone2Uniform : public BlackoilPhases { public: - void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); - void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg, + 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; @@ -46,5 +50,31 @@ namespace Opm UniformTableLinear pcog_; 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 krw_; + NonuniformTableLinear krow_; + NonuniformTableLinear pcow_; + NonuniformTableLinear krg_; + NonuniformTableLinear krog_; + NonuniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + } // namespace Opm #endif // SATFUNCSTONE2_HPP diff --git a/opm/core/fluid/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp index 760f934d..10ddf23d 100644 --- a/opm/core/fluid/SaturationPropsBasic.cpp +++ b/opm/core/fluid/SaturationPropsBasic.cpp @@ -29,64 +29,64 @@ namespace Opm namespace { - struct KrFunConstant - { - double kr(double) - { - return 1.0; - } - double dkrds(double) - { - return 0.0; - } - }; + struct KrFunConstant + { + double kr(double) + { + return 1.0; + } + double dkrds(double) + { + return 0.0; + } + }; - struct KrFunLinear - { - double kr(double s) - { - return s; - } - double dkrds(double) - { - return 1.0; - } - }; + struct KrFunLinear + { + double kr(double s) + { + return s; + } + double dkrds(double) + { + return 1.0; + } + }; - struct KrFunQuadratic - { - double kr(double s) - { - return s*s; - } - double dkrds(double s) - { - return 2.0*s; - } - }; + struct KrFunQuadratic + { + double kr(double s) + { + return s*s; + } + double dkrds(double s) + { + return 2.0*s; + } + }; - template - static inline void evalAllKrDeriv(const int n, const int np, - const double* s, double* kr, double* dkrds, Fun fun) - { - if (dkrds == 0) { + template + static inline void evalAllKrDeriv(const int n, const int np, + const double* s, double* kr, double* dkrds, Fun fun) + { + if (dkrds == 0) { // #pragma omp parallel for - for (int i = 0; i < n*np; ++i) { - kr[i] = fun.kr(s[i]); - } - return; - } + for (int i = 0; i < n*np; ++i) { + kr[i] = fun.kr(s[i]); + } + return; + } // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); - for (int phase = 0; phase < np; ++phase) { - kr[i*np + phase] = fun.kr(s[i*np + phase]); - // Only diagonal elements in derivative. - dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]); - } - } - } + for (int i = 0; i < n; ++i) { + std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); + for (int phase = 0; phase < np; ++phase) { + kr[i*np + phase] = fun.kr(s[i*np + phase]); + // Only diagonal elements in derivative. + dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]); + } + } + } } // anon namespace @@ -109,25 +109,25 @@ namespace Opm /// Initialize from parameters. void SaturationPropsBasic::init(const parameter::ParameterGroup& param) { - int num_phases = param.getDefault("num_phases", 2); - if (num_phases > 2 || num_phases < 1) { - THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases); - } + int num_phases = param.getDefault("num_phases", 2); + if (num_phases > 2 || num_phases < 1) { + THROW("SaturationPropsBasic::init() illegal 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("Linear")); - if (rpf == "Constant") { - relperm_func_ = Constant; - if(num_phases!=1){ - THROW("Constant relperm with more than one phase???"); - } - } else if (rpf == "Linear") { - relperm_func_ = Linear; - } else if (rpf == "Quadratic") { - relperm_func_ = Quadratic; - } else { - THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf); - } + //std::string rpf = param.getDefault("relperm_func", std::string("Unset")); + std::string rpf = param.getDefault("relperm_func", std::string("Linear")); + if (rpf == "Constant") { + relperm_func_ = Constant; + if(num_phases!=1){ + THROW("Constant relperm with more than one phase???"); + } + } else if (rpf == "Linear") { + relperm_func_ = Linear; + } else if (rpf == "Quadratic") { + relperm_func_ = Quadratic; + } else { + THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf); + } } @@ -136,7 +136,7 @@ namespace Opm /// \return P, the number of phases. int SaturationPropsBasic::numPhases() const { - return num_phases_; + return num_phases_; } @@ -152,29 +152,29 @@ namespace Opm /// m_{ij} = \frac{dkr_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void SaturationPropsBasic::relperm(const int n, - const double* s, - double* kr, - double* dkrds) const + const double* s, + double* kr, + double* dkrds) const { - switch (relperm_func_) { - case Constant: - { - evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant()); - break; - } - case Linear: - { - evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear()); - break; - } - case Quadratic: - { - evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic()); - break; - } - default: - THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_); - } + switch (relperm_func_) { + case Constant: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant()); + break; + } + case Linear: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear()); + break; + } + case Quadratic: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic()); + break; + } + default: + THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_); + } } @@ -190,13 +190,13 @@ namespace Opm /// m_{ij} = \frac{dpc_i}{ds^j}, /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) void SaturationPropsBasic::capPress(const int n, - const double* /*s*/, - double* pc, - double* dpcds) const + const double* /*s*/, + double* pc, + double* dpcds) const { - std::fill(pc, pc + num_phases_*n, 0.0); + std::fill(pc, pc + num_phases_*n, 0.0); 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] smax Array of nP maximum s values, array must be valid before calling. void SaturationPropsBasic::satRange(const int n, - double* smin, - double* smax) const + double* smin, + double* smax) const { - std::fill(smin, smin + num_phases_*n, 0.0); - std::fill(smax, smax + num_phases_*n, 1.0); + std::fill(smin, smin + num_phases_*n, 0.0); + std::fill(smax, smax + num_phases_*n, 1.0); } diff --git a/opm/core/fluid/SaturationPropsBasic.hpp b/opm/core/fluid/SaturationPropsBasic.hpp index 789e0a35..2f90672d 100644 --- a/opm/core/fluid/SaturationPropsBasic.hpp +++ b/opm/core/fluid/SaturationPropsBasic.hpp @@ -40,16 +40,16 @@ namespace Opm SaturationPropsBasic(); /// Initialize from parameters. - /// The following parameters are accepted (defaults): - /// num_phases (2) Must be 1 or 2. - /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". void init(const parameter::ParameterGroup& param); - enum RelPermFunc { Constant, Linear, Quadratic }; + enum RelPermFunc { Constant, Linear, Quadratic }; /// Initialize from arguments a basic Saturation property. void init(const int num_phases, - const RelPermFunc& relperm_func) + const RelPermFunc& relperm_func) { num_phases_ = num_phases; relperm_func_ = relperm_func; @@ -86,18 +86,18 @@ namespace Opm double* pc, 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[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 satRange(const int n, - double* smin, - double* smax) const; + void satRange(const int n, + double* smin, + double* smax) const; private: - int num_phases_; - RelPermFunc relperm_func_; + int num_phases_; + RelPermFunc relperm_func_; }; diff --git a/opm/core/fluid/SaturationPropsFromDeck.cpp b/opm/core/fluid/SaturationPropsFromDeck.cpp index 3481b874..83c942c5 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/opm/core/fluid/SaturationPropsFromDeck.cpp @@ -19,7 +19,6 @@ #include #include -#include #include #include #include @@ -27,236 +26,9 @@ namespace Opm { - /// Default constructor. - SaturationPropsFromDeck::SaturationPropsFromDeck() - { - } + // This file should be removed in the future. + // 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& 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& 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 diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index bd70d802..0a0599ec 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -19,9 +19,10 @@ #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED + +#include #include #include -#include #include #include #include @@ -32,23 +33,31 @@ struct UnstructuredGrid; 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 SaturationPropsFromDeck : public SaturationPropsInterface { public: /// Default constructor. SaturationPropsFromDeck(); /// Initialize from deck and grid. - /// \param deck Deck input parser - /// \param grid Grid to which property object applies, needed for the + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. - void init(const EclipseGridParser& deck, - const UnstructuredGrid& grid); - + /// \param[in] samples Number of uniform sample points for saturation tables. + /// NOTE: samples will only be used with the SatFuncSetUniform template argument. void init(const EclipseGridParser& deck, const UnstructuredGrid& grid, - const parameter::ParameterGroup& param); + const int samples); /// \return P, the number of phases. int numPhases() const; @@ -83,22 +92,23 @@ namespace Opm double* pc, 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[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 satRange(const int n, + void satRange(const int n, const int* cells, - double* smin, - double* smax) const; + double* smin, + double* smax) const; private: PhaseUsage phase_usage_; - typedef SatFuncSimple satfunc_t; - std::vector satfuncset_; + std::vector satfuncset_; std::vector 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 +#include #endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED diff --git a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp new file mode 100644 index 00000000..f7a80b45 --- /dev/null +++ b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp @@ -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 . +*/ + +#ifndef OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED +#define OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED + + +#include +#include +#include +#include + +namespace Opm +{ + + + // ----------- Methods of SaturationPropsFromDeck --------- + + + /// Default constructor. + template + SaturationPropsFromDeck::SaturationPropsFromDeck() + { + } + + /// Initialize from deck. + template + void SaturationPropsFromDeck::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& 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 + 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 ...) + template + 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 ...) + template + 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. + template + 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. + template + const typename SaturationPropsFromDeck::Funcs& + SaturationPropsFromDeck::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 diff --git a/opm/core/fluid/SaturationPropsInterface.hpp b/opm/core/fluid/SaturationPropsInterface.hpp new file mode 100644 index 00000000..30f277cc --- /dev/null +++ b/opm/core/fluid/SaturationPropsInterface.hpp @@ -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 . +*/ + +#ifndef OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED +#define OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED + +#include + + +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 diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp index 39b4ba05..4a0004ee 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -43,14 +44,14 @@ namespace Opm { typedef std::vector > > table_t; // If we need multiple regions, this class and the SinglePvt* classes must change. - region_number_ = 0; + region_number_ = 0; phase_usage_ = phaseUsageFromDeck(deck); - // Surface densities. Accounting for different orders in eclipse and our code. - if (deck.hasField("DENSITY")) { - const std::vector& d = deck.getDENSITY().densities_[region_number_]; - enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; + // Surface densities. Accounting for different orders in eclipse and our code. + if (deck.hasField("DENSITY")) { + const std::vector& d = deck.getDENSITY().densities_[region_number_]; + enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; if (phase_usage_.phase_used[Aqua]) { densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water]; } @@ -60,9 +61,9 @@ namespace Opm if (phase_usage_.phase_used[Liquid]) { densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil]; } - } else { - THROW("Input is missing DENSITY\n"); - } + } else { + THROW("Input is missing DENSITY\n"); + } // Set the properties. props_.resize(phase_usage_.num_phases); @@ -78,7 +79,11 @@ namespace Opm // Oil PVT if (phase_usage_.phase_used[Liquid]) { 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")) { props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_)); } else if (deck.hasField("PVCDO")) { @@ -87,10 +92,14 @@ namespace Opm THROW("Input is missing PVDO or PVTO\n"); } } - // Gas PVT + // Gas PVT if (phase_usage_.phase_used[Vapour]) { 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")) { props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_)); } else { diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 191c6fa8..8633fd9e 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp @@ -47,7 +47,13 @@ namespace Opm BlackoilPvtProperties(); /// 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. int numPhases() const; @@ -64,7 +70,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities() const; /// Viscosity as a function of p and z. void mu(const int n, @@ -105,11 +111,11 @@ namespace Opm PhaseUsage phase_usage_; - int region_number_; + int region_number_; std::vector > props_; - double densities_[MaxNumPhases]; + double densities_[MaxNumPhases]; mutable std::vector data1_; mutable std::vector data2_; }; diff --git a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp b/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp index 557bc1e8..9226d0d8 100644 --- a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp +++ b/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp @@ -40,12 +40,12 @@ namespace Opm public: typedef std::vector > table_t; - SinglePvtConstCompr(const table_t& pvtw) + SinglePvtConstCompr(const table_t& pvtw) { - const int region_number = 0; - if (pvtw.size() != 1) { - THROW("More than one PVD-region"); - } + const int region_number = 0; + if (pvtw.size() != 1) { + THROW("More than one PVD-region"); + } ref_press_ = pvtw[region_number][0]; ref_B_ = pvtw[region_number][1]; comp_ = pvtw[region_number][2]; @@ -53,7 +53,7 @@ namespace Opm visc_comp_ = pvtw[region_number][4]; } - SinglePvtConstCompr(double visc) + SinglePvtConstCompr(double visc) : ref_press_(0.0), ref_B_(1.0), comp_(0.0), @@ -62,7 +62,7 @@ namespace Opm { } - virtual ~SinglePvtConstCompr() + virtual ~SinglePvtConstCompr() { } diff --git a/opm/core/fluid/blackoil/SinglePvtDead.cpp b/opm/core/fluid/blackoil/SinglePvtDead.cpp index 1a0993be..34aeaf6e 100644 --- a/opm/core/fluid/blackoil/SinglePvtDead.cpp +++ b/opm/core/fluid/blackoil/SinglePvtDead.cpp @@ -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). @@ -17,8 +17,8 @@ along with OPM. If not, see . */ + #include -#include #include // Extra includes for debug dumping of tables. @@ -33,25 +33,25 @@ namespace Opm // Member functions //------------------------------------------------------------------------- /// Constructor - SinglePvtDead::SinglePvtDead(const table_t& pvd_table, const int samples) + SinglePvtDead::SinglePvtDead(const table_t& pvd_table) { - const int region_number = 0; - if (pvd_table.size() != 1) { - THROW("More than one PVT-region"); - } + 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(); + // Copy data + const int sz = pvd_table[region_number][0].size(); std::vector press(sz); std::vector B_inv(sz); std::vector visc(sz); - for (int i = 0; i < sz; ++i) { + 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_); + } + one_over_B_ = NonuniformTableLinear(press, B_inv); + viscosity_ = NonuniformTableLinear(press, visc); // Dumping the created tables. // static int count = 0; diff --git a/opm/core/fluid/blackoil/SinglePvtDead.hpp b/opm/core/fluid/blackoil/SinglePvtDead.hpp index 90a7d695..db86b39d 100644 --- a/opm/core/fluid/blackoil/SinglePvtDead.hpp +++ b/opm/core/fluid/blackoil/SinglePvtDead.hpp @@ -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). @@ -22,7 +22,7 @@ #include -#include +#include #include namespace Opm @@ -36,9 +36,9 @@ namespace Opm class SinglePvtDead : public SinglePvtInterface { public: - typedef std::vector > > table_t; - SinglePvtDead(const table_t& pvd_table, const int samples = 16); - virtual ~SinglePvtDead(); + typedef std::vector > > table_t; + SinglePvtDead(const table_t& pvd_table); + virtual ~SinglePvtDead(); /// Viscosity as a function of p and z. virtual void mu(const int n, @@ -72,12 +72,12 @@ namespace Opm double* output_R, double* output_dRdp) const; private: - // PVT properties of dry gas or dead oil - UniformTableLinear one_over_B_; - UniformTableLinear viscosity_; + // PVT properties of dry gas or dead oil + NonuniformTableLinear one_over_B_; + NonuniformTableLinear viscosity_; }; } -#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED +#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED diff --git a/opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp b/opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp new file mode 100644 index 00000000..132ffa51 --- /dev/null +++ b/opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp @@ -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 . +*/ + +#include +#include +#include + +// Extra includes for debug dumping of tables. +// #include +// #include +// #include + +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 press(sz); + std::vector B_inv(sz); + std::vector 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(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); + } + +} diff --git a/opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp b/opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp new file mode 100644 index 00000000..18d7ecfe --- /dev/null +++ b/opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp @@ -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 . +*/ + +#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED +#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED + + +#include +#include +#include + +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 > > 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 one_over_B_; + UniformTableLinear viscosity_; + }; + +} + +#endif // OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED + diff --git a/opm/core/fluid/blackoil/SinglePvtInterface.hpp b/opm/core/fluid/blackoil/SinglePvtInterface.hpp index d9369d33..7bcbd44b 100644 --- a/opm/core/fluid/blackoil/SinglePvtInterface.hpp +++ b/opm/core/fluid/blackoil/SinglePvtInterface.hpp @@ -32,7 +32,7 @@ namespace Opm public: SinglePvtInterface(); - virtual ~SinglePvtInterface(); + virtual ~SinglePvtInterface(); /// \param[in] num_phases The number of active phases. /// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases diff --git a/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp b/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp index 2526d34c..a6785467 100644 --- a/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp +++ b/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp @@ -1,13 +1,13 @@ //=========================================================================== -// -// File: MiscibilityLiveGas.cpp -// -// Created: Wed Feb 10 09:21:53 2010 -// +// +// File: MiscibilityLiveGas.cpp +// +// Created: Wed Feb 10 09:21:53 2010 +// // Author: Bjørn Spjelkavik -// +// // Revision: $Id$ -// +// //=========================================================================== /* Copyright 2010 SINTEF ICT, Applied Mathematics. @@ -47,37 +47,37 @@ namespace Opm /// Constructor SinglePvtLiveGas::SinglePvtLiveGas(const table_t& pvtg) { - // GAS, PVTG - const int region_number = 0; - if (pvtg.size() != 1) { - THROW("More than one PVD-region"); - } - saturated_gas_table_.resize(4); - const int sz = pvtg[region_number].size(); - for (int k=0; k<4; ++k) { - saturated_gas_table_[k].resize(sz); - } + // GAS, PVTG + const int region_number = 0; + if (pvtg.size() != 1) { + THROW("More than one PVD-region"); + } + saturated_gas_table_.resize(4); + const int sz = pvtg[region_number].size(); + for (int k=0; k<4; ++k) { + saturated_gas_table_[k].resize(sz); + } - for (int i=0; i saturated_gas_table_[0][ltp]) { - return linearInterpolationExtrap(undersat_gas_tables_[ltp][0], - undersat_gas_tables_[ltp][item], - maxR); - } + // Extrapolate from last table section + int ltp = saturated_gas_table_[0].size() - 1; + if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) { + return linearInterpolationExtrap(undersat_gas_tables_[ltp][0], + undersat_gas_tables_[ltp][item], + maxR); + } - // Interpolate between table sections - double w = (press - saturated_gas_table_[0][is]) / - (saturated_gas_table_[0][is+1] - - saturated_gas_table_[0][is]); - if (undersat_gas_tables_[is][0].size() < 2) { - double val = saturated_gas_table_[item][is] + - w*(saturated_gas_table_[item][is+1] - - saturated_gas_table_[item][is]); - return val; - } - double val1 = - linearInterpolationExtrap(undersat_gas_tables_[is][0], - undersat_gas_tables_[is][item], - maxR); - double val2 = - linearInterpolationExtrap(undersat_gas_tables_[is+1][0], - undersat_gas_tables_[is+1][item], - maxR); - double val = val1 + w*(val2 - val1); - return val; - } - } + // Interpolate between table sections + double w = (press - saturated_gas_table_[0][is]) / + (saturated_gas_table_[0][is+1] - + saturated_gas_table_[0][is]); + if (undersat_gas_tables_[is][0].size() < 2) { + double val = saturated_gas_table_[item][is] + + w*(saturated_gas_table_[item][is+1] - + saturated_gas_table_[item][is]); + return val; + } + double val1 = + linearInterpolationExtrap(undersat_gas_tables_[is][0], + undersat_gas_tables_[is][item], + maxR); + double val2 = + linearInterpolationExtrap(undersat_gas_tables_[is+1][0], + undersat_gas_tables_[is+1][item], + maxR); + double val = val1 + w*(val2 - val1); + return val; + } + } } diff --git a/opm/core/fluid/blackoil/SinglePvtLiveGas.hpp b/opm/core/fluid/blackoil/SinglePvtLiveGas.hpp index 3eff52f4..c7761551 100644 --- a/opm/core/fluid/blackoil/SinglePvtLiveGas.hpp +++ b/opm/core/fluid/blackoil/SinglePvtLiveGas.hpp @@ -33,10 +33,10 @@ namespace Opm class SinglePvtLiveGas : public SinglePvtInterface { public: - typedef std::vector > > table_t; + typedef std::vector > > table_t; - SinglePvtLiveGas(const table_t& pvto); - virtual ~SinglePvtLiveGas(); + SinglePvtLiveGas(const table_t& pvto); + virtual ~SinglePvtLiveGas(); /// Viscosity as a function of p and z. virtual void mu(const int n, @@ -76,14 +76,14 @@ namespace Opm double evalR(double press, const double* surfvol) const; void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const; - // item: 1=>B 2=>mu; - double miscible_gas(const double press, + // item: 1=>B 2=>mu; + double miscible_gas(const double press, const double* surfvol, const int item, - const bool deriv = false) const; - // PVT properties of wet gas (with vaporised oil) - std::vector > saturated_gas_table_; - std::vector > > undersat_gas_tables_; + const bool deriv = false) const; + // PVT properties of wet gas (with vaporised oil) + std::vector > saturated_gas_table_; + std::vector > > undersat_gas_tables_; }; diff --git a/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp b/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp index 2a393839..990ae114 100644 --- a/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp +++ b/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp @@ -38,122 +38,122 @@ namespace Opm /// Constructor SinglePvtLiveOil::SinglePvtLiveOil(const table_t& pvto) { - // OIL, PVTO - const int region_number = 0; - if (pvto.size() != 1) { - THROW("More than one PVD-region"); - } - saturated_oil_table_.resize(4); - const int sz = pvto[region_number].size(); - for (int k=0; k<4; ++k) { - saturated_oil_table_[k].resize(sz); - } - for (int i=0; i 1) { - iPrev = i; - continue; - } - - bool flagPrev = (iPrev >= 0); - bool flagNext = true; - if (iNext < i) { - iPrev = iNext; - flagPrev = true; - iNext = i+1; - while (undersat_oil_tables_[iNext][0].size() < 2) { - ++iNext; - } - } - double slopePrevBinv = 0.0; - double slopePrevVisc = 0.0; - double slopeNextBinv = 0.0; - double slopeNextVisc = 0.0; - while (flagPrev || flagNext) { - double pressure0 = undersat_oil_tables_[i][0].back(); - double pressure = 1.0e47; - if (flagPrev) { - std::vector::iterator itPrev = upper_bound(undersat_oil_tables_[iPrev][0].begin(), - undersat_oil_tables_[iPrev][0].end(),pressure0+1.); - if (itPrev == undersat_oil_tables_[iPrev][0].end()) { - --itPrev; // Extrapolation ... - } else if (itPrev == undersat_oil_tables_[iPrev][0].begin()) { - ++itPrev; - } - if (itPrev == undersat_oil_tables_[iPrev][0].end()-1) { - flagPrev = false; // Last data set for "prev" ... - } - double dPPrev = *itPrev - *(itPrev-1); - pressure = *itPrev; - 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; - slopePrevVisc = (undersat_oil_tables_[iPrev][2][index] - undersat_oil_tables_[iPrev][2][index-1])/dPPrev; - } - if (flagNext) { - std::vector::iterator itNext = upper_bound(undersat_oil_tables_[iNext][0].begin(), - undersat_oil_tables_[iNext][0].end(),pressure0+1.); - if (itNext == undersat_oil_tables_[iNext][0].end()) { - --itNext; // Extrapolation ... - } else if (itNext == undersat_oil_tables_[iNext][0].begin()) { - ++itNext; - } - if (itNext == undersat_oil_tables_[iNext][0].end()-1) { - flagNext = false; // Last data set for "next" ... - } - double dPNext = *itNext - *(itNext-1); - if (flagPrev) { - pressure = std::min(pressure,*itNext); - } else { - pressure = *itNext; - } - 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; - slopeNextVisc = (undersat_oil_tables_[iNext][2][index] - undersat_oil_tables_[iNext][2][index-1])/dPNext; - } - double dP = pressure - pressure0; - if (iPrev >= 0) { - double w = (saturated_oil_table_[3][i] - 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][1].push_back(undersat_oil_tables_[i][1].back() + - dP*(slopePrevBinv+w*(slopeNextBinv-slopePrevBinv))); - undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back() + - dP*(slopePrevVisc+w*(slopeNextVisc-slopePrevVisc))); - } else { - 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][2].push_back(undersat_oil_tables_[i][2].back()+dP*slopeNextVisc); - } - } - } + // OIL, PVTO + const int region_number = 0; + if (pvto.size() != 1) { + THROW("More than one PVD-region"); + } + saturated_oil_table_.resize(4); + const int sz = pvto[region_number].size(); + for (int k=0; k<4; ++k) { + saturated_oil_table_[k].resize(sz); + } + for (int i=0; i 1) { + iPrev = i; + continue; + } + + bool flagPrev = (iPrev >= 0); + bool flagNext = true; + if (iNext < i) { + iPrev = iNext; + flagPrev = true; + iNext = i+1; + while (undersat_oil_tables_[iNext][0].size() < 2) { + ++iNext; + } + } + double slopePrevBinv = 0.0; + double slopePrevVisc = 0.0; + double slopeNextBinv = 0.0; + double slopeNextVisc = 0.0; + while (flagPrev || flagNext) { + double pressure0 = undersat_oil_tables_[i][0].back(); + double pressure = 1.0e47; + if (flagPrev) { + std::vector::iterator itPrev = upper_bound(undersat_oil_tables_[iPrev][0].begin(), + undersat_oil_tables_[iPrev][0].end(),pressure0+1.); + if (itPrev == undersat_oil_tables_[iPrev][0].end()) { + --itPrev; // Extrapolation ... + } else if (itPrev == undersat_oil_tables_[iPrev][0].begin()) { + ++itPrev; + } + if (itPrev == undersat_oil_tables_[iPrev][0].end()-1) { + flagPrev = false; // Last data set for "prev" ... + } + double dPPrev = *itPrev - *(itPrev-1); + pressure = *itPrev; + 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; + slopePrevVisc = (undersat_oil_tables_[iPrev][2][index] - undersat_oil_tables_[iPrev][2][index-1])/dPPrev; + } + if (flagNext) { + std::vector::iterator itNext = upper_bound(undersat_oil_tables_[iNext][0].begin(), + undersat_oil_tables_[iNext][0].end(),pressure0+1.); + if (itNext == undersat_oil_tables_[iNext][0].end()) { + --itNext; // Extrapolation ... + } else if (itNext == undersat_oil_tables_[iNext][0].begin()) { + ++itNext; + } + if (itNext == undersat_oil_tables_[iNext][0].end()-1) { + flagNext = false; // Last data set for "next" ... + } + double dPNext = *itNext - *(itNext-1); + if (flagPrev) { + pressure = std::min(pressure,*itNext); + } else { + pressure = *itNext; + } + 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; + slopeNextVisc = (undersat_oil_tables_[iNext][2][index] - undersat_oil_tables_[iNext][2][index-1])/dPNext; + } + double dP = pressure - pressure0; + if (iPrev >= 0) { + double w = (saturated_oil_table_[3][i] - 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][1].push_back(undersat_oil_tables_[i][1].back() + + dP*(slopePrevBinv+w*(slopeNextBinv-slopePrevBinv))); + undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back() + + dP*(slopePrevVisc+w*(slopeNextVisc-slopePrevVisc))); + } else { + 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][2].push_back(undersat_oil_tables_[i][2].back()+dP*slopeNextVisc); + } + } + } } /// Destructor. @@ -238,30 +238,30 @@ namespace Opm double SinglePvtLiveOil::evalB(double press, const double* surfvol) const { // 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, double& Bval, double& dBdpval) const { - Bval = evalB(press, surfvol); - dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true); + Bval = evalB(press, surfvol); + dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true); } double SinglePvtLiveOil::evalR(double press, const double* surfvol) const { if (surfvol[phase_pos_[Vapour]] == 0.0) { return 0.0; - } - double Rval = linearInterpolationExtrap(saturated_oil_table_[0], - saturated_oil_table_[3], press); - double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; - if (Rval < maxR ) { // Saturated case - return Rval; - } else { - return maxR; // Undersaturated case - } + } + double Rval = linearInterpolationExtrap(saturated_oil_table_[0], + saturated_oil_table_[3], press); + double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; + if (Rval < maxR ) { // Saturated case + return Rval; + } else { + return maxR; // Undersaturated case + } } void SinglePvtLiveOil::evalRDeriv(const double press, const double* surfvol, @@ -272,19 +272,19 @@ namespace Opm dRdpval = 0.0; return; } - Rval = linearInterpolationExtrap(saturated_oil_table_[0], + Rval = linearInterpolationExtrap(saturated_oil_table_[0], saturated_oil_table_[3], press); - double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; - if (Rval < maxR ) { + double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; + if (Rval < maxR ) { // Saturated case - dRdpval = linearInterpolDerivative(saturated_oil_table_[0], - saturated_oil_table_[3], - press); - } else { + dRdpval = linearInterpolDerivative(saturated_oil_table_[0], + saturated_oil_table_[3], + press); + } else { // Undersaturated case Rval = maxR; - dRdpval = 0.0; - } + dRdpval = 0.0; + } } @@ -293,57 +293,57 @@ namespace Opm const int item, const bool deriv) const { - int section; - double Rval = linearInterpolationExtrap(saturated_oil_table_[0], + int section; + double Rval = linearInterpolationExtrap(saturated_oil_table_[0], saturated_oil_table_[3], press, section); - double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; - if (deriv) { - if (Rval < maxR ) { // Saturated case - return linearInterpolDerivative(saturated_oil_table_[0], - saturated_oil_table_[item], - press); - } 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 + double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]]; + if (deriv) { + if (Rval < maxR ) { // Saturated case + return linearInterpolDerivative(saturated_oil_table_[0], + saturated_oil_table_[item], + press); + } 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]); + 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; - } - } + 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); + 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 diff --git a/opm/core/fluid/blackoil/SinglePvtLiveOil.hpp b/opm/core/fluid/blackoil/SinglePvtLiveOil.hpp index f8be6353..cef9595e 100644 --- a/opm/core/fluid/blackoil/SinglePvtLiveOil.hpp +++ b/opm/core/fluid/blackoil/SinglePvtLiveOil.hpp @@ -34,10 +34,10 @@ namespace Opm class SinglePvtLiveOil : public SinglePvtInterface { public: - typedef std::vector > > table_t; + typedef std::vector > > table_t; - SinglePvtLiveOil(const table_t& pvto); - virtual ~SinglePvtLiveOil(); + SinglePvtLiveOil(const table_t& pvto); + virtual ~SinglePvtLiveOil(); /// Viscosity as a function of p and z. virtual void mu(const int n, @@ -77,15 +77,15 @@ namespace Opm double evalR(double press, const double* surfvol) const; void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const; - // item: 1=>1/B 2=>mu; - double miscible_oil(const double press, + // item: 1=>1/B 2=>mu; + double miscible_oil(const double press, const double* surfvol, const int item, - const bool deriv = false) const; + const bool deriv = false) const; - // PVT properties of live oil (with dissolved gas) - std::vector > saturated_oil_table_; - std::vector > > undersat_oil_tables_; + // PVT properties of live oil (with dissolved gas) + std::vector > saturated_oil_table_; + std::vector > > undersat_oil_tables_; }; } diff --git a/opm/core/utility/NonuniformTableLinear.hpp b/opm/core/utility/NonuniformTableLinear.hpp new file mode 100644 index 00000000..71334460 --- /dev/null +++ b/opm/core/utility/NonuniformTableLinear.hpp @@ -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 . +*/ + +#ifndef OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED +#define OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED + +#include +#include +#include +#include + +#include +#include + +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 + 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& x_values, + const std::vector& y_values); + + /// @brief Get the domain. + /// @return the domain as a pair of doubles. + std::pair domain(); + + /// @brief Rescale the domain. + /// @param new_domain the new domain as a pair of doubles. + void rescaleDomain(std::pair 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 x_values_; + std::vector y_values_; + mutable std::vector x_values_reversed_; + mutable std::vector 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 + 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 + inline + NonuniformTableLinear + ::NonuniformTableLinear() + : left_(ClosestValue), right_(ClosestValue) + { + } + + template + inline + NonuniformTableLinear + ::NonuniformTableLinear(const std::vector& x_values, + const std::vector& y_values) + : x_values_(x_values), y_values_(y_values), + left_(ClosestValue), right_(ClosestValue) + { + ASSERT(isNondecreasing(x_values.begin(), x_values.end())); + } + + template + inline std::pair + NonuniformTableLinear + ::domain() + { + return std::make_pair(x_values_[0], x_values_.back()); + } + + template + inline void + NonuniformTableLinear + ::rescaleDomain(std::pair 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 + inline double + NonuniformTableLinear + ::operator()(const double x) const + { + return Opm::linearInterpolation(x_values_, y_values_, x); + } + + template + inline double + NonuniformTableLinear + ::derivative(const double x) const + { + return Opm::linearInterpolationDerivative(x_values_, y_values_, x); + } + + template + inline double + NonuniformTableLinear + ::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 + inline bool + NonuniformTableLinear + ::operator==(const NonuniformTableLinear& other) const + { + return x_values_ == other.x_values_ + && y_values_ == other.y_values_ + && left_ == other.left_ + && right_ == other.right_; + } + + template + inline void + NonuniformTableLinear + ::setLeftPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + left_ = rp; + } + + template + inline void + NonuniformTableLinear + ::setRightPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + right_ = rp; + } + +} // namespace Opm + +#endif // OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED diff --git a/tests/bo_resprop_test.cpp b/tests/bo_resprop_test.cpp index 5bf3aea1..02357e75 100644 --- a/tests/bo_resprop_test.cpp +++ b/tests/bo_resprop_test.cpp @@ -43,7 +43,7 @@ int main(int argc, char** argv) UnstructuredGrid grid; grid.number_of_cells = 1; grid.global_cell = NULL; - Opm::BlackoilPropertiesFromDeck props(deck, grid); + Opm::BlackoilPropertiesFromDeck props(deck, grid, param); const int n = 1; double p[n] = { 150e5 };