Merge branch 'master' into reorder_tof

opm/core/fluid/BlackoilPropertiesFromDeck.cpp

@@ -18,6 +18,7 @@
 */
 
 #include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 
 namespace Opm
 {
@@ -26,11 +27,59 @@ namespace Opm
                                                            const UnstructuredGrid& grid)
     {
         rock_.init(deck, grid);
-        pvt_.init(deck);
+        pvt_.init(deck, 200);
-        satprops_.init(deck, grid);
+        SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
-	if (pvt_.numPhases() != satprops_.numPhases()) {
+            = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
-	    THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
+        satprops_.reset(ptr);
-		  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
+        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,
+                                                           const UnstructuredGrid& grid,
+                                                           const parameter::ParameterGroup& param)
+    {
+        rock_.init(deck, grid);
+        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<std::string>("threephase_model", "simple");
+        bool use_stone2 = (threephase_model == "stone2");
+        if (sat_samples > 1) {
+            if (use_stone2) {
+                SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, grid, sat_samples);
+            } else {
+                SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, grid, sat_samples);
+            }
+        } else {
+            if (use_stone2) {
+                SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, grid, sat_samples);
+            } else {
+                SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, grid, sat_samples);
+            }
+        }
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
         }
     }
 
@@ -235,7 +284,7 @@ namespace Opm
                                              double* kr,
                                              double* dkrds) const
     {
-        satprops_.relperm(n, s, cells, kr, dkrds);
+        satprops_->relperm(n, s, cells, kr, dkrds);
     }
 
 
@@ -254,7 +303,7 @@ namespace Opm
                                               double* pc,
                                               double* dpcds) const
     {
-        satprops_.capPress(n, s, cells, pc, dpcds);
+        satprops_->capPress(n, s, cells, pc, dpcds);
     }
 
 
@@ -270,7 +319,7 @@ namespace Opm
                                               double* smin,
                                               double* smax) const
     {
-	satprops_.satRange(n, cells, smin, smax);
+        satprops_->satRange(n, cells, smin, smax);
     }
 
 

opm/core/fluid/BlackoilPropertiesFromDeck.hpp

@@ -26,6 +26,8 @@
 #include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
 #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <boost/scoped_ptr.hpp>
 
 struct UnstructuredGrid;
 
@@ -38,13 +40,28 @@ namespace Opm
     {
     public:
         /// Initialize from deck and grid.
-        /// \param  deck         Deck input parser
+        /// \param[in]  deck     Deck input parser
-        /// \param  grid         Grid to which property object applies, needed for the 
+        /// \param[in]  grid     Grid to which property object applies, needed for the
         ///                      mapping from cell indices (typically from a processed grid)
         ///                      to logical cartesian indices consistent with the deck.
         BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
                                    const UnstructuredGrid& grid);
 
+        /// Initialize from deck, grid and parameters.
+        /// \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.
+        /// \param[in]  param    Parameters. Accepted parameters include:
+        ///                        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);
+
         /// Destructor.
         virtual ~BlackoilPropertiesFromDeck();
 
@@ -165,7 +182,7 @@ namespace Opm
     private:
         RockFromDeck rock_;
         BlackoilPvtProperties pvt_;
-        SaturationPropsFromDeck satprops_;
+        boost::scoped_ptr<SaturationPropsInterface> satprops_;
         mutable std::vector<double> B_;
         mutable std::vector<double> dB_;
         mutable std::vector<double> R_;

opm/core/fluid/IncompPropertiesFromDeck.cpp

@@ -31,7 +31,7 @@ namespace Opm
     {
         rock_.init(deck, grid);
         pvt_.init(deck);
-        satprops_.init(deck, grid);
+        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() << ").");

opm/core/fluid/IncompPropertiesFromDeck.hpp

@@ -135,7 +135,7 @@ namespace Opm
     private:
         RockFromDeck rock_;
         PvtPropertiesIncompFromDeck pvt_;
-        SaturationPropsFromDeck satprops_;
+        SaturationPropsFromDeck<SatFuncStone2Uniform> satprops_;
     };
 
 

opm/core/fluid/SaturationPropsFromDeck.hpp

@@ -20,9 +20,12 @@
 #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
 #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
 
+#include <opm/core/fluid/SaturationPropsInterface.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
-#include <opm/core/utility/UniformTableLinear.hpp>
 #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+#include <opm/core/fluid/SatFuncStone2.hpp>
+#include <opm/core/fluid/SatFuncSimple.hpp>
 #include <vector>
 
 struct UnstructuredGrid;
@@ -30,19 +33,31 @@ struct UnstructuredGrid;
 namespace Opm
 {
 
-    class SaturationPropsFromDeck : public BlackoilPhases
+
+
+    /// Interface to saturation functions from deck.
+    /// Possible values for template argument (for now):
+    ///   SatFuncSetStone2Nonuniform,
+    ///   SatFuncSetStone2Uniform.
+    ///   SatFuncSetSimpleNonuniform,
+    ///   SatFuncSetSimpleUniform.
+    template <class SatFuncSet>
+    class SaturationPropsFromDeck : public SaturationPropsInterface
     {
     public:
         /// Default constructor.
         SaturationPropsFromDeck();
 
         /// Initialize from deck and grid.
-        /// \param  deck         Deck input parser
+        /// \param[in]  deck     Deck input parser
-        /// \param  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]  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 UnstructuredGrid& grid,
+                  const int samples);
 
         /// \return   P, the number of phases.
         int numPhases() const;
@@ -88,30 +103,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<double> krw_;
-            UniformTableLinear<double> krow_;
-            UniformTableLinear<double> pcow_;
-            UniformTableLinear<double> krg_;
-            UniformTableLinear<double> krog_;
-            UniformTableLinear<double> pcog_;
-            double krocw_; // = krow_(s_wc)
-        };
         std::vector<SatFuncSet> satfuncset_;
         std::vector<int> cell_to_func_; // = SATNUM - 1
 
-        const SatFuncSet& funcForCell(const int cell) const;
+        typedef SatFuncSet Funcs;
+
+        const Funcs& funcForCell(const int cell) const;
     };
 
 
@@ -119,6 +116,7 @@ namespace Opm
 } // namespace Opm
 
 
+#include <opm/core/fluid/SaturationPropsFromDeck_impl.hpp>
 
 
 #endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
+#define OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
+
+
+#include <opm/core/utility/UniformTableLinear.hpp>
+#include <opm/core/utility/NonuniformTableLinear.hpp>
+#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
+#include <opm/core/grid.h>
+
+namespace Opm
+{
+
+
+    // ----------- Methods of SaturationPropsFromDeck ---------
+
+
+    /// Default constructor.
+    template <class SatFuncSet>
+    SaturationPropsFromDeck<SatFuncSet>::SaturationPropsFromDeck()
+    {
+    }
+
+    /// Initialize from deck.
+    template <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::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<int>& 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 <class SatFuncSet>
+    int SaturationPropsFromDeck<SatFuncSet>::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 <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::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 <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::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 <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::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 <class SatFuncSet>
+    const typename SaturationPropsFromDeck<SatFuncSet>::Funcs&
+    SaturationPropsFromDeck<SatFuncSet>::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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
+#define OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
+
+#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+
+
+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

opm/core/fluid/blackoil/BlackoilPvtProperties.hpp

@@ -47,7 +47,13 @@ namespace Opm
         BlackoilPvtProperties();
 
         /// Initialize from deck.
-	void init(const EclipseGridParser& deck);
+        /// \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;

opm/core/pressure/CompressibleTpfa.hpp

@@ -90,14 +90,13 @@ namespace Opm
         bool singularPressure() const;
 
     private:
-        void computePerSolveDynamicData(const double dt,
+        virtual void computePerSolveDynamicData(const double dt,
                                                 const BlackoilState& state,
                                                 const WellState& well_state);
-        void computeWellPotentials(const BlackoilState& state);
         void computePerIterationDynamicData(const double dt,
                                             const BlackoilState& state,
                                             const WellState& well_state);
-        void computeCellDynamicData(const double dt,
+        virtual void computeCellDynamicData(const double dt,
                                             const BlackoilState& state,
                                             const WellState& well_state);
         void computeFaceDynamicData(const double dt,
@@ -114,6 +113,8 @@ namespace Opm
         double incrementNorm() const;
         void computeResults(BlackoilState& state,
                             WellState& well_state) const;
+    protected:
+        void computeWellPotentials(const BlackoilState& state);
 
         // ------ Data that will remain unmodified after construction. ------
         const UnstructuredGrid& grid_;

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -124,7 +124,7 @@ namespace Opm
     //
     // [[ incompressible was: r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) ]]
     //
-    //     r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s) )
+    //     r(s) = s - B*z0 + s*(poro - poro0)/poro0 + dt/pv*( influx + outflux*f(s) )
     //
     // @@@ What about the source term
     //

opm/core/utility/initState_impl.hpp

@@ -512,11 +512,11 @@ namespace Opm
                            State& state)
     {
         const int num_phases = props.numPhases();
-        if (num_phases != 2) {
-            THROW("initStateFromDeck(): currently handling only two-phase scenarios.");
-        }
         state.init(grid, num_phases);
         if (deck.hasField("EQUIL")) {
+            if (num_phases != 2) {
+                THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
+            }
             // Set saturations depending on oil-water contact.
             const EQUIL& equil= deck.getEQUIL();
             if (equil.equil.size() != 1) {
@@ -535,12 +535,28 @@ namespace Opm
             const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
             const std::vector<double>& 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];
                 }
+            } else if (num_phases == 3) {
+                if (!deck.hasField("SGAS")) {
+                    THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found).");
+                }
+                const std::vector<double>& 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];
+                    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.");
         }

opm/core/utility/miscUtilities.cpp

@@ -70,6 +70,25 @@ namespace Opm
         }
     }
 
+    /// @brief Computes porosity of all cells in a grid, with rock compressibility effects.
+    /// @param[in]  grid               a grid
+    /// @param[in]  porosity_standard  array of grid.number_of_cells porosity values (at standard conditions)
+    /// @param[in]  rock_comp          rock compressibility properties
+    /// @param[in]  pressure           pressure by cell
+    /// @param[out] porosity           porosity (at reservoir condition)
+    void computePorosity(const UnstructuredGrid& grid,
+                         const double* porosity_standard,
+                         const RockCompressibility& rock_comp,
+                         const std::vector<double>& pressure,
+                         std::vector<double>& porosity)
+    {
+        int num_cells = grid.number_of_cells;
+        porosity.resize(num_cells);
+        for (int i = 0; i < num_cells; ++i) {
+            porosity[i] = porosity_standard[i]*rock_comp.poroMult(pressure[i]);
+        }
+    }
+
 
     /// @brief Computes total saturated volumes over all grid cells.
     /// @param[in]  pv        the pore volume by cell.
@@ -574,8 +593,10 @@ namespace Opm
     {
         int nw = well_bhp.size();
         ASSERT(nw == wells.number_of_wells);
-        if (props.numPhases() != 2) {
+        int np = props.numPhases();
-            THROW("WellReport for now assumes two phase flow.");
+        const int max_np = 3;
+        if (np > max_np) {
+            THROW("WellReport for now assumes #phases <= " << max_np);
         }
         const double* visc = props.viscosity();
         std::vector<double> data_now;
@@ -586,7 +607,8 @@ namespace Opm
             double well_rate_total = 0.0;
             double well_rate_water = 0.0;
             for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) {
-                const double perf_rate = well_perfrates[perf]*(unit::day/unit::second);
+                const double perf_rate = unit::convert::to(well_perfrates[perf],
+                                                           unit::cubic(unit::meter)/unit::day);
                 well_rate_total += perf_rate;
                 if (perf_rate > 0.0) {
                     // Injection.
@@ -594,11 +616,14 @@ namespace Opm
                 } else {
                     // Production.
                     const int cell = wells.well_cells[perf];
-                    double mob[2];
+                    double mob[max_np];
                     props.relperm(1, &saturation[2*cell], &cell, mob, 0);
-                    mob[0] /= visc[0];
+                    double tmob = 0;
-                    mob[1] /= visc[1];
+                    for(int i = 0; i < np; ++i) {
-                    const double fracflow = mob[0]/(mob[0] + mob[1]);
+                        mob[i] /= visc[i];
+                        tmob += mob[i];
+                    }
+                    const double fracflow = mob[0]/tmob;
                     well_rate_water += perf_rate*fracflow;
                 }
             }
@@ -627,8 +652,10 @@ namespace Opm
         // TODO: refactor, since this is almost identical to the other push().
         int nw = well_bhp.size();
         ASSERT(nw == wells.number_of_wells);
-        if (props.numPhases() != 2) {
+        int np = props.numPhases();
-            THROW("WellReport for now assumes two phase flow.");
+        const int max_np = 3;
+        if (np > max_np) {
+            THROW("WellReport for now assumes #phases <= " << max_np);
         }
         std::vector<double> data_now;
         data_now.reserve(1 + 3*nw);
@@ -638,7 +665,8 @@ namespace Opm
             double well_rate_total = 0.0;
             double well_rate_water = 0.0;
             for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) {
-                const double perf_rate = well_perfrates[perf]*(unit::day/unit::second);
+                const double perf_rate = unit::convert::to(well_perfrates[perf],
+                                                           unit::cubic(unit::meter)/unit::day);
                 well_rate_total += perf_rate;
                 if (perf_rate > 0.0) {
                     // Injection.
@@ -646,13 +674,16 @@ namespace Opm
                 } else {
                     // Production.
                     const int cell = wells.well_cells[perf];
-                    double mob[2];
+                    double mob[max_np];
-                    props.relperm(1, &s[2*cell], &cell, mob, 0);
+                    props.relperm(1, &s[np*cell], &cell, mob, 0);
-                    double visc[2];
+                    double visc[max_np];
-                    props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0);
+                    props.viscosity(1, &p[cell], &z[np*cell], &cell, visc, 0);
-                    mob[0] /= visc[0];
+                    double tmob = 0;
-                    mob[1] /= visc[1];
+                    for(int i = 0; i < np; ++i) {
-                    const double fracflow = mob[0]/(mob[0] + mob[1]);
+                        mob[i] /= visc[i];
+                        tmob += mob[i];
+                    }
+                    const double fracflow = mob[0]/(tmob);
                     well_rate_water += perf_rate*fracflow;
                 }
             }

opm/core/utility/miscUtilities.hpp

@@ -44,7 +44,7 @@ namespace Opm
 
     /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
     /// @param[in]  grid      a grid
-    /// @param[in]  porosity  array of grid.number_of_cells porosity values
+    /// @param[in]  porosity  array of grid.number_of_cells porosity values (at reference pressure)
     /// @param[in]  rock_comp rock compressibility properties
     /// @param[in]  pressure  pressure by cell
     /// @param[out] porevol   the pore volume by cell.
@@ -54,6 +54,17 @@ namespace Opm
                            const std::vector<double>& pressure,
                            std::vector<double>& porevol);
 
+    /// @brief Computes porosity of all cells in a grid, with rock compressibility effects.
+    /// @param[in]  grid               a grid
+    /// @param[in]  porosity_standard  array of grid.number_of_cells porosity values (at reference presure)
+    /// @param[in]  rock_comp          rock compressibility properties
+    /// @param[in]  pressure           pressure by cell
+    /// @param[out] porosity           porosity (at reservoir condition)
+    void computePorosity(const UnstructuredGrid& grid,
+                         const double* porosity_standard,
+                         const RockCompressibility& rock_comp,
+                         const std::vector<double>& pressure,
+                         std::vector<double>& porosity);
 
     /// @brief Computes total saturated volumes over all grid cells.
     /// @param[in]  pv        the pore volume by cell.