Added class SinglePvtDead, add parameter 'props_use_spline' to simulators.

Recall that the class that used to be called SinglePvtDead has
been renamed to SinglePvtDeadSpline. If 'props_use_spline' is true,
that class is used (this is the default), which makes a monotone
spline that is uniformly, densely sampled. The new class simply
uses linear interpolation in the input tables.

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			\

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.

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);

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 ("

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();

opm/core/fluid/blackoil/BlackoilPvtProperties.cpp

@@ -20,6 +20,7 @@
 
 
 #include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
+#include <opm/core/fluid/blackoil/SinglePvtDead.hpp>
 #include <opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp>
 #include <opm/core/fluid/blackoil/SinglePvtLiveOil.hpp>
 #include <opm/core/fluid/blackoil/SinglePvtLiveGas.hpp>
@@ -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<std::vector<std::vector<double> > > 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 {

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;

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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <opm/core/fluid/blackoil/SinglePvtDead.hpp>
+#include <algorithm>
+
+// Extra includes for debug dumping of tables.
+// #include <boost/lexical_cast.hpp>
+// #include <string>
+// #include <fstream>
+
+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<double> press(sz);
+        std::vector<double> B_inv(sz);
+        std::vector<double> visc(sz);
+	for (int i = 0; i < sz; ++i) {
+            press[i] = pvd_table[region_number][0][i];
+            B_inv[i] = 1.0 / pvd_table[region_number][1][i];
+            visc[i]  = pvd_table[region_number][2][i];
+	}
+        one_over_B_ = NonuniformTableLinear<double>(press, B_inv);
+        viscosity_ = NonuniformTableLinear<double>(press, visc);
+
+        // Dumping the created tables.
+//         static int count = 0;
+//         std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
+//         os.precision(15);
+//         os << "1/B\n\n" << one_over_B_
+//            << "\n\nvisc\n\n" << viscosity_ << std::endl;
+    }
+
+    // Destructor
+    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);
+    }
+
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SINGLEPVTDEAD_HEADER_INCLUDED
+#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED
+
+
+#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
+#include <opm/core/utility/NonuniformTableLinear.hpp>
+#include <vector>
+
+namespace Opm
+{
+
+    /// Class for immiscible dead oil and dry gas.
+    /// For all the virtual methods, the following apply: p and z
+    /// are expected to be of size n and n*num_phases, respectively.
+    /// Output arrays shall be of size n, and must be valid before
+    /// calling the method.
+    class SinglePvtDead : public SinglePvtInterface
+    {
+    public:
+	typedef std::vector<std::vector<std::vector<double> > > 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<double> one_over_B_;
+        NonuniformTableLinear<double> viscosity_;
+    };
+
+}
+
+
+#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
+#define OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
+
+#include <cmath>
+#include <exception>
+#include <vector>
+#include <utility>
+
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/linearInterpolation.hpp>
+
+namespace Opm
+{
+
+
+    /// @brief Exception used for domain errors.
+    struct OutsideDomainException : public std::exception {};
+
+
+    /// @brief This class uses linear interpolation to compute the value
+    ///        (and its derivative) of a function f sampled at possibly
+    ///         nonuniform points.
+    /// @tparam T the range type of the function (should be an algebraic ring type)
+    template<typename T>
+    class NonuniformTableLinear
+    {
+    public:
+        /// @brief Default constructor.
+        NonuniformTableLinear();
+
+        /// @brief Construct from vectors of x and y values.
+        /// @param x_values vector of domain values
+        /// @param y_values vector of corresponding range values.
+        NonuniformTableLinear(const std::vector<double>& x_values,
+                              const std::vector<T>& y_values);
+
+        /// @brief Get the domain.
+        /// @return the domain as a pair of doubles.
+        std::pair<double, double> domain();
+
+        /// @brief Rescale the domain.
+        /// @param new_domain the new domain as a pair of doubles.
+        void rescaleDomain(std::pair<double, double> new_domain);
+
+        /// @brief Evaluate the value at x.
+        /// @param x a domain value
+        /// @return f(x)
+        double operator()(const double x) const;
+
+        /// @brief Evaluate the derivative at x.
+        /// @param x a domain value
+        /// @return f'(x)
+        double derivative(const double x) const;
+
+        /// @brief Evaluate the inverse at y. Requires T to be a double.
+        /// @param y a range value
+        /// @return f^{-1}(y)
+        double inverse(const double y) const;
+
+        /// @brief Equality operator.
+        /// @param other another NonuniformTableLinear.
+        /// @return true if they are represented exactly alike.
+        bool operator==(const NonuniformTableLinear& other) const;
+
+        /// @brief Policies for how to behave when trying to evaluate outside the domain.
+        enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2};
+
+        /// @brief Sets the behavioural policy for evaluation to the left of the domain.
+        /// @param rp the policy
+        void setLeftPolicy(RangePolicy rp);
+
+        /// @brief Sets the behavioural policy for evaluation to the right of the domain.
+        /// @param rp the policy
+        void setRightPolicy(RangePolicy rp);
+
+    protected:
+        std::vector<double> x_values_;
+        std::vector<T> y_values_;
+        mutable std::vector<T> x_values_reversed_;
+        mutable std::vector<T> y_values_reversed_;
+        RangePolicy left_;
+        RangePolicy right_;
+    };
+
+
+    // A utility function
+    /// @brief Detect if a sequence is nondecreasing.
+    /// @tparam FI a forward iterator whose value type has operator< defined.
+    /// @param beg start of sequence
+    /// @param end one-beyond-end of sequence
+    /// @return false if there exists two consecutive values (v1, v2) in the sequence
+    ///         for which v2 < v1, else returns true.
+    template <typename FI>
+    bool isNondecreasing(const FI beg, const FI end)
+    {
+        if (beg == end) return true;
+        FI it = beg;
+        ++it;
+        FI prev = beg;
+        for (; it != end; ++it, ++prev) {
+            if (*it < *prev) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+
+
+    // Member implementations.
+
+    template<typename T>
+    inline
+    NonuniformTableLinear<T>
+    ::NonuniformTableLinear()
+        : left_(ClosestValue), right_(ClosestValue)
+    {
+    }
+
+    template<typename T>
+    inline
+    NonuniformTableLinear<T>
+    ::NonuniformTableLinear(const std::vector<double>& x_values,
+                            const std::vector<T>& y_values)
+        : x_values_(x_values), y_values_(y_values),
+          left_(ClosestValue), right_(ClosestValue)
+    {
+        ASSERT(isNondecreasing(x_values.begin(), x_values.end()));
+    }
+
+    template<typename T>
+    inline std::pair<double, double>
+    NonuniformTableLinear<T>
+    ::domain()
+    {
+        return std::make_pair(x_values_[0], x_values_.back());
+    }
+
+    template<typename T>
+    inline void
+    NonuniformTableLinear<T>
+    ::rescaleDomain(std::pair<double, double> new_domain)
+    {
+        const double a = x_values_[0];
+        const double b = x_values_.back();
+        const double c = new_domain.first;
+        const double d = new_domain.second;
+        // x in [a, b] -> x in [c, d]
+        for (int i = 0; i < int(x_values_.size()); ++i) {
+            x_values_[i] = (x_values_[i] - a)*(d - c)/(b - a) + c;
+        }
+    }
+
+    template<typename T>
+    inline double
+    NonuniformTableLinear<T>
+    ::operator()(const double x) const
+    {
+        return Opm::linearInterpolation(x_values_, y_values_, x);
+    }
+
+    template<typename T>
+    inline double
+    NonuniformTableLinear<T>
+    ::derivative(const double x) const
+    {
+        return Opm::linearInterpolationDerivative(x_values_, y_values_, x);
+    }
+
+    template<typename T>
+    inline double
+    NonuniformTableLinear<T>
+    ::inverse(const double y) const
+    {
+        if (y_values_.front() < y_values_.back()) {
+            return Opm::linearInterpolation(y_values_, x_values_, y);
+        } else {
+            if (y_values_reversed_.empty()) {
+                y_values_reversed_ = y_values_;
+                std::reverse(y_values_reversed_.begin(), y_values_reversed_.end());
+                ASSERT(isNondecreasing(y_values_reversed_.begin(), y_values_reversed_.end()));
+                x_values_reversed_ = x_values_;
+                std::reverse(x_values_reversed_.begin(), x_values_reversed_.end());
+            }
+            return Opm::linearInterpolation(y_values_reversed_, x_values_reversed_, y);
+        }
+    }
+
+    template<typename T>
+    inline bool
+    NonuniformTableLinear<T>
+    ::operator==(const NonuniformTableLinear<T>& other) const
+    {
+        return x_values_ == other.x_values_
+            && y_values_ == other.y_values_
+            && left_ == other.left_
+            && right_ == other.right_;
+    }
+
+    template<typename T>
+    inline void
+    NonuniformTableLinear<T>
+    ::setLeftPolicy(RangePolicy rp)
+    {
+        if (rp != ClosestValue) {
+            THROW("Only ClosestValue RangePolicy implemented.");
+        }
+        left_ = rp;
+    }
+
+    template<typename T>
+    inline void
+    NonuniformTableLinear<T>
+    ::setRightPolicy(RangePolicy rp)
+    {
+        if (rp != ClosestValue) {
+            THROW("Only ClosestValue RangePolicy implemented.");
+        }
+        right_ = rp;
+    }
+
+} // namespace Opm
+
+#endif // OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED

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 };