diff --git a/Makefile.am b/Makefile.am index c30a6544..cb3d794e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -51,6 +51,7 @@ opm/core/fluid/RockFromDeck.cpp \ opm/core/fluid/SaturationPropsBasic.cpp \ opm/core/fluid/SaturationPropsFromDeck.cpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \ +opm/core/fluid/blackoil/SinglePvtDead.cpp \ opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp \ opm/core/fluid/blackoil/SinglePvtInterface.cpp \ opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \ @@ -151,6 +152,7 @@ opm/core/fluid/SimpleFluid2p.hpp \ opm/core/fluid/blackoil/BlackoilPhases.hpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \ opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \ +opm/core/fluid/blackoil/SinglePvtDead.hpp \ opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp \ opm/core/fluid/blackoil/SinglePvtInterface.hpp \ opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \ @@ -226,6 +228,7 @@ opm/core/utility/ColumnExtract.hpp \ opm/core/utility/ErrorMacros.hpp \ opm/core/utility/Factory.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \ +opm/core/utility/NonuniformTableLinear.hpp \ opm/core/utility/RootFinders.hpp \ opm/core/utility/SparseTable.hpp \ opm/core/utility/SparseVector.hpp \ diff --git a/examples/sim_2p_comp_reorder.cpp b/examples/sim_2p_comp_reorder.cpp index 9f2b516b..705812c5 100644 --- a/examples/sim_2p_comp_reorder.cpp +++ b/examples/sim_2p_comp_reorder.cpp @@ -94,7 +94,8 @@ main(int argc, char** argv) // Grid init grid.reset(new GridManager(*deck)); // Rock and fluid init - props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid())); + bool props_use_spline = param.getDefault("props_use_spline", true); + props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), props_use_spline)); // check_well_controls = param.getDefault("check_well_controls", false); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Rock compressibility. diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index f7e0172f..bde2ecf9 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -175,7 +175,8 @@ main(int argc, char** argv) // Grid init grid.reset(new Opm::GridManager(deck)); // Rock and fluid init - props.reset(new Opm::BlackoilPropertiesFromDeck(deck, *grid->c_grid())); + bool props_use_spline = param.getDefault("props_use_spline", true); + props.reset(new BlackoilPropertiesFromDeck(deck, *grid->c_grid(), props_use_spline)); // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); check_well_controls = param.getDefault("check_well_controls", false); diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 004e6f88..e7244d87 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -23,10 +23,11 @@ namespace Opm { BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid) + const UnstructuredGrid& grid, + const bool use_spline) { rock_.init(deck, grid); - pvt_.init(deck); + pvt_.init(deck, use_spline); satprops_.init(deck, grid); if (pvt_.numPhases() != satprops_.numPhases()) { THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index a2a6c3de..c434b7f4 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -43,7 +43,8 @@ namespace Opm /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid, + const bool use_spline); /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp index c86724eb..dd8ce720 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp @@ -20,6 +20,7 @@ #include +#include #include #include #include @@ -39,7 +40,7 @@ namespace Opm } - void BlackoilPvtProperties::init(const EclipseGridParser& deck) + void BlackoilPvtProperties::init(const EclipseGridParser& deck, const bool use_spline) { typedef std::vector > > table_t; // If we need multiple regions, this class and the SinglePvt* classes must change. @@ -78,7 +79,11 @@ namespace Opm // Oil PVT if (phase_usage_.phase_used[Liquid]) { if (deck.hasField("PVDO")) { - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_)); + if (use_spline) { + props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_)); + } else { + props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_)); + } } else if (deck.hasField("PVTO")) { props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_)); } else if (deck.hasField("PVCDO")) { @@ -90,7 +95,11 @@ namespace Opm // Gas PVT if (phase_usage_.phase_used[Vapour]) { if (deck.hasField("PVDG")) { - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_)); + if (use_spline) { + props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_)); + } else { + props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_)); + } } else if (deck.hasField("PVTG")) { props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_)); } else { diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 627f266b..ca193a4d 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp @@ -47,7 +47,7 @@ namespace Opm BlackoilPvtProperties(); /// Initialize from deck. - void init(const EclipseGridParser& deck); + void init(const EclipseGridParser& deck, const bool use_spline); /// Number of active phases. int numPhases() const; diff --git a/opm/core/fluid/blackoil/SinglePvtDead.cpp b/opm/core/fluid/blackoil/SinglePvtDead.cpp new file mode 100644 index 00000000..f9bcd090 --- /dev/null +++ b/opm/core/fluid/blackoil/SinglePvtDead.cpp @@ -0,0 +1,127 @@ +/* + 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 + +// Extra includes for debug dumping of tables. +// #include +// #include +// #include + +namespace Opm +{ + + //------------------------------------------------------------------------ + // Member functions + //------------------------------------------------------------------------- + + /// Constructor + SinglePvtDead::SinglePvtDead(const table_t& pvd_table) + { + const int region_number = 0; + if (pvd_table.size() != 1) { + THROW("More than one PVT-region"); + } + + // Copy data + const int sz = pvd_table[region_number][0].size(); + std::vector press(sz); + std::vector B_inv(sz); + std::vector visc(sz); + for (int i = 0; i < sz; ++i) { + press[i] = pvd_table[region_number][0][i]; + B_inv[i] = 1.0 / pvd_table[region_number][1][i]; + visc[i] = pvd_table[region_number][2][i]; + } + one_over_B_ = NonuniformTableLinear(press, B_inv); + viscosity_ = NonuniformTableLinear(press, visc); + + // Dumping the created tables. +// static int count = 0; +// std::ofstream os((std::string("dump-") + boost::lexical_cast(count++)).c_str()); +// os.precision(15); +// os << "1/B\n\n" << one_over_B_ +// << "\n\nvisc\n\n" << viscosity_ << std::endl; + } + + // Destructor + SinglePvtDead::~SinglePvtDead() + { + } + + + + void SinglePvtDead::mu(const int n, + const double* p, + const double* /*z*/, + double* output_mu) const + { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + output_mu[i] = viscosity_(p[i]); + } + } + + void SinglePvtDead::B(const int n, + const double* p, + const double* /*z*/, + double* output_B) const + { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + output_B[i] = 1.0/one_over_B_(p[i]); + } + } + + void SinglePvtDead::dBdp(const int n, + const double* p, + const double* /*z*/, + double* output_B, + double* output_dBdp) const + { + B(n, p, 0, output_B); +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + double Bg = output_B[i]; + output_dBdp[i] = -Bg*Bg*one_over_B_.derivative(p[i]); + } + } + + + void SinglePvtDead::R(const int n, + const double* /*p*/, + const double* /*z*/, + double* output_R) const + { + std::fill(output_R, output_R + n, 0.0); + } + + void SinglePvtDead::dRdp(const int n, + const double* /*p*/, + const double* /*z*/, + double* output_R, + double* output_dRdp) const + { + std::fill(output_R, output_R + n, 0.0); + std::fill(output_dRdp, output_dRdp + n, 0.0); + } + +} diff --git a/opm/core/fluid/blackoil/SinglePvtDead.hpp b/opm/core/fluid/blackoil/SinglePvtDead.hpp new file mode 100644 index 00000000..629053fe --- /dev/null +++ b/opm/core/fluid/blackoil/SinglePvtDead.hpp @@ -0,0 +1,84 @@ +/* + 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_SINGLEPVTDEAD_HEADER_INCLUDED +#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED + + +#include +#include +#include + +namespace Opm +{ + + /// Class for immiscible dead oil and dry gas. + /// For all the virtual methods, the following apply: p and z + /// are expected to be of size n and n*num_phases, respectively. + /// Output arrays shall be of size n, and must be valid before + /// calling the method. + class SinglePvtDead : public SinglePvtInterface + { + public: + typedef std::vector > > table_t; + + SinglePvtDead(const table_t& pvd_table); + virtual ~SinglePvtDead(); + + /// Viscosity as a function of p and z. + virtual void mu(const int n, + const double* p, + const double* z, + double* output_mu) const; + + /// Formation volume factor as a function of p and z. + virtual void B(const int n, + const double* p, + const double* z, + double* output_B) const; + + /// Formation volume factor and p-derivative as functions of p and z. + virtual void dBdp(const int n, + const double* p, + const double* z, + double* output_B, + double* output_dBdp) const; + + /// Solution factor as a function of p and z. + virtual void R(const int n, + const double* p, + const double* z, + double* output_R) const; + + /// Solution factor and p-derivative as functions of p and z. + virtual void dRdp(const int n, + const double* p, + const double* z, + double* output_R, + double* output_dRdp) const; + private: + // PVT properties of dry gas or dead oil + NonuniformTableLinear one_over_B_; + NonuniformTableLinear viscosity_; + }; + +} + + +#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED diff --git a/opm/core/utility/NonuniformTableLinear.hpp b/opm/core/utility/NonuniformTableLinear.hpp new file mode 100644 index 00000000..71334460 --- /dev/null +++ b/opm/core/utility/NonuniformTableLinear.hpp @@ -0,0 +1,244 @@ +/* + Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010, 2011, 2012 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED +#define OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED + +#include +#include +#include +#include + +#include +#include + +namespace Opm +{ + + + /// @brief Exception used for domain errors. + struct OutsideDomainException : public std::exception {}; + + + /// @brief This class uses linear interpolation to compute the value + /// (and its derivative) of a function f sampled at possibly + /// nonuniform points. + /// @tparam T the range type of the function (should be an algebraic ring type) + template + class NonuniformTableLinear + { + public: + /// @brief Default constructor. + NonuniformTableLinear(); + + /// @brief Construct from vectors of x and y values. + /// @param x_values vector of domain values + /// @param y_values vector of corresponding range values. + NonuniformTableLinear(const std::vector& x_values, + const std::vector& y_values); + + /// @brief Get the domain. + /// @return the domain as a pair of doubles. + std::pair domain(); + + /// @brief Rescale the domain. + /// @param new_domain the new domain as a pair of doubles. + void rescaleDomain(std::pair new_domain); + + /// @brief Evaluate the value at x. + /// @param x a domain value + /// @return f(x) + double operator()(const double x) const; + + /// @brief Evaluate the derivative at x. + /// @param x a domain value + /// @return f'(x) + double derivative(const double x) const; + + /// @brief Evaluate the inverse at y. Requires T to be a double. + /// @param y a range value + /// @return f^{-1}(y) + double inverse(const double y) const; + + /// @brief Equality operator. + /// @param other another NonuniformTableLinear. + /// @return true if they are represented exactly alike. + bool operator==(const NonuniformTableLinear& other) const; + + /// @brief Policies for how to behave when trying to evaluate outside the domain. + enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2}; + + /// @brief Sets the behavioural policy for evaluation to the left of the domain. + /// @param rp the policy + void setLeftPolicy(RangePolicy rp); + + /// @brief Sets the behavioural policy for evaluation to the right of the domain. + /// @param rp the policy + void setRightPolicy(RangePolicy rp); + + protected: + std::vector x_values_; + std::vector y_values_; + mutable std::vector x_values_reversed_; + mutable std::vector y_values_reversed_; + RangePolicy left_; + RangePolicy right_; + }; + + + // A utility function + /// @brief Detect if a sequence is nondecreasing. + /// @tparam FI a forward iterator whose value type has operator< defined. + /// @param beg start of sequence + /// @param end one-beyond-end of sequence + /// @return false if there exists two consecutive values (v1, v2) in the sequence + /// for which v2 < v1, else returns true. + template + bool isNondecreasing(const FI beg, const FI end) + { + if (beg == end) return true; + FI it = beg; + ++it; + FI prev = beg; + for (; it != end; ++it, ++prev) { + if (*it < *prev) { + return false; + } + } + return true; + } + + + + // Member implementations. + + template + inline + NonuniformTableLinear + ::NonuniformTableLinear() + : left_(ClosestValue), right_(ClosestValue) + { + } + + template + inline + NonuniformTableLinear + ::NonuniformTableLinear(const std::vector& x_values, + const std::vector& y_values) + : x_values_(x_values), y_values_(y_values), + left_(ClosestValue), right_(ClosestValue) + { + ASSERT(isNondecreasing(x_values.begin(), x_values.end())); + } + + template + inline std::pair + NonuniformTableLinear + ::domain() + { + return std::make_pair(x_values_[0], x_values_.back()); + } + + template + inline void + NonuniformTableLinear + ::rescaleDomain(std::pair new_domain) + { + const double a = x_values_[0]; + const double b = x_values_.back(); + const double c = new_domain.first; + const double d = new_domain.second; + // x in [a, b] -> x in [c, d] + for (int i = 0; i < int(x_values_.size()); ++i) { + x_values_[i] = (x_values_[i] - a)*(d - c)/(b - a) + c; + } + } + + template + inline double + NonuniformTableLinear + ::operator()(const double x) const + { + return Opm::linearInterpolation(x_values_, y_values_, x); + } + + template + inline double + NonuniformTableLinear + ::derivative(const double x) const + { + return Opm::linearInterpolationDerivative(x_values_, y_values_, x); + } + + template + inline double + NonuniformTableLinear + ::inverse(const double y) const + { + if (y_values_.front() < y_values_.back()) { + return Opm::linearInterpolation(y_values_, x_values_, y); + } else { + if (y_values_reversed_.empty()) { + y_values_reversed_ = y_values_; + std::reverse(y_values_reversed_.begin(), y_values_reversed_.end()); + ASSERT(isNondecreasing(y_values_reversed_.begin(), y_values_reversed_.end())); + x_values_reversed_ = x_values_; + std::reverse(x_values_reversed_.begin(), x_values_reversed_.end()); + } + return Opm::linearInterpolation(y_values_reversed_, x_values_reversed_, y); + } + } + + template + inline bool + NonuniformTableLinear + ::operator==(const NonuniformTableLinear& other) const + { + return x_values_ == other.x_values_ + && y_values_ == other.y_values_ + && left_ == other.left_ + && right_ == other.right_; + } + + template + inline void + NonuniformTableLinear + ::setLeftPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + left_ = rp; + } + + template + inline void + NonuniformTableLinear + ::setRightPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + right_ = rp; + } + +} // namespace Opm + +#endif // OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED diff --git a/tests/bo_resprop_test.cpp b/tests/bo_resprop_test.cpp index 5bf3aea1..41baccbb 100644 --- a/tests/bo_resprop_test.cpp +++ b/tests/bo_resprop_test.cpp @@ -43,7 +43,8 @@ int main(int argc, char** argv) UnstructuredGrid grid; grid.number_of_cells = 1; grid.global_cell = NULL; - Opm::BlackoilPropertiesFromDeck props(deck, grid); + bool props_use_spline = param.getDefault("props_use_spline", true); + Opm::BlackoilPropertiesFromDeck props(deck, grid, props_use_spline); const int n = 1; double p[n] = { 150e5 };