diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp index be8211c15..2496d98ef 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 d879036ef..44a76013b 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 1f34e89e2..0eba8ffb5 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -19,6 +19,7 @@ #include #include + namespace Opm { @@ -26,15 +27,19 @@ namespace Opm const UnstructuredGrid& grid, bool init_rock) { - if (init_rock){ + if (init_rock){ 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, @@ -45,13 +50,56 @@ namespace Opm if(init_rock){ rock_.init(deck, grid); } - const int samples = param.getDefault("dead_tab_size", 1025); - pvt_.init(deck, samples); - satprops_.init(deck, grid, 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() << ")."); + const int pvt_samples = param.getDefault("pvt_tab_size", 200); + pvt_.init(deck, pvt_samples); + + // 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"); + if (sat_samples > 1) { + if (threephase_model == "stone2") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "simple") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else { + THROW("Unknown threephase_model: " << threephase_model); + } + } else { + if (threephase_model == "stone2") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "simple") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else { + THROW("Unknown threephase_model: " << threephase_model); + } + } + + if (pvt_.numPhases() != satprops_->numPhases()) { + THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); } } @@ -256,7 +304,7 @@ namespace Opm double* kr, double* dkrds) const { - satprops_.relperm(n, s, cells, kr, dkrds); + satprops_->relperm(n, s, cells, kr, dkrds); } @@ -275,7 +323,7 @@ namespace Opm double* pc, double* dpcds) const { - satprops_.capPress(n, s, cells, pc, dpcds); + satprops_->capPress(n, s, cells, pc, dpcds); } @@ -291,7 +339,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 782407dc7..f78e2cd71 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -27,6 +27,7 @@ #include #include #include +#include struct UnstructuredGrid; @@ -40,7 +41,7 @@ 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, @@ -48,12 +49,15 @@ namespace Opm /// 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, @@ -164,9 +168,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. @@ -179,7 +183,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 dd8a857fb..501c8a406 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 0ce3c88aa..ab68017c1 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 b90b0876f..72c39b1b0 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 570154bf6..9aa690980 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 68623e5cc..6680eaa49 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 a5025e8a1..d3aaaa7b0 100644 --- a/opm/core/fluid/IncompPropertiesInterface.hpp +++ b/opm/core/fluid/IncompPropertiesInterface.hpp @@ -62,7 +62,7 @@ namespace Opm /// \return Array of P viscosity values. virtual const double* viscosity() const = 0; - /// Densities of fluid phases at surface conditions. + /// Densities of fluid phases at reservoir conditions. /// \return Array of P density values. virtual const double* density() const = 0; @@ -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 3d7d7bb2d..5220502ad 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 d1a3e9f3d..3e30e807c 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 62cc1b231..34a4ba002 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 acbabaa49..01bbeed1d 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.hpp b/opm/core/fluid/RockBasic.hpp index 965637807..2adc3ffd6 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 104861013..af65d852b 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 7ba8903fd..ec7db9a6a 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 65027e425..63e71a6e4 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/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp index 760f934d4..10ddf23dd 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 789e0a351..2f90672d8 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.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index bd70d8029..ea2acb342 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -19,12 +19,14 @@ #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED + +#include #include #include -#include #include #include #include +#include #include struct UnstructuredGrid; @@ -32,23 +34,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 +93,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 +117,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 000000000..f7a80b453 --- /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 000000000..30f277cc1 --- /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.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 191c6fa89..8633fd9ee 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/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 3c4435e8e..a3f9583ee 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -26,6 +26,10 @@ #include +// Silence compatibility warning from DUNE headers since we don't use +// the deprecated member anyway (in this compilation unit) +#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1 + // TODO: clean up includes. #include #include diff --git a/opm/core/newwells.h b/opm/core/newwells.h index f16675b6c..20cfce61d 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -235,6 +235,19 @@ void destroy_wells(struct Wells *W); +/** + * Create a deep-copy (i.e., clone) of an existing Wells object, including its + * controls. + * + * @param[in] W Existing Wells object. + * @return Complete clone of the input object. Dispose of resources using + * function destroy_wells() when no longer needed. Returns @c NULL in case of + * allocation failure. + */ +struct Wells * +clone_wells(const struct Wells *W); + + #ifdef __cplusplus } #endif diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 954e07fa4..8edb41a0d 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -36,6 +36,7 @@ #include #include #include +#include namespace Opm { @@ -122,7 +123,7 @@ namespace Opm WellState& well_state) { const int nc = grid_.number_of_cells; - const int nw = wells_->number_of_wells; + const int nw = (wells_ != 0) ? wells_->number_of_wells : 0; // Set up dynamic data. computePerSolveDynamicData(dt, state, well_state); @@ -446,30 +447,54 @@ namespace Opm /// Compute per-iteration dynamic properties for wells. void CompressibleTpfa::computeWellDynamicData(const double /*dt*/, const BlackoilState& /*state*/, - const WellState& /*well_state*/) + const WellState& well_state) { // These are the variables that get computed by this function: // // std::vector wellperf_A_; // std::vector wellperf_phasemob_; const int np = props_.numPhases(); - const int nw = wells_->number_of_wells; - const int nperf = wells_->well_connpos[nw]; + const int nw = (wells_ != 0) ? wells_->number_of_wells : 0; + const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0; wellperf_A_.resize(nperf*np*np); wellperf_phasemob_.resize(nperf*np); // The A matrix is set equal to the perforation grid cells' - // matrix, for both injectors and producers. + // matrix for producers, computed from bhp and injection + // component fractions from // The mobilities are set equal to the perforation grid cells' - // mobilities, for both injectors and producers. + // mobilities for producers. + std::vector mu(np); for (int w = 0; w < nw; ++w) { + bool producer = (wells_->type[w] == PRODUCER); + const double* comp_frac = &wells_->comp_frac[np*w]; for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) { const int c = wells_->well_cells[j]; - const double* cA = &cell_A_[np*np*c]; double* wpA = &wellperf_A_[np*np*j]; - std::copy(cA, cA + np*np, wpA); - const double* cM = &cell_phasemob_[np*c]; double* wpM = &wellperf_phasemob_[np*j]; - std::copy(cM, cM + np, wpM); + if (producer) { + const double* cA = &cell_A_[np*np*c]; + std::copy(cA, cA + np*np, wpA); + const double* cM = &cell_phasemob_[np*c]; + std::copy(cM, cM + np, wpM); + } else { + const double bhp = well_state.bhp()[w]; + double perf_p = bhp; + for (int phase = 0; phase < np; ++phase) { + perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase]; + } + // Hack warning: comp_frac is used as a component + // surface-volume variable in calls to matrix() and + // viscosity(), but as a saturation in the call to + // relperm(). This is probably ok as long as injectors + // only inject pure fluids. + props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL); + props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL); + ASSERT(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6); + props_.relperm (1, comp_frac, &c, wpM , NULL); + for (int phase = 0; phase < np; ++phase) { + wpM[phase] /= mu[phase]; + } + } } } } @@ -487,9 +512,9 @@ namespace Opm const double* z = &state.surfacevol()[0]; UnstructuredGrid* gg = const_cast(&grid_); CompletionData completion_data; - completion_data.gpot = &wellperf_gpot_[0]; - completion_data.A = &wellperf_A_[0]; - completion_data.phasemob = &wellperf_phasemob_[0]; + completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0; + completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0; + completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(wells_); wells_tmp.data = &completion_data; @@ -574,9 +599,9 @@ namespace Opm { UnstructuredGrid* gg = const_cast(&grid_); CompletionData completion_data; - completion_data.gpot = const_cast(&wellperf_gpot_[0]); - completion_data.A = const_cast(&wellperf_A_[0]); - completion_data.phasemob = const_cast(&wellperf_phasemob_[0]); + completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast(&wellperf_gpot_[0]) : 0; + completion_data.A = ! wellperf_A_.empty() ? const_cast(&wellperf_A_[0]) : 0; + completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast(&wellperf_phasemob_[0]) : 0; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(wells_); wells_tmp.data = &completion_data; @@ -584,6 +609,9 @@ namespace Opm forces.wells = &wells_tmp; forces.src = NULL; + double* wpress = ! well_state.bhp ().empty() ? & well_state.bhp ()[0] : 0; + double* wflux = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0; + cfs_tpfa_res_flux(gg, &forces, props_.numPhases(), @@ -592,9 +620,9 @@ namespace Opm &face_phasemob_[0], &face_gravcap_[0], &state.pressure()[0], - &well_state.bhp()[0], + wpress, &state.faceflux()[0], - &well_state.perfRates()[0]); + wflux); cfs_tpfa_res_fpress(gg, props_.numPhases(), &htrans_[0], @@ -604,6 +632,23 @@ namespace Opm &state.pressure()[0], &state.faceflux()[0], &state.facepressure()[0]); + + // Compute well perforation pressures (not done by the C code). + if (wells_ != 0) { + const int nw = wells_->number_of_wells; + const int np = props_.numPhases(); + for (int w = 0; w < nw; ++w) { + const double* comp_frac = &wells_->comp_frac[np*w]; + for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) { + const double bhp = well_state.bhp()[w]; + double perf_p = bhp; + for (int phase = 0; phase < np; ++phase) { + perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase]; + } + well_state.perfPress()[j] = perf_p; + } + } + } } } // namespace Opm diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index eb8dc8f34..bcbacec28 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid *G , nnu = G->number_of_cells; nwperf = 0; - if (wells != NULL) { assert (wells->W != NULL); + if ((wells != NULL) && (wells->W != NULL)) { nnu += wells->W->number_of_wells; nwperf = wells->W->well_connpos[ wells->W->number_of_wells ]; } @@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid *G , /* ---------------------------------------------------------------------- */ { int f, c1, c2, w, i, nc, nnu; + int wells_present; size_t nnz; struct CSRMatrix *A; nc = nnu = G->number_of_cells; - if (wells != NULL) { - assert (wells->W != NULL); + + wells_present = (wells != NULL) && (wells->W != NULL); + if (wells_present) { nnu += wells->W->number_of_wells; } @@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid *G , } } - if (wells != NULL) { + if (wells_present) { /* Well <-> cell connections */ struct Wells *W = wells->W; @@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid *G , } } - if (wells != NULL) { + if (wells_present) { /* Fill well <-> cell connections */ struct Wells *W = wells->W; @@ -741,6 +743,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt, } +/* ---------------------------------------------------------------------- */ +static void +welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h, + double *res, double *w2c, double *w2w) +/* ---------------------------------------------------------------------- */ +{ + int p; + double fwi; + const double *dpflux_w; + + /* Sum reservoir phase flux derivatives set by + * compute_darcyflux_and_deriv(). */ + dpflux_w = h->pimpl->flux_work + (1 * np); + for (p = 0, fwi = 0.0; p < np; p++) { + fwi += dpflux_w[ p ]; + } + + *res = 0.0; + *w2c = 0.0; + *w2w = fwi; +} + + /* ---------------------------------------------------------------------- */ static void welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h, @@ -795,7 +820,34 @@ welleq_coeff_resv(int np, struct cfs_tpfa_res_data *h, /* ---------------------------------------------------------------------- */ static void -assemble_completion_to_well(int w, int c, int nc, int np, +welleq_coeff_surfrate(int i, int np, struct cfs_tpfa_res_data *h, + struct WellControls *ctrl, + double *res, double *w2c, double *w2w) +/* ---------------------------------------------------------------------- */ +{ + int p; + double f; + const double *pflux, *dpflux_w, *dpflux_c, *distr; + + pflux = h->pimpl->compflux_p + (i * (1 * np)); + dpflux_w = h->pimpl->compflux_deriv_p + (i * (2 * np)); + dpflux_c = dpflux_w + (1 * (1 * np)); + distr = ctrl->distr + (ctrl->current * (1 * np)); + + *res = *w2c = *w2w = 0.0; + for (p = 0; p < np; p++) { + f = distr[ p ]; + + *res += f * pflux [ p ]; + *w2w += f * dpflux_w[ p ]; + *w2c += f * dpflux_c[ p ]; + } +} + + +/* ---------------------------------------------------------------------- */ +static void +assemble_completion_to_well(int i, int w, int c, int nc, int np, double pw, double dt, struct cfs_tpfa_res_wells *wells, struct cfs_tpfa_res_data *h ) @@ -815,23 +867,25 @@ assemble_completion_to_well(int w, int c, int nc, int np, ctrl = W->ctrls[ w ]; if (ctrl->current < 0) { - assert (0); /* Shut wells currently not supported */ + /* Interpreting a negative current control index to mean a shut well */ + welleq_coeff_shut(np, h, &res, &w2c, &w2w); } + else { + switch (ctrl->type[ ctrl->current ]) { + case BHP : + welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ], + h, &res, &w2c, &w2w); + break; - switch (ctrl->type[ ctrl->current ]) { - case BHP : - welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ], - h, &res, &w2c, &w2w); - break; + case RESERVOIR_RATE: + assert (W->number_of_phases == np); + welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w); + break; - case RESERVOIR_RATE: - assert (W->number_of_phases == np); - welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w); - break; - - case SURFACE_RATE: - assert (0); /* Surface rate currently not supported */ - break; + case SURFACE_RATE: + welleq_coeff_surfrate(i, np, h, ctrl, &res, &w2c, &w2w); + break; + } } /* Assemble completion contributions */ @@ -854,7 +908,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , struct cfs_tpfa_res_data *h ) { int w, i, c, np, np2, nc; - int is_neumann; + int is_neumann, is_open; double pw, dp; double *WI, *gpot, *pmobp; const double *Ac, *dAc; @@ -876,6 +930,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , for (w = i = 0; w < W->number_of_wells; w++) { pw = wpress[ w ]; + is_open = W->ctrls[w]->current >= 0; for (; i < W->well_connpos[w + 1]; i++, gpot += np, pmobp += np) { @@ -888,14 +943,16 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , init_completion_contrib(i, np, Ac, dAc, h->pimpl); - assemble_completion_to_cell(c, nc + w, np, dt, h); + if (is_open) { + assemble_completion_to_cell(c, nc + w, np, dt, h); + } /* Prepare for RESV controls */ compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot, h->pimpl->flux_work, h->pimpl->flux_work + np); - assemble_completion_to_well(w, c, nc, np, pw, dt, wells, h); + assemble_completion_to_well(i, w, c, nc, np, pw, dt, wells, h); } ctrl = W->ctrls[ w ]; @@ -1127,8 +1184,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , nf = G->number_of_faces; nwperf = 0; - if (wells != NULL) { - assert (wells->W != NULL); + if ((wells != NULL) && (wells->W != NULL)) { nwperf = wells->W->well_connpos[ wells->W->number_of_wells ]; } @@ -1194,7 +1250,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , assemble_cell_contrib(G, c, h); } - if ((forces != NULL) && (forces->wells != NULL)) { + if ((forces != NULL) && + (forces->wells != NULL) && + (forces->wells->W != NULL)) { compute_well_compflux_and_deriv(forces->wells, cq->nphases, cpress, wpress, h->pimpl); @@ -1297,8 +1355,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G , { compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux); - if ((forces != NULL) && (forces->wells != NULL) && - (wpress != NULL) && (wflux != NULL)) { + if ((forces != NULL) && + (forces->wells != NULL) && + (forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) { compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux); } diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 59d9f0602..981ccd9ba 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -37,12 +37,20 @@ namespace Opm if (wells) { const int nw = wells->number_of_wells; bhp_.resize(nw); - // Initialize bhp to be pressure in first perforation cell. + // Initialize bhp to be target pressure + // if bhp-controlled well, otherwise set + // to pressure in first perforation cell. for (int w = 0; w < nw; ++w) { - const int cell = wells->well_cells[wells->well_connpos[w]]; - bhp_[w] = state.pressure()[cell]; + const WellControls* ctrl = wells->ctrls[w]; + if (ctrl->type[ctrl->current] == BHP) { + bhp_[w] = ctrl->target[ctrl->current]; + } else { + const int cell = wells->well_cells[wells->well_connpos[w]]; + bhp_[w] = state.pressure()[cell]; + } } - perfrates_.resize(wells->well_connpos[nw]); + perfrates_.resize(wells->well_connpos[nw], 0.0); + perfpress_.resize(wells->well_connpos[nw], -1e100); } } @@ -54,9 +62,14 @@ namespace Opm std::vector& perfRates() { return perfrates_; } const std::vector& perfRates() const { return perfrates_; } + /// One pressure per well connection. + std::vector& perfPress() { return perfpress_; } + const std::vector& perfPress() const { return perfpress_; } + private: std::vector bhp_; std::vector perfrates_; + std::vector perfpress_; }; } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 01882053b..bad712630 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -152,8 +152,8 @@ namespace Opm B_cell = 1.0/tm.A_[np*np*cell + 0]; double src_flux = -tm.source_[cell]; bool src_is_inflow = src_flux < 0.0; - influx = src_is_inflow ? B_cell*src_flux : 0.0; - outflux = !src_is_inflow ? B_cell*src_flux : 0.0; + influx = src_is_inflow ? B_cell* src_flux : 0.0; + outflux = !src_is_inflow ? src_flux : 0.0; comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell]; dtpv = tm.dt_/tm.porevolume0_[cell]; for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { @@ -349,7 +349,7 @@ namespace Opm gf[1] = gravflux[pos]; } s0 = tm.saturation_[cell]; - dtpv = tm.dt_/tm.porevolume0_[cell]; + dtpv = tm.dt_/tm.porevolume_[cell]; } double operator()(double s) const @@ -380,8 +380,8 @@ namespace Opm { double sat[2] = { s, 1.0 - s }; props_.relperm(1, sat, &cell, mob, 0); - mob[0] /= visc_[0]; - mob[1] /= visc_[1]; + mob[0] /= visc_[2*cell + 0]; + mob[1] /= visc_[2*cell + 1]; } @@ -407,7 +407,7 @@ namespace Opm { // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] - // But b_w,i * rho_w,S = rho_w,i, which we conmpute with a call to props_.density(). + // But b_w,i * rho_w,S = rho_w,i, which we compute with a call to props_.density(). // We assume that we already have stored T_ij in trans_. // We also assume that the A_ matrices are updated from an earlier call to solve(). const int nc = grid_.number_of_cells; @@ -437,16 +437,15 @@ namespace Opm void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector& cells, - const int pos, - const double* gravflux) + const int pos, + const double* gravflux) { const int cell = cells[pos]; GravityResidual res(*this, cells, pos, gravflux); if (std::fabs(res(saturation_[cell])) > tol_) { int iters_used; - saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); + saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); } - saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]); mobility(saturation_[cell], cell, &mob_[2*cell]); } @@ -506,8 +505,6 @@ namespace Opm void TransportModelCompressibleTwophase::solveGravity(const std::vector >& columns, - const double* pressure, - const double* porevolume0, const double dt, std::vector& saturation, std::vector& surfacevol) @@ -522,18 +519,22 @@ namespace Opm cells[c] = c; } mob_.resize(2*nc); - std::vector boths; - Opm::toBothSat(saturation, boths); - props_.relperm(cells.size(), &boths[0], &cells[0], &mob_[0], 0); - props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); - for (int c = 0; c < nc; ++c) { - mob_[2*c + 0] /= visc_[2*c + 0]; - mob_[2*c + 1] /= visc_[2*c + 1]; + + // props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0); + + // props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); + // for (int c = 0; c < nc; ++c) { + // mob_[2*c + 0] /= visc_[2*c + 0]; + // mob_[2*c + 1] /= visc_[2*c + 1]; + // } + + const int np = props_.numPhases(); + for (int cell = 0; cell < nc; ++cell) { + mobility(saturation_[cell], cell, &mob_[np*cell]); } // Set up other variables. - porevolume0_ = porevolume0; dt_ = dt; toWaterSat(saturation, saturation_); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 4ee93e0af..147d7da10 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -74,13 +74,10 @@ namespace Opm /// vertical stack, that do not interact with other columns (for /// gravity segregation. /// \param[in] columns Vector of cell-columns. - /// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] dt Time step. /// \param[in, out] saturation Phase saturations. /// \param[in, out] surfacevol Surface volume densities for each phase. void solveGravity(const std::vector >& columns, - const double* pressure, - const double* porevolume0, const double dt, std::vector& saturation, std::vector& surfacevol); diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index f8e138a34..8a90ae34f 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g std::vector sequence(grid.number_of_cells); std::vector components(grid.number_of_cells + 1); int ncomponents; + time::StopWatch clock; + clock.start(); compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents); + clock.stop(); + std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl; // Invoke appropriate solve method for each interdependent component. for (int comp = 0; comp < ncomponents; ++comp) { diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp new file mode 100644 index 000000000..8d93f2480 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -0,0 +1,122 @@ +/* + 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 +#include +#include +#include + +namespace Opm +{ + + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid) + : grid_(grid) + { + } + + + + + /// Solve for time-of-flight. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. + /// \param[out] tof Array of time-of-flight values. + void TransportModelTracerTof::solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + std::vector& tof) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + THROW("Sources do not sum to zero: " << cum_src); + } +#endif + tof.resize(grid_.number_of_cells); + std::fill(tof.begin(), tof.end(), 0.0); + tof_ = &tof[0]; + reorderAndTransport(grid_, darcyflux); + } + + + + + void TransportModelTracerTof::solveSingleCell(const int cell) + { + // Compute flux terms. + // Sources have zero tof, and therefore do not contribute + // to upwind_term. Sinks on the other hand, must be added + // to the downwind_flux (note sign change resulting from + // different sign conventions: pos. source is injection, + // pos. flux is outflow). + double upwind_term = 0.0; + double downwind_flux = std::max(-source_[cell], 0.0); + for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { + int f = grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == grid_.face_cells[2*f]) { + flux = darcyflux_[f]; + other = grid_.face_cells[2*f+1]; + } else { + flux =-darcyflux_[f]; + other = grid_.face_cells[2*f]; + } + // Add flux to upwind_term or downwind_flux, if interior. + if (other != -1) { + if (flux < 0.0) { + upwind_term += flux*tof_[other]; + } else { + downwind_flux += flux; + } + } + } + + // Compute tof. + tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux; + } + + + + + void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells) + { + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + for (int i = 0; i < num_cells; ++i) { + solveSingleCell(cells[i]); + } + } + + + + +} // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp new file mode 100644 index 000000000..270ac328a --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -0,0 +1,74 @@ +/* + 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_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED + +#include +#include +#include +#include +struct UnstructuredGrid; + +namespace Opm +{ + + class IncompPropertiesInterface; + + /// Implements a first-order finite volume solver for + /// (single-phase) time-of-flight using reordering. + /// The equation solved is: + /// v \cdot \grad\tau = \phi + /// where v is the fluid velocity, \tau is time-of-flight and + /// \phi is the porosity. This is a boundary value problem, where + /// \tau is specified to be zero on all inflow boundaries. + class TransportModelTracerTof : public TransportModelInterface + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTof(const UnstructuredGrid& grid); + + /// Solve for time-of-flight. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. + /// \param[out] tof Array of time-of-flight values. + void solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + std::vector& tof); + + private: + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + + private: + const UnstructuredGrid& grid_; + const double* darcyflux_; // one flux per grid face + const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell + double* tof_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp new file mode 100644 index 000000000..404452e0e --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -0,0 +1,657 @@ +/* + 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 +#include +#include +#include +#include +#include + +namespace Opm +{ + + + // --------------- Helpers for TransportModelTracerTofDiscGal --------------- + + + + /// A class providing discontinuous Galerkin basis functions. + struct DGBasis + { + static int numBasisFunc(const int dimensions, + const int degree) + { + switch (dimensions) { + case 1: + return degree + 1; + case 2: + return (degree + 2)*(degree + 1)/2; + case 3: + return (degree + 3)*(degree + 2)*(degree + 1)/6; + default: + THROW("Dimensions must be 1, 2 or 3."); + } + } + + /// Evaluate all nonzero basis functions at x, + /// writing to f_x. The array f_x must have + /// size numBasisFunc(grid.dimensions, degree). + /// + /// The basis functions are the following + /// Degree 0: 1. + /// Degree 1: x - xc, y - yc, z - zc etc. + /// Further degrees await development. + static void eval(const UnstructuredGrid& grid, + const int cell, + const int degree, + const double* x, + double* f_x) + { + const int dim = grid.dimensions; + const double* cc = grid.cell_centroids + dim*cell; + // Note intentional fallthrough in this switch statement! + switch (degree) { + case 1: + for (int ix = 0; ix < dim; ++ix) { + f_x[1 + ix] = x[ix] - cc[ix]; + } + case 0: + f_x[0] = 1; + break; + default: + THROW("Maximum degree is 1 for now."); + } + } + + /// Evaluate gradients of all nonzero basis functions at x, + /// writing to grad_f_x. The array grad_f_x must have size + /// numBasisFunc(grid.dimensions, degree) * grid.dimensions. + /// The components of the first basis function + /// gradient come before the components of the second etc. + static void evalGrad(const UnstructuredGrid& grid, + const int /*cell*/, + const int degree, + const double* /*x*/, + double* grad_f_x) + { + const int dim = grid.dimensions; + const int num_basis = numBasisFunc(dim, degree); + std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); + if (degree > 1) { + THROW("Maximum degree is 1 for now."); + } else if (degree == 1) { + for (int ix = 0; ix < dim; ++ix) { + grad_f_x[dim*(ix + 1) + ix] = 1.0; + } + } + } + }; + + + + static void cross(const double* a, const double* b, double* res) + { + res[0] = a[1]*b[2] - a[2]*b[1]; + res[1] = a[2]*b[0] - a[0]*b[2]; + res[2] = a[0]*b[1] - a[1]*b[0]; + } + + + + + static double triangleArea3d(const double* p0, + const double* p1, + const double* p2) + { + double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; + double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; + double cr[3]; + cross(a, b, cr); + return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]); + } + + + + + /// Calculates the determinant of a 3 x 3 matrix, represented as + /// three three-dimensional arrays. + static double determinantOf(const double* a0, + const double* a1, + const double* a2) + { + return + a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - + a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + + a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); + } + + + + + /// Computes the volume of a tetrahedron consisting of 4 vertices + /// with 3-dimensional coordinates + static double tetVolume(const double* p0, + const double* p1, + const double* p2, + const double* p3) + { + double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; + double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; + double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] }; + return std::fabs(determinantOf(a, b, c) / 6.0); + } + + + + + /// A class providing numerical quadrature for cells. + /// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). + /// Note that this class does multiply weights by cell volume, + /// so weights always sum to cell volume. + /// Degree 1 method: + /// Midpoint (centroid) method. + /// n = 1, w_0 = cell volume, x_0 = cell centroid + /// Degree 2 method: + /// Based on subdivision of each cell face into triangles + /// with the face centroid as a common vertex, and then + /// subdividing the cell into tetrahedra with the cell + /// centroid as a common vertex. Then apply the tetrahedron + /// rule with the following 4 nodes (uniform weights): + /// a = 0.138196601125010515179541316563436 + /// x_i has all barycentric coordinates = a, except for + /// the i'th coordinate which is = 1 - 3a. + /// This rule is from http://nines.cs.kuleuven.be/ecf, + /// it is the second degree 2 4-point rule for tets, + /// referenced to Stroud(1971). + /// The tetrahedra are numbered T_{i,j}, and are given by the + /// cell centroid C, the face centroid FC_i, and two nodes + /// of face i: FN_{i,j}, FN_{i,j+1}. + class CellQuadrature + { + public: + CellQuadrature(const UnstructuredGrid& grid, + const int cell, + const int degree) + : grid_(grid), cell_(cell), degree_(degree) + { + if (degree > 2) { + THROW("CellQuadrature exact for polynomial degrees > 1 not implemented."); + } + if (degree == 2) { + // Prepare subdivision. + } + } + + int numQuadPts() const + { + if (degree_ < 2) { + return 1; + } + // Degree 2 case. + int sumnodes = 0; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + } + return 4*sumnodes; + } + + void quadPtCoord(const int index, double* coord) const + { + const int dim = grid_.dimensions; + const double* cc = grid_.cell_centroids + dim*cell_; + if (degree_ < 2) { + std::copy(cc, cc + dim, coord); + return; + } + // Degree 2 case. + int tetindex = index / 4; + const int subindex = index % 4; + const double* nc = grid_.node_coordinates; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + if (nfn <= tetindex) { + // Our tet is not associated with this face. + tetindex -= nfn; + continue; + } + const double* fc = grid_.face_centroids + dim*face; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; + const int node0 = fnodes[tetindex]; + const int node1 = fnodes[(tetindex + 1) % nfn]; + const double* n0c = nc + dim*node0; + const double* n1c = nc + dim*node1; + const double a = 0.138196601125010515179541316563436; + // Barycentric coordinates of our point in the tet. + double baryc[4] = { a, a, a, a }; + baryc[subindex] = 1.0 - 3.0*a; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd]; + } + return; + } + THROW("Should never reach this point."); + } + + double quadPtWeight(const int index) const + { + if (degree_ < 2) { + return grid_.cell_volumes[cell_]; + } + // Degree 2 case. + const int dim = grid_.dimensions; + const double* cc = grid_.cell_centroids + dim*cell_; + int tetindex = index / 4; + const double* nc = grid_.node_coordinates; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + if (nfn <= tetindex) { + // Our tet is not associated with this face. + tetindex -= nfn; + continue; + } + const double* fc = grid_.face_centroids + dim*face; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; + const int node0 = fnodes[tetindex]; + const int node1 = fnodes[(tetindex + 1) % nfn]; + const double* n0c = nc + dim*node0; + const double* n1c = nc + dim*node1; + return 0.25*tetVolume(cc, fc, n0c, n1c); + } + THROW("Should never reach this point."); + } + + private: + const UnstructuredGrid& grid_; + const int cell_; + const int degree_; + }; + + + + + + /// A class providing numerical quadrature for faces. + /// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). + /// Note that this class does multiply weights by face area, + /// so weights always sum to face area. + /// Degree 1 method: + /// Midpoint (centroid) method. + /// n = 1, w_0 = face area, x_0 = face centroid + /// Degree 2 method: + /// Based on subdivision of the face into triangles, + /// with the centroid as a common vertex, and the triangle + /// edge midpoint rule. + /// Triangle i consists of the centroid C, nodes N_i and N_{i+1}. + /// Its area is A_i. + /// n = 2 * nn (nn = num nodes in face) + /// For i = 0..(nn-1): + /// w_i = 1/3 A_i. + /// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i + /// x_i = (N_i + N_{i+1})/2 + /// x_{nn+i} = (C + N_i)/2 + /// All N and A indices are interpreted cyclic, modulus nn. + class FaceQuadrature + { + public: + FaceQuadrature(const UnstructuredGrid& grid, + const int face, + const int degree) + : grid_(grid), face_(face), degree_(degree) + { + if (grid_.dimensions != 3) { + THROW("FaceQuadrature only implemented for 3D case."); + } + if (degree_ > 2) { + THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented."); + } + } + + int numQuadPts() const + { + if (degree_ < 2) { + return 1; + } + // Degree 2 case. + return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]); + } + + void quadPtCoord(const int index, double* coord) const + { + const int dim = grid_.dimensions; + const double* fc = grid_.face_centroids + dim*face_; + if (degree_ < 2) { + std::copy(fc, fc + dim, coord); + return; + } + // Degree 2 case. + const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; + const double* nc = grid_.node_coordinates; + if (index < nn) { + // Boundary edge midpoint. + const int node0 = fnodes[index]; + const int node1 = fnodes[(index + 1)%nn]; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]); + } + } else { + // Interiour edge midpoint. + // Recall that index is now in [nn, 2*nn). + const int node = fnodes[index - nn]; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]); + } + } + } + + double quadPtWeight(const int index) const + { + if (degree_ < 2) { + return grid_.face_areas[face_]; + } + // Degree 2 case. + const int dim = grid_.dimensions; + const double* fc = grid_.face_centroids + dim*face_; + const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; + const double* nc = grid_.node_coordinates; + if (index < nn) { + // Boundary edge midpoint. + const int node0 = fnodes[index]; + const int node1 = fnodes[(index + 1)%nn]; + const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); + return area / 3.0; + } else { + // Interiour edge midpoint. + // Recall that index is now in [nn, 2*nn). + const int node0 = fnodes[(index - 1) % nn]; + const int node1 = fnodes[index - nn]; + const int node2 = fnodes[(index + 1) % nn]; + const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); + const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc); + return (area0 + area1) / 3.0; + } + } + + private: + const UnstructuredGrid& grid_; + const int face_; + const int degree_; + }; + + + + + + // --------------- Methods of TransportModelTracerTofDiscGal --------------- + + + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] use_cvi If true, use corner point velocity interpolation. + /// Otherwise, use the basic constant interpolation. + TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, + const bool use_cvi) + : grid_(grid), + coord_(grid.dimensions), + velocity_(grid.dimensions) + { + if (use_cvi) { + velocity_interpolation_.reset(new VelocityInterpolationECVI(grid)); + } else { + velocity_interpolation_.reset(new VelocityInterpolationConstant(grid)); + } + } + + + + + /// Solve for time-of-flight. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. + /// \param[in] degree Polynomial degree of DG basis functions used. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + const int degree, + std::vector& tof_coeff) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + THROW("Sources do not sum to zero: " << cum_src); + } +#endif + degree_ = degree; + const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + tof_coeff.resize(num_basis*grid_.number_of_cells); + std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); + tof_coeff_ = &tof_coeff[0]; + rhs_.resize(num_basis); + jac_.resize(num_basis*num_basis); + orig_jac_.resize(num_basis*num_basis); + basis_.resize(num_basis); + basis_nb_.resize(num_basis); + grad_basis_.resize(num_basis*grid_.dimensions); + velocity_interpolation_->setupFluxes(darcyflux); + reorderAndTransport(grid_, darcyflux); + } + + + + + void TransportModelTracerTofDiscGal::solveSingleCell(const int cell) + { + // Residual: + // For each cell K, basis function b_j (spanning V_h), + // writing the solution u_h|K = \sum_i c_i b_i + // Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx + // + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds + // - \int_K \phi b_j + // This is linear in c_i, so we do not need any nonlinear iterations. + // We assemble the jacobian and the right-hand side. The residual is + // equal to Res = Jac*c - rhs, and we compute rhs directly. + + const int dim = grid_.dimensions; + const int num_basis = DGBasis::numBasisFunc(dim, degree_); + + std::fill(rhs_.begin(), rhs_.end(), 0.0); + std::fill(jac_.begin(), jac_.end(), 0.0); + + // Compute cell residual contribution. + // Note: Assumes that \int_K b_j = 0 for all j > 0 + rhs_[0] += porevolume_[cell]; + + // Compute upstream residual contribution. + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + int upstream_cell = -1; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face+1]; + } else { + flux = -darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face]; + } + if (upstream_cell < 0) { + // This is an outer boundary. Assumed tof = 0 on inflow, so no contribution. + continue; + } + if (flux >= 0.0) { + // This is an outflow boundary. + continue; + } + // Do quadrature over the face to compute + // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds + // (where u_h^{ext} is the upstream unknown (tof)). + const double normal_velocity = flux / grid_.face_areas[face]; + FaceQuadrature quad(grid_, face, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]); + const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(), + tof_coeff_ + num_basis*upstream_cell, 0.0); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j]; + } + } + } + + // Compute cell jacobian contribution. We use Fortran ordering + // for jac_, i.e. rows cycling fastest. + { + CellQuadrature quad(grid_, cell, 2*degree_ - 1); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // b_i (v \cdot \grad b_j) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]); + velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + for (int dd = 0; dd < dim; ++dd) { + jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd]; + } + } + } + } + } + + // Compute downstream jacobian contribution from faces. + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux <= 0.0) { + // This is an inflow boundary. + continue; + } + // Do quadrature over the face to compute + // \int_{\partial K} b_i (v(x) \cdot n) b_j ds + const double normal_velocity = flux / grid_.face_areas[face]; + FaceQuadrature quad(grid_, face, 2*degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // u^ext flux B (B = {b_j}) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j]; + } + } + } + } + + // Compute downstream jacobian contribution from sink terms. + // Contribution from inflow sources would be + // similar to the contribution from upstream faces, but + // it is zero since we let all external inflow be associated + // with a zero tof. + if (source_[cell] < 0.0) { + // A sink. + const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. + const double flux_density = flux / grid_.cell_volumes[cell]; + // Do quadrature over the cell to compute + // \int_{K} b_i flux b_j dx + CellQuadrature quad(grid_, cell, 2*degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j]; + } + } + } + } + + // Solve linear equation. + MAT_SIZE_T n = num_basis; + MAT_SIZE_T nrhs = 1; + MAT_SIZE_T lda = num_basis; + std::vector piv(num_basis); + MAT_SIZE_T ldb = num_basis; + MAT_SIZE_T info = 0; + orig_jac_ = jac_; + dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); + if (info != 0) { + // Print the local matrix and rhs. + std::cerr << "Failed solving single-cell system Ax = b in cell " << cell + << " with A = \n"; + for (int row = 0; row < n; ++row) { + for (int col = 0; col < n; ++col) { + std::cerr << " " << orig_jac_[row + n*col]; + } + std::cerr << '\n'; + } + std::cerr << "and b = \n"; + for (int row = 0; row < n; ++row) { + std::cerr << " " << rhs_[row] << '\n'; + } + THROW("Lapack error: " << info << " encountered in cell " << cell); + } + // The solution ends up in rhs_, so we must copy it. + std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); + } + + + + + void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells) + { + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + for (int i = 0; i < num_cells; ++i) { + solveSingleCell(cells[i]); + } + } + + + + +} // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp new file mode 100644 index 000000000..9a2e3b9c2 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -0,0 +1,105 @@ +/* + 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_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED + +#include +#include +#include +#include +#include + +struct UnstructuredGrid; + +namespace Opm +{ + + class IncompPropertiesInterface; + class VelocityInterpolationInterface; + + /// Implements a discontinuous Galerkin solver for + /// (single-phase) time-of-flight using reordering. + /// The equation solved is: + /// v \cdot \grad\tau = \phi + /// where v is the fluid velocity, \tau is time-of-flight and + /// \phi is the porosity. This is a boundary value problem, where + /// \tau is specified to be zero on all inflow boundaries. + /// The user may specify the polynomial degree of the basis function space + /// used, but only degrees 0 and 1 are supported so far. + class TransportModelTracerTofDiscGal : public TransportModelInterface + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] use_cvi If true, use corner point velocity interpolation. + /// Otherwise, use the basic constant interpolation. + TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, + const bool use_cvi); + + + /// Solve for time-of-flight. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. + /// \param[in] degree Polynomial degree of DG basis functions used. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + void solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + const int degree, + std::vector& tof_coeff); + + private: + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + + private: + // Disable copying and assignment. + TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&); + TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&); + + const UnstructuredGrid& grid_; + boost::shared_ptr velocity_interpolation_; + const double* darcyflux_; // one flux per grid face + const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell + int degree_; + double* tof_coeff_; + std::vector rhs_; // single-cell right-hand-side + std::vector jac_; // single-cell jacobian + std::vector orig_jac_; // single-cell jacobian (copy) + // Below: storage for quantities needed by solveSingleCell(). + std::vector coord_; + std::vector basis_; + std::vector basis_nb_; + std::vector grad_basis_; + std::vector velocity_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 7d01ade0e..575e18a41 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include namespace Opm @@ -512,15 +513,23 @@ namespace Opm State& state) { const int num_phases = props.numPhases(); + const PhaseUsage pu = phaseUsageFromDeck(deck); + if (num_phases != pu.num_phases) { + THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, " + "found " << pu.num_phases << " phases in deck."); + } state.init(grid, num_phases); if (deck.hasField("EQUIL")) { if (num_phases != 2) { THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); } + if (pu.phase_used[BlackoilPhases::Vapour]) { + THROW("initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas)."); + } // Set saturations depending on oil-water contact. const EQUIL& equil= deck.getEQUIL(); if (equil.equil.size() != 1) { - THROW("No region support yet."); + THROW("initStateFromDeck(): No region support yet."); } const double woc = equil.equil[0].water_oil_contact_depth_; initWaterOilContact(grid, props, woc, WaterBelow, state); @@ -528,37 +537,64 @@ namespace Opm const double datum_z = equil.equil[0].datum_depth_; const double datum_p = equil.equil[0].datum_depth_pressure_; initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); - } else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) { - // Set saturations from SWAT, pressure from PRESSURE. + } else if (deck.hasField("PRESSURE")) { + // Set saturations from SWAT/SGAS, pressure from PRESSURE. std::vector& s = state.saturation(); std::vector& p = state.pressure(); - const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); const int num_cells = grid.number_of_cells; if (num_phases == 2) { - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; - s[2*c] = sw_deck[c_deck]; - s[2*c + 1] = 1.0 - s[2*c]; - p[c] = p_deck[c_deck]; + if (!pu.phase_used[BlackoilPhases::Aqua]) { + // oil-gas: we require SGAS + if (!deck.hasField("SGAS")) { + THROW("initStateFromDeck(): missing SGAS keyword in 2-phase init"); + } + const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + gpos] = sg_deck[c_deck]; + s[2*c + opos] = 1.0 - sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + // water-oil or water-gas: we require SWAT + if (!deck.hasField("SWAT")) { + THROW("initStateFromDeck(): missing SWAT keyword in 2-phase init"); + } + const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int nwpos = (wpos + 1) % 2; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + wpos] = sw_deck[c_deck]; + s[2*c + nwpos] = 1.0 - sw_deck[c_deck]; + p[c] = p_deck[c_deck]; + } } } else if (num_phases == 3) { - if (!deck.hasField("SGAS")) { - THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found)."); + const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS"); + if (!has_swat_sgas) { + THROW("initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init."); } + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); for (int c = 0; c < num_cells; ++c) { int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; - s[3*c] = sw_deck[c_deck]; - s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); - s[3*c + 2] = sg_deck[c_deck]; + s[3*c + wpos] = sw_deck[c_deck]; + s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); + s[3*c + gpos] = sg_deck[c_deck]; p[c] = p_deck[c_deck]; } } else { THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); } } else { - THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); + THROW("initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS."); } // Finally, init face pressures. diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index aeafbfca7..73e10ac90 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -374,7 +374,10 @@ namespace Opm if (perf_rate > 0.0) { // perf_rate is a total inflow rate, we want a water rate. if (wells->type[w] != INJECTOR) { - std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl; + std::cout << "**** Warning: crossflow in well " + << w << " perf " << perf - wells->well_connpos[w] + << " ignored. Rate was " + << perf_rate/Opm::unit::day << " m^3/day." << std::endl; perf_rate = 0.0; } else { ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6); @@ -550,6 +553,7 @@ namespace Opm { const int np = wells.number_of_phases; const int nw = wells.number_of_wells; + ASSERT(int(flow_rates_per_well_cell.size()) == wells.well_connpos[nw]); phase_flow_per_well.resize(nw * np); for (int wix = 0; wix < nw; ++wix) { for (int phase = 0; phase < np; ++phase) { diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 6f4f75565..ae806e9a4 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -21,7 +21,10 @@ #include #include #include +#include #include +#include +#include #include #include #include @@ -31,53 +34,74 @@ namespace Opm { - /// @brief Computes injected and produced volumes of all phases. + /// @brief Computes injected and produced surface volumes of all phases. /// Note 1: assumes that only the first phase is injected. /// Note 2: assumes that transport has been done with an /// implicit method, i.e. that the current state /// gives the mobilities used for the preceding timestep. - /// @param[in] props fluid and rock properties. - /// @param[in] p pressure (one value per cell) - /// @param[in] z surface-volume values (for all P phases) - /// @param[in] s saturation values (for all P phases) - /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. - /// @param[in] dt timestep used - /// @param[out] injected must point to a valid array with P elements, - /// where P = s.size()/src.size(). - /// @param[out] produced must also point to a valid array with P elements. + /// Note 3: Gives surface volume values, not reservoir volumes + /// (as the incompressible version of the function does). + /// Also, assumes that transport_src is given in surface volumes + /// for injector terms! + /// @param[in] props fluid and rock properties. + /// @param[in] state state variables (pressure, sat, surfvol) + /// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const BlackoilPropertiesInterface& props, - const std::vector& press, - const std::vector& z, - const std::vector& s, - const std::vector& src, + const BlackoilState& state, + const std::vector& transport_src, const double dt, double* injected, double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); + const int num_cells = transport_src.size(); + if (props.numCells() != num_cells) { + THROW("Size of transport_src vector does not match number of cells in props."); } + const int np = props.numPhases(); + if (int(state.saturation().size()) != num_cells*np) { + THROW("Sizes of state vectors do not match number of cells."); + } + const std::vector& press = state.pressure(); + const std::vector& s = state.saturation(); + const std::vector& z = state.surfacevol(); std::fill(injected, injected + np, 0.0); std::fill(produced, produced + np, 0.0); std::vector visc(np); std::vector mob(np); + std::vector A(np*np); + std::vector prod_resv_phase(np); + std::vector prod_surfvol(np); for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; + if (transport_src[c] > 0.0) { + // Inflowing transport source is a surface volume flux + // for the first phase. + injected[0] += transport_src[c]*dt; + } else if (transport_src[c] < 0.0) { + // Outflowing transport source is a total reservoir + // volume flux. + const double flux = -transport_src[c]*dt; const double* sat = &s[np*c]; props.relperm(1, sat, &c, &mob[0], 0); props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); + props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0); double totmob = 0.0; for (int p = 0; p < np; ++p) { mob[p] /= visc[p]; totmob += mob[p]; } + std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0); for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; + prod_resv_phase[p] = (mob[p]/totmob)*flux; + for (int q = 0; q < np; ++q) { + prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p]; + } + } + for (int p = 0; p < np; ++p) { + produced[p] += prod_surfvol[p]; } } } @@ -251,4 +275,58 @@ namespace Opm } } + + /// Compute two-phase transport source terms from well terms. + /// Note: Unlike the incompressible version of this function, + /// this version computes surface volume injection rates, + /// production rates are still total reservoir volumes. + /// \param[in] props Fluid and rock properties. + /// \param[in] wells Wells data structure. + /// \param[in] well_state Well pressures and fluxes. + /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign: + /// (+) positive inflow of first (water) phase (surface volume), + /// (-) negative total outflow of both phases (reservoir volume). + void computeTransportSource(const BlackoilPropertiesInterface& props, + const Wells* wells, + const WellState& well_state, + std::vector& transport_src) + { + int nc = props.numCells(); + transport_src.clear(); + transport_src.resize(nc, 0.0); + // Well contributions. + if (wells) { + const int nw = wells->number_of_wells; + const int np = wells->number_of_phases; + if (np != 2) { + THROW("computeTransportSource() requires a 2 phase case."); + } + std::vector A(np*np); + for (int w = 0; w < nw; ++w) { + const double* comp_frac = wells->comp_frac + np*w; + for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { + const int perf_cell = wells->well_cells[perf]; + double perf_rate = well_state.perfRates()[perf]; + if (perf_rate > 0.0) { + // perf_rate is a total inflow reservoir rate, we want a surface water rate. + if (wells->type[w] != INJECTOR) { + std::cout << "**** Warning: crossflow in well " + << w << " perf " << perf - wells->well_connpos[w] + << " ignored. Reservoir rate was " + << perf_rate/Opm::unit::day << " m^3/day." << std::endl; + perf_rate = 0.0; + } else { + ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6); + perf_rate *= comp_frac[0]; // Water reservoir volume rate. + props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0); + perf_rate *= A[0]; // Water surface volume rate. + } + } + transport_src[perf_cell] += perf_rate; + } + } + } + } + + } // namespace Opm diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 565acc6ef..f4949ce2d 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -22,36 +22,40 @@ #include -struct UnstructuredGrid; +struct Wells; namespace Opm { class BlackoilPropertiesInterface; + class BlackoilState; + class WellState; - /// @brief Computes injected and produced volumes of all phases. + + /// @brief Computes injected and produced surface volumes of all phases. /// Note 1: assumes that only the first phase is injected. /// Note 2: assumes that transport has been done with an /// implicit method, i.e. that the current state /// gives the mobilities used for the preceding timestep. - /// @param[in] props fluid and rock properties. - /// @param[in] p pressure (one value per cell) - /// @param[in] z surface-volume values (for all P phases) - /// @param[in] s saturation values (for all P phases) - /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. - /// @param[in] dt timestep used - /// @param[out] injected must point to a valid array with P elements, - /// where P = s.size()/src.size(). - /// @param[out] produced must also point to a valid array with P elements. + /// Note 3: Gives surface volume values, not reservoir volumes + /// (as the incompressible version of the function does). + /// Also, assumes that transport_src is given in surface volumes + /// for injector terms! + /// @param[in] props fluid and rock properties. + /// @param[in] state state variables (pressure, sat, surfvol) + /// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const BlackoilPropertiesInterface& props, - const std::vector& p, - const std::vector& z, - const std::vector& s, - const std::vector& src, + const BlackoilState& state, + const std::vector& transport_src, const double dt, double* injected, double* produced); + /// @brief Computes total mobility for a set of saturation values. /// @param[in] props rock and fluid properties /// @param[in] cells cells with which the saturation values are associated @@ -66,6 +70,7 @@ namespace Opm const std::vector& s, std::vector& totmob); + /// @brief Computes total mobility and omega for a set of saturation values. /// @param[in] props rock and fluid properties /// @param[in] cells cells with which the saturation values are associated @@ -131,6 +136,22 @@ namespace Opm const double* saturation, double* surfacevol); + + /// Compute two-phase transport source terms from well terms. + /// Note: Unlike the incompressible version of this function, + /// this version computes surface volume injection rates, + /// production rates are still total reservoir volumes. + /// \param[in] props Fluid and rock properties. + /// \param[in] wells Wells data structure. + /// \param[in] well_state Well pressures and fluxes. + /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign: + /// (+) positive inflow of first (water) phase (surface volume), + /// (-) negative total outflow of both phases (reservoir volume). + void computeTransportSource(const BlackoilPropertiesInterface& props, + const Wells* wells, + const WellState& well_state, + std::vector& transport_src); + } // namespace Opm #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED diff --git a/opm/core/utility/writeVtkData.hpp b/opm/core/utility/writeVtkData.hpp index 200b2b261..96cdfac22 100644 --- a/opm/core/utility/writeVtkData.hpp +++ b/opm/core/utility/writeVtkData.hpp @@ -26,25 +26,23 @@ #include #include #include +#include struct UnstructuredGrid; namespace Opm { - /// Intended to map strings (giving the output field names) to data. - typedef std::map*> DataMap; - /// Vtk output for cartesian grids. void writeVtkData(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const DataMap& data, - std::ostream& os); + const std::tr1::array& cell_size, + const DataMap& data, + std::ostream& os); /// Vtk output for general grids. void writeVtkData(const UnstructuredGrid& grid, - const DataMap& data, - std::ostream& os); + const DataMap& data, + std::ostream& os); } // namespace Opm #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index afa17f12d..8ffb30273 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -572,6 +572,7 @@ namespace Opm break; } const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase); + const double total_reinjected = - total_produced; // Production negative, injection positive const double my_guide_rate = injectionGuideRate(true); for (size_t i = 0; i < children_.size(); ++i) { // Apply for all children. @@ -580,11 +581,11 @@ namespace Opm const double children_guide_rate = children_[i]->injectionGuideRate(true); #ifdef DIRTY_WELLCTRL_HACK children_[i]->applyInjGroupControl(InjectionSpecification::RESV, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_, false); #else children_[i]->applyInjGroupControl(InjectionSpecification::RATE, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_, false); #endif } @@ -600,15 +601,15 @@ namespace Opm if (phaseUsage().phase_used[BlackoilPhases::Vapour]) { total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour); } - - const double my_guide_rate = injectionGuideRate(true); + const double total_reinjected = - total_produced; // Production negative, injection positive + const double my_guide_rate = injectionGuideRate(true); for (size_t i = 0; i < children_.size(); ++i) { // Apply for all children. // Note, we do _not_ want to call the applyProdGroupControl in this object, // as that would check if we're under group control, something we're not. const double children_guide_rate = children_[i]->injectionGuideRate(true); children_[i]->applyInjGroupControl(InjectionSpecification::RESV, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_, + (children_guide_rate / my_guide_rate) * total_reinjected * injSpec().voidage_replacment_fraction_, false); } @@ -760,16 +761,17 @@ namespace Opm wells_->ctrls[self_index_]->current = ~ wells_->ctrls[self_index_]->current; } } - + std::pair WellNode::getWorstOffending(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase, ProductionSpecification::ControlMode mode) { const int np = phaseUsage().num_phases; const int index = self_index_*np; - return std::make_pair(this, rateByMode(&well_reservoirrates_phase[index], - &well_surfacerates_phase[index], - mode)); + return std::pair(this, + rateByMode(&well_reservoirrates_phase[index], + &well_surfacerates_phase[index], + mode)); } void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode, @@ -781,7 +783,7 @@ namespace Opm && (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) { return; } - if (!wells_->type[self_index_] == INJECTOR) { + if (wells_->type[self_index_] != INJECTOR) { ASSERT(target == 0.0); return; } @@ -858,12 +860,12 @@ namespace Opm std::cout << "Returning" << std::endl; return; } - if (!wells_->type[self_index_] == PRODUCER) { + if (wells_->type[self_index_] != PRODUCER) { ASSERT(target == 0.0); return; } // We're a producer, so we need to negate the input - double ntarget = target; + double ntarget = -target; double distr[3] = { 0.0, 0.0, 0.0 }; const int* phase_pos = phaseUsage().phase_pos; diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 79387242d..f4c7f0692 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -230,6 +230,14 @@ namespace Opm + /// Construct from existing wells object. + WellsManager::WellsManager(struct Wells* W) + : w_(clone_wells(W)) + { + } + + + /// Construct wells from deck. WellsManager::WellsManager(const Opm::EclipseGridParser& deck, const UnstructuredGrid& grid, @@ -548,7 +556,7 @@ namespace Opm cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; } else if (wci_line.injector_type_[0] == 'G') { if (!pu.phase_used[BlackoilPhases::Vapour]) { - THROW("Water phase not used, yet found water-injecting well."); + THROW("Gas phase not used, yet found gas-injecting well."); } cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; } @@ -720,32 +728,31 @@ namespace Opm #endif if (deck.hasField("WELOPEN")) { - const WELOPEN& welopen = deck.getWELOPEN(); - - for (size_t i = 0; i < welopen.welopen.size(); ++i) { - WelopenLine line = welopen.welopen[i]; - std::string wellname = line.well_; - std::map::const_iterator it = well_names_to_index.find(wellname); - if (it == well_names_to_index.end()) { - THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); - } - int index = it->second; - if (line.openshutflag_ == "SHUT") { - // We currently don't care if the well is open or not. - /// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }? - //ASSERT(w_->ctrls[index]->current < 0); - } else if (line.openshutflag_ == "OPEN") { - //ASSERT(w_->ctrls[index]->current >= 0); - } else { - THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT."); - } - - // We revert back to it's original control. - // Note that this is OK as ~~ = id. - w_->ctrls[index]->current = ~w_->ctrls[index]->current; - - - } + const WELOPEN& welopen = deck.getWELOPEN(); + for (size_t i = 0; i < welopen.welopen.size(); ++i) { + WelopenLine line = welopen.welopen[i]; + std::string wellname = line.well_; + std::map::const_iterator it = well_names_to_index.find(wellname); + if (it == well_names_to_index.end()) { + THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); + } + const int index = it->second; + if (line.openshutflag_ == "SHUT") { + int& cur_ctrl = w_->ctrls[index]->current; + if (cur_ctrl >= 0) { + cur_ctrl = ~cur_ctrl; + } + ASSERT(w_->ctrls[index]->current < 0); + } else if (line.openshutflag_ == "OPEN") { + int& cur_ctrl = w_->ctrls[index]->current; + if (cur_ctrl < 0) { + cur_ctrl = ~cur_ctrl; + } + ASSERT(w_->ctrls[index]->current >= 0); + } else { + THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT."); + } + } } // Build the well_collection_ well group hierarchy. diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 650eaf371..19289c5d9 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -42,8 +42,12 @@ namespace Opm public: /// Default constructor -- no wells. WellsManager(); - /// Construct from mrst type output. - /// Wellmanger is not properly initialized + + /// Construct from existing wells object. + /// WellsManager is not properly initialised in the sense that the logic to + /// manage control switching does not exist. + /// + /// @param[in] W Existing wells object. WellsManager(struct Wells* W); /// Construct from input deck and grid. diff --git a/tests/test_wells.cpp b/tests/test_wells.cpp new file mode 100644 index 000000000..4f8ad7bb7 --- /dev/null +++ b/tests/test_wells.cpp @@ -0,0 +1,205 @@ +/* + 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 + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellsModuleTest +#include + +#include + +#include +#include +#include + +BOOST_AUTO_TEST_CASE(Construction) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr W(create_wells(nphases, nwells, nperfs), + destroy_wells); + + if (W) { + int cells[] = { 0, 9 }; + double WI = 1.0; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0], + &WI, "INJECTOR", W.get()); + + const double pfrac[] = { 0.0, 0.0 }; + const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1], + &WI, "PRODUCER", W.get()); + + if (ok0 && ok1) { + BOOST_CHECK_EQUAL(W->number_of_phases, nphases); + BOOST_CHECK_EQUAL(W->number_of_wells , nwells ); + + BOOST_CHECK_EQUAL(W->well_connpos[0], 0); + BOOST_CHECK_EQUAL(W->well_connpos[1], 1); + BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs); + + BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]); + BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]); + + BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI); + BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI); + + using std::string; + BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR")); + BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER")); + } + } +} + + +BOOST_AUTO_TEST_CASE(Controls) +{ + const int nphases = 2; + const int nwells = 1; + const int nperfs = 2; + + boost::shared_ptr W(create_wells(nphases, nwells, nperfs), + destroy_wells); + + if (W) { + int cells[] = { 0 , 9 }; + double WI [] = { 1.0, 1.0 }; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0], + &WI[0], "INJECTOR", W.get()); + + if (ok) { + const double distr[] = { 1.0, 0.0 }; + const bool ok1 = append_well_controls(BHP, 1, &distr[0], + 0, W.get()); + const bool ok2 = append_well_controls(SURFACE_RATE, 1, + &distr[0], 0, W.get()); + + if (ok1 && ok2) { + WellControls* ctrls = W->ctrls[0]; + + BOOST_CHECK_EQUAL(ctrls->num , 2); + BOOST_CHECK_EQUAL(ctrls->current, -1); + + set_current_control(0, 0, W.get()); + BOOST_CHECK_EQUAL(ctrls->current, 0); + + set_current_control(0, 1, W.get()); + BOOST_CHECK_EQUAL(ctrls->current, 1); + + BOOST_CHECK_EQUAL(ctrls->type[0], BHP); + BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE); + + BOOST_CHECK_EQUAL(ctrls->target[0], 1.0); + BOOST_CHECK_EQUAL(ctrls->target[1], 1.0); + } + } + } +} + + +BOOST_AUTO_TEST_CASE(Copy) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr W1(create_wells(nphases, nwells, nperfs), + destroy_wells); + boost::shared_ptr W2; + + if (W1) { + int cells[] = { 0, 9 }; + const double WI = 1.0; + const double ifrac[] = { 1.0, 0.0 }; + + const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0], + &WI, "INJECTOR", W1.get()); + + const double pfrac[] = { 0.0, 0.0 }; + const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1], + &WI, "PRODUCER", W1.get()); + + bool ok = ok0 && ok1; + for (int w = 0; ok && (w < W1->number_of_wells); ++w) { + const double distr[] = { 1.0, 0.0 }; + const bool okc1 = append_well_controls(BHP, 1, &distr[0], + w, W1.get()); + const bool okc2 = append_well_controls(SURFACE_RATE, 1, + &distr[0], w, + W1.get()); + + ok = okc1 && okc2; + } + + if (ok) { + W2.reset(clone_wells(W1.get()), destroy_wells); + } + } + + if (W2) { + BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases); + BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells ); + BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] ); + + for (int w = 0; w < W1->number_of_wells; ++w) { + using std::string; + BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w])); + BOOST_CHECK_EQUAL( W2->type[w] , W1->type[w] ); + + BOOST_CHECK_EQUAL(W2->well_connpos[w + 1], + W1->well_connpos[w + 1]); + + for (int j = W1->well_connpos[w]; + j < W1->well_connpos[w + 1]; ++j) { + BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]); + BOOST_CHECK_EQUAL(W2->WI [j], W1->WI [j]); + } + + BOOST_CHECK(W1->ctrls[w] != 0); + BOOST_CHECK(W2->ctrls[w] != 0); + + WellControls* c1 = W1->ctrls[w]; + WellControls* c2 = W2->ctrls[w]; + + BOOST_CHECK_EQUAL(c2->num , c1->num ); + BOOST_CHECK_EQUAL(c2->current, c1->current); + + for (int c = 0; c < c1->num; ++c) { + BOOST_CHECK_EQUAL(c2->type [c], c1->type [c]); + BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]); + + for (int p = 0; p < W1->number_of_phases; ++p) { + BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p], + c1->distr[c*W1->number_of_phases + p]); + } + } + } + } +}