diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index e7244d87c..fdc9519cf 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -28,10 +28,20 @@ namespace Opm { rock_.init(deck, grid); pvt_.init(deck, use_spline); - satprops_.init(deck, grid); - if (pvt_.numPhases() != satprops_.numPhases()) { + // Unfortunate lack of pointer smartness here... + if (use_spline) { + SaturationPropsFromDeck<>* ptr = new SaturationPropsFromDeck<>(); + ptr->init(deck, grid); + satprops_.reset(ptr); + } else { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + ptr->init(deck, grid); + satprops_.reset(ptr); + } + 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_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); } } @@ -236,7 +246,7 @@ namespace Opm double* kr, double* dkrds) const { - satprops_.relperm(n, s, cells, kr, dkrds); + satprops_->relperm(n, s, cells, kr, dkrds); } @@ -255,7 +265,7 @@ namespace Opm double* pc, double* dpcds) const { - satprops_.capPress(n, s, cells, pc, dpcds); + satprops_->capPress(n, s, cells, pc, dpcds); } @@ -271,7 +281,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 c434b7f4b..60dc1e754 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include struct UnstructuredGrid; @@ -166,7 +167,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/IncompPropertiesFromDeck.hpp b/opm/core/fluid/IncompPropertiesFromDeck.hpp index 68623e5cc..5820bb668 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.hpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.hpp @@ -135,7 +135,7 @@ namespace Opm private: RockFromDeck rock_; PvtPropertiesIncompFromDeck pvt_; - SaturationPropsFromDeck satprops_; + SaturationPropsFromDeck<> satprops_; }; diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index b83a12711..a978bbac0 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -20,8 +20,8 @@ #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED +#include #include -#include #include #include @@ -30,7 +30,23 @@ struct UnstructuredGrid; namespace Opm { - class SaturationPropsFromDeck : public BlackoilPhases + + /// Class storing saturation functions in a uniform table, + /// densely sampled from a monotone spline, + /// using linear interpolation. + class SatFuncSetUniform; + + /// Class storing saturation functions in a nonuniform table + /// using linear interpolation. + class SatFuncSetNonuniform; + + + /// Interface to saturation functions from deck. + /// Possible values for template argument: + /// SatFuncSetNonuniform, + /// SatFuncSetUniform (default). + template + class SaturationPropsFromDeck : public SaturationPropsInterface { public: /// Default constructor. @@ -88,30 +104,12 @@ namespace Opm private: PhaseUsage phase_usage_; - class SatFuncSet - { - public: - void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); - void evalKr(const double* s, double* kr) const; - void evalKrDeriv(const double* s, double* kr, double* dkrds) const; - void evalPc(const double* s, double* pc) const; - void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - UniformTableLinear krw_; - UniformTableLinear krow_; - UniformTableLinear pcow_; - UniformTableLinear krg_; - UniformTableLinear krog_; - UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - }; std::vector satfuncset_; std::vector cell_to_func_; // = SATNUM - 1 - const SatFuncSet& funcForCell(const int cell) const; + typedef SatFuncSet Funcs; + + const Funcs& funcForCell(const int cell) const; }; @@ -119,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..f64f12bcf --- /dev/null +++ b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp @@ -0,0 +1,271 @@ +/* + 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 +{ + + /// Class storing saturation functions in a uniform table, + /// densely sampled from a monotone spline, + /// using linear interpolation. + class SatFuncSetUniform : public BlackoilPhases + { + public: + void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); + void evalKr(const double* s, double* kr) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds) const; + void evalPc(const double* s, double* pc) const; + void evalPcDeriv(const double* s, double* pc, double* dpcds) const; + double smin_[PhaseUsage::MaxNumPhases]; + double smax_[PhaseUsage::MaxNumPhases]; + private: + PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. + UniformTableLinear krw_; + UniformTableLinear krow_; + UniformTableLinear pcow_; + UniformTableLinear krg_; + UniformTableLinear krog_; + UniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + + + + + /// Class storing saturation functions in a nonuniform table + /// using linear interpolation. + class SatFuncSetNonuniform : public BlackoilPhases + { + public: + void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); + void evalKr(const double* s, double* kr) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds) const; + void evalPc(const double* s, double* pc) const; + void evalPcDeriv(const double* s, double* pc, double* dpcds) const; + double smin_[PhaseUsage::MaxNumPhases]; + double smax_[PhaseUsage::MaxNumPhases]; + private: + PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. + NonuniformTableLinear krw_; + NonuniformTableLinear krow_; + NonuniformTableLinear pcow_; + NonuniformTableLinear krg_; + NonuniformTableLinear krog_; + NonuniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + + + + // ----------- Methods of SaturationPropsFromDeck --------- + + + /// Default constructor. + template + SaturationPropsFromDeck::SaturationPropsFromDeck() + { + } + + /// Initialize from deck. + template + void SaturationPropsFromDeck::init(const EclipseGridParser& deck, + const UnstructuredGrid& grid) + { + phase_usage_ = phaseUsageFromDeck(deck); + + // Extract input data. + // Oil phase should be active. + if (!phase_usage_.phase_used[Liquid]) { + THROW("SaturationPropsFromDeck::init() -- oil phase must be active."); + } + + // Obtain SATNUM, if it exists, and create cell_to_func_. + // Otherwise, let the cell_to_func_ mapping be just empty. + int satfuncs_expected = 1; + if (deck.hasField("SATNUM")) { + const std::vector& satnum = deck.getIntegerValue("SATNUM"); + satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); + const int num_cells = grid.number_of_cells; + cell_to_func_.resize(num_cells); + const int* gc = grid.global_cell; + for (int cell = 0; cell < num_cells; ++cell) { + const int deck_pos = (gc == NULL) ? cell : gc[cell]; + cell_to_func_[cell] = satnum[deck_pos] - 1; + } + } + + // Find number of tables, check for consistency. + enum { Uninitialized = -1 }; + int num_tables = Uninitialized; + if (phase_usage_.phase_used[Aqua]) { + const SWOF::table_t& swof_table = deck.getSWOF().swof_; + num_tables = swof_table.size(); + if (num_tables < satfuncs_expected) { + THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected); + } + } + if (phase_usage_.phase_used[Vapour]) { + const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; + int num_sgof_tables = sgof_table.size(); + if (num_sgof_tables < satfuncs_expected) { + THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected); + } + if (num_tables == Uninitialized) { + num_tables = num_sgof_tables; + } else if (num_tables != num_sgof_tables) { + THROW("Inconsistent number of tables in SWOF and SGOF."); + } + } + + // Initialize tables. + satfuncset_.resize(num_tables); + for (int table = 0; table < num_tables; ++table) { + satfuncset_[table].init(deck, table, phase_usage_); + } + } + + + + + /// \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..e55995677 --- /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