Merge remote-tracking branch 'upstream/master'

opm/core/io/eclipse/EclipseGridParser.cpp

@@ -558,13 +558,15 @@ void EclipseGridParser::convertToSI()
                    key == "LAMEMOD"  || key == "SHEARMOD" || key == "POISSONMOD" ||
                    key == "PWAVEMOD" || key == "MULTPV"   || key == "PWAVEMOD" ||
                    key == "SGAS"     || key == "SWAT"     || key == "SOIL"     ||
-                   key == "RS"       || key == "SWCR"     || key == "SWL"      || 
+                                        key == "SWCR"     || key == "SWL"      ||
                    key == "SWU"      || key == "SOWCR"    || key == "KRW"      || 
                    key == "KRWR"     || key == "KRORW"    || key == "KRO") {
             unit = 1.0;
             do_convert = false; // Dimensionless keywords...
         } else if (key == "PRESSURE") {
             unit = units_.pressure;
+        } else if (key == "RS") {
+            unit = units_.gasvol_s / units_.liqvol_s;
         } else if (key == "MAPAXES") {
             MESSAGE("Not applying units to MAPAXES yet!");
             unit = 1.0;

opm/core/pressure/CompressibleTpfa.cpp

@@ -429,7 +429,7 @@ namespace Opm
             for (int phase = 0; phase < np; ++phase) {
                 int upwindc = -1;
                 if (c[0] >=0 && c[1] >= 0) {
-                    upwindc = (pot[0] < pot[1]) ? c[1] : c[0];
+                    upwindc = (pot[0][phase] < pot[1][phase]) ? c[1] : c[0];
                 } else {
                     upwindc = (c[0] >= 0) ? c[0] : c[1];
                 }

opm/core/props/pvt/SinglePvtConstCompr.hpp

@@ -31,8 +31,10 @@ namespace Opm
 {
 
     /// Class for constant compressible phases (PVTW or PVCDO).
-    /// For all the virtual methods, the following apply: p and z
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
+    /// or pressure (p) and gas resolution factor (r).
+    /// For all the virtual methods, the following apply: p, r and z
+    /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
     class SinglePvtConstCompr : public SinglePvtInterface
@@ -83,6 +85,29 @@ namespace Opm
             }
         }
 
+        virtual void mu(const int n,
+                        const double* p,
+                        const double* /*r*/,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const
+        {
+            if (visc_comp_) {
+                // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    // Computing a polynomial approximation to the exponential.
+                    double x = -visc_comp_*(p[i] - ref_press_);
+                    double d = (1.0 + x + 0.5*x*x);
+                    output_mu[i] = viscosity_/d;
+                    output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
+                }
+            } else {
+                std::fill(output_mu, output_mu + n, viscosity_);
+                std::fill(output_dmudp, output_dmudp + n, 0.0);
+            }
+            std::fill(output_dmudr, output_dmudr + n, 0.0);
+        }
+
         virtual void B(const int n,
                        const double* p,
                        const double* /*z*/,
@@ -120,6 +145,43 @@ namespace Opm
             }
         }
 
+        virtual void b(const int n,
+                               const double* p,
+                               const double* /*r*/,
+                               double* output_b,
+                               double* output_dbdp,
+                               double* output_dbdr) const
+                {
+                    if (comp_) {
+        // #pragma omp parallel for
+                        for (int i = 0; i < n; ++i) {
+                            // Computing a polynomial approximation to the exponential.
+                            double x = comp_*(p[i] - ref_press_);
+                            double d = (1.0 + x + 0.5*x*x);
+
+                            // b = 1/B = d/ref_B_B;
+                            output_b[i] = d/ref_B_;
+                            output_dbdp[i] = (1 + x) * comp_/ref_B_;
+                        }
+                    } else {
+
+                        std::fill(output_b, output_b + n, 1/ref_B_);
+                        std::fill(output_dbdp, output_dbdp + n, 0.0);
+                    }
+                    std::fill(output_dbdr, output_dbdr + n, 0.0);
+
+                }
+
+        virtual void rbub(const int n,
+                             const double* /*p*/,
+                             double* output_rbub,
+                             double* output_drbubdp) const
+
+              {
+                  std::fill(output_rbub, output_rbub + n, 0.0);
+                  std::fill(output_drbubdp, output_drbubdp + n, 0.0);
+              }
+
         virtual void R(const int n,
                        const double* /*p*/,
                        const double* /*z*/,

opm/core/props/pvt/SinglePvtDead.cpp

@@ -44,14 +44,14 @@ namespace Opm
         // 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> b(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];
+            b[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);
+        b_ = NonuniformTableLinear<double>(press, b);
         viscosity_ = NonuniformTableLinear<double>(press, visc);
 
         // Dumping the created tables.
@@ -80,14 +80,31 @@ namespace Opm
         }
     }
 
+    void SinglePvtDead::mu(const int n,
+                               const double* p,
+                               const double* /*r*/,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+        {
+    // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                output_mu[i] = viscosity_(p[i]);
+                output_dmudp[i] = viscosity_.derivative(p[i]);
+            }
+            std::fill(output_dmudr, output_dmudr + n, 0.0);
+
+        }
+
     void SinglePvtDead::B(const int n,
                           const double* p,
                           const double* /*z*/,
                           double* output_B) const
     {
 // #pragma omp parallel for
+        // B = 1/b
         for (int i = 0; i < n; ++i) {
-            output_B[i] = 1.0/one_over_B_(p[i]);
+            output_B[i] = 1.0/b_(p[i]);
         }
     }
 
@@ -101,10 +118,36 @@ namespace Opm
 // #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]);
+            output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
         }
     }
 
+    void SinglePvtDead::b(const int n,
+                              const double* p,
+                              const double* /*r*/,
+                              double* output_b,
+                              double* output_dbdp,
+                              double* output_dbdr) const
+
+        {
+    // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                output_b[i] = b_(p[i]);
+                output_dbdp[i] = b_.derivative(p[i]);
+
+            }
+            std::fill(output_dbdr, output_dbdr + n, 0.0);
+
+        }
+
+    void SinglePvtDead::rbub(const int n,
+                             const double* /*p*/,
+                             double* output_rbub,
+                             double* output_drbubdp) const
+    {
+        std::fill(output_rbub, output_rbub + n, 0.0);
+        std::fill(output_drbubdp, output_drbubdp + n, 0.0);
+    }
 
     void SinglePvtDead::R(const int n,
                           const double* /*p*/,

opm/core/props/pvt/SinglePvtDead.hpp

@@ -29,8 +29,10 @@ 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.
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
+    /// or pressure (p) and gas resolution factor (r).
+    /// For all the virtual methods, the following apply: p, r and z
+    /// are expected to be of size n, 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
@@ -46,6 +48,14 @@ namespace Opm
                         const double* z,
                         double* output_mu) const;
 
+        /// Viscosity and its derivatives as a function of p and r.
+        virtual void mu(const int n,
+                        const double* p,
+                        const double* r,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const;
+
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
                        const double* p,
@@ -59,6 +69,22 @@ namespace Opm
                           double* output_B,
                           double* output_dBdp) const;
 
+        /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+        virtual void b(const int n,
+                       const double* p,
+                       const double* r,
+                       double* output_b,
+                       double* output_dbdp,
+                       double* output_dbdr) const;
+
+
+
+        /// Gas resolution and its derivatives at bublepoint as a function of p.
+        virtual void rbub(const int n,
+                          const double* p,
+                          double* output_rbub,
+                          double* output_drbubdp) const;
+
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
                        const double* p,
@@ -73,7 +99,7 @@ namespace Opm
                           double* output_dRdp) const;
     private:
         // PVT properties of dry gas or dead oil
-        NonuniformTableLinear<double> one_over_B_;
+        NonuniformTableLinear<double> b_;
         NonuniformTableLinear<double> viscosity_;
     };
 

opm/core/props/pvt/SinglePvtDeadSpline.cpp

@@ -52,7 +52,7 @@ namespace Opm
             B_inv[i] = 1.0 / pvd_table[region_number][1][i];
             visc[i]  = pvd_table[region_number][2][i];
         }
-        buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
+        buildUniformMonotoneTable(press, B_inv, samples, b_);
         buildUniformMonotoneTable(press, visc, samples, viscosity_);
 
         // Dumping the created tables.
@@ -81,6 +81,23 @@ namespace Opm
         }
     }
 
+    void SinglePvtDeadSpline::mu(const int n,
+                               const double* p,
+                               const double* /*r*/,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+        {
+    // #pragma omp parallel for
+
+            for (int i = 0; i < n; ++i) {
+                output_mu[i] = viscosity_(p[i]);
+                output_dmudp[i] = viscosity_.derivative(p[i]);
+            }
+            std::fill(output_dmudr, output_dmudr + n, 0.0);
+
+        }
+
     void SinglePvtDeadSpline::B(const int n,
                                 const double* p,
                                 const double* /*z*/,
@@ -88,7 +105,7 @@ namespace Opm
     {
 // #pragma omp parallel for
         for (int i = 0; i < n; ++i) {
-            output_B[i] = 1.0/one_over_B_(p[i]);
+            output_B[i] = 1.0/b_(p[i]);
         }
     }
 
@@ -102,10 +119,36 @@ namespace Opm
 // #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]);
+            output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
         }
     }
 
+    void SinglePvtDeadSpline::b(const int n,
+                              const double* p,
+                              const double* /*r*/,
+                              double* output_b,
+                              double* output_dbdp,
+                              double* output_dbdr) const
+
+        {
+    // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                output_b[i] = b_(p[i]);
+                output_dbdp[i] = b_.derivative(p[i]);
+
+            }
+            std::fill(output_dbdr, output_dbdr + n, 0.0);
+
+        }
+
+    void SinglePvtDeadSpline::rbub(const int n,
+                             const double* /*p*/,
+                             double* output_rbub,
+                             double* output_drbubdp) const
+    {
+        std::fill(output_rbub, output_rbub + n, 0.0);
+        std::fill(output_drbubdp, output_drbubdp + n, 0.0);
+    }
 
     void SinglePvtDeadSpline::R(const int n,
                                 const double* /*p*/,

opm/core/props/pvt/SinglePvtDeadSpline.hpp

@@ -29,8 +29,10 @@ 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.
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
+    /// or pressure (p) and gas resolution factor (r).
+    /// For all the virtual methods, the following apply: p, r and z
+    /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
     class SinglePvtDeadSpline : public SinglePvtInterface
@@ -47,6 +49,14 @@ namespace Opm
                         const double* z,
                         double* output_mu) const;
 
+        /// Viscosity and its derivatives as a function of p and r.
+        virtual void mu(const int n,
+                        const double* p,
+                        const double* r,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const;
+
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
                        const double* p,
@@ -60,6 +70,20 @@ namespace Opm
                           double* output_B,
                           double* output_dBdp) const;
 
+        /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+        virtual void b(const int n,
+                       const double* p,
+                       const double* r,
+                       double* output_b,
+                       double* output_dbdp,
+                       double* output_dbdr) const;
+
+        /// Gas resolution and its derivatives at bublepoint as a function of p.
+        virtual void rbub(const int n,
+                          const double* p,
+                          double* output_rbub,
+                          double* output_drbubdp) const;
+
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
                        const double* p,
@@ -74,7 +98,7 @@ namespace Opm
                           double* output_dRdp) const;
     private:
         // PVT properties of dry gas or dead oil
-        UniformTableLinear<double> one_over_B_;
+        UniformTableLinear<double> b_;
         UniformTableLinear<double> viscosity_;
     };
 

opm/core/props/pvt/SinglePvtInterface.hpp

@@ -42,8 +42,10 @@ namespace Opm
         ///                          arbitrary two-phase and three-phase situations.
         void setPhaseConfiguration(const int num_phases, const int* phase_pos);
 
-        /// For all the virtual methods, the following apply: p and z
-        /// are expected to be of size n and n*num_phases, respectively.
+        /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
+        /// or pressure (p) and gas resolution factor (r).
+        /// For all the virtual methods, the following apply: p, r and z
+        /// are expected to be of size n, size n and n*num_phases, respectively.
         /// Output arrays shall be of size n, and must be valid before
         /// calling the method.
 
@@ -53,6 +55,14 @@ namespace Opm
                         const double* z,
                         double* output_mu) const = 0;
 
+        /// Viscosity as a function of p and r.
+        virtual void mu(const int n,
+                              const double* p,
+                              const double* r,
+                              double* output_mu,
+                              double* output_dmudp,
+                              double* output_dmudr) const = 0;
+
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
                        const double* p,
@@ -66,6 +76,21 @@ namespace Opm
                           double* output_B,
                           double* output_dBdp) const = 0;
 
+        /// The inverse of the volume factor b = 1 / B as a function of p and r.
+        virtual void b(const int n,
+                          const double* p,
+                          const double* r,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dpdr) const = 0;
+
+        /// Gas resolution at bublepoint as a function of pressure
+        virtual void rbub(const int n,
+                          const double* p,
+                          double* output_rbub,
+                          double* output_drbubdp) const = 0;
+
+
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
                        const double* p,
@@ -78,6 +103,8 @@ namespace Opm
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const = 0;
+
+
     protected:
         int num_phases_;
         int phase_pos_[MaxNumPhases];

opm/core/props/pvt/SinglePvtLiveGas.cpp

@@ -98,6 +98,17 @@ namespace Opm
         }
     }
 
+    /// Viscosity and its derivatives as a function of p and r.
+    void SinglePvtLiveGas::mu(const int n,
+                               const double* p,
+                               const double* r,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+    {
+        THROW("The new fluid interface not yet implemented");
+    }
+
 
     /// Formation volume factor as a function of p and z.
     void SinglePvtLiveGas::B(const int n,
@@ -126,6 +137,26 @@ namespace Opm
         }
     }
 
+    /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+    void SinglePvtLiveGas::b(const int n,
+                          const double* p,
+                          const double* r,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dbdr) const
+
+    {
+        THROW("The new fluid interface not yet implemented");
+    }
+
+    /// Gas resolution and its derivatives at bublepoint as a function of p.
+    void SinglePvtLiveGas::rbub(const int n,
+                             const double* p,
+                             double* output_rbub,
+                             double* output_drbubdp) const
+    {
+        THROW("The new fluid interface not yet implemented");
+    }
 
     /// Solution factor as a function of p and z.
     void SinglePvtLiveGas::R(const int n,

opm/core/props/pvt/SinglePvtLiveGas.hpp

@@ -26,8 +26,10 @@
 namespace Opm
 {
     /// Class for miscible wet gas (with vaporized oil in vapour phase).
-    /// For all the virtual methods, the following apply: p and z
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
+    /// or pressure (p) and gas resolution factor (r).
+    /// For all the virtual methods, the following apply: p, r and z
+    /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
     class SinglePvtLiveGas : public SinglePvtInterface
@@ -44,6 +46,14 @@ namespace Opm
                         const double* z,
                         double* output_mu) const;
 
+        /// Viscosity and its derivatives as a function of p and r.
+        virtual void mu(const int n,
+                        const double* p,
+                        const double* r,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const;
+
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
                        const double* p,
@@ -57,6 +67,22 @@ namespace Opm
                           double* output_B,
                           double* output_dBdp) const;
 
+        /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+        virtual void b(const int n,
+                       const double* p,
+                       const double* r,
+                       double* output_b,
+                       double* output_dbdp,
+                       double* output_dbdr) const;
+
+
+
+        /// Gas resolution and its derivatives at bublepoint as a function of p.
+        virtual void rbub(const int n,
+                          const double* p,
+                          double* output_rbub,
+                          double* output_drbubdp) const;
+
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
                        const double* p,

opm/core/props/pvt/SinglePvtLiveOil.cpp

@@ -175,6 +175,23 @@ namespace Opm
         }
     }
 
+    /// Viscosity and its derivatives as a function of p and r.
+    void SinglePvtLiveOil::mu(const int n,
+                               const double* p,
+                               const double* r,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    output_mu[i] = miscible_oil(p[i], r[i], 2, 0);
+                    output_dmudp[i] = miscible_oil(p[i], r[i], 2, 1);
+                    output_dmudr[i] = miscible_oil(p[i], r[i], 2, 2);
+
+                }
+    }
+
 
     /// Formation volume factor as a function of p and z.
     void SinglePvtLiveOil::B(const int n,
@@ -203,6 +220,37 @@ namespace Opm
         }
     }
 
+    void SinglePvtLiveOil::b(const int n,
+                          const double* p,
+                          const double* r,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dbdr) const
+
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    output_b[i] = miscible_oil(p[i], r[i], 1, 0);
+                    output_dbdp[i] = miscible_oil(p[i], r[i], 1, 1);
+                    output_dbdr[i] = miscible_oil(p[i], r[i], 1, 2);
+
+                }
+    }
+
+    void SinglePvtLiveOil::rbub(const int n,
+                             const double* p,
+                             double* output_rbub,
+                             double* output_drbubdp) const
+    {
+
+        for (int i = 0; i < n; ++i) {
+            output_rbub[i] = linearInterpolation(saturated_oil_table_[0],
+                    saturated_oil_table_[3],p[i]);
+            output_drbubdp[i] = linearInterpolationDerivative(saturated_oil_table_[0],
+                    saturated_oil_table_[3],p[i]);
+
+        }
+    }
 
     /// Solution factor as a function of p and z.
     void SinglePvtLiveOil::R(const int n,
@@ -347,4 +395,84 @@ namespace Opm
         }
     }
 
+    double SinglePvtLiveOil::miscible_oil(const double press,
+                                          const double r,
+                                          const int item,
+                                          const int deriv) const
+    {
+        int section;
+        double Rval = linearInterpolation(saturated_oil_table_[0],
+                                                saturated_oil_table_[3],
+                                                press, section);
+        // derivative with respect to frist component (pressure)
+        if (deriv == 1) {
+            if (Rval < r ) {  // Saturated case
+                return linearInterpolationDerivative(saturated_oil_table_[0],
+                                                saturated_oil_table_[item],
+                                                press);
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[3], r);
+                double w = (r - saturated_oil_table_[3][is]) /
+                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
+                ASSERT(undersat_oil_tables_[is][0].size() >= 2);
+                ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolationDerivative(undersat_oil_tables_[is][0],
+                                             undersat_oil_tables_[is][item],
+                                             press);
+                double val2 =
+                    linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
+                                             undersat_oil_tables_[is+1][item],
+                                             press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+            // derivative with respect to second component (r)
+        } else if (deriv == 2)  {
+            if (Rval < r ) {  // Saturated case
+                return 0;
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[3], r);
+                ASSERT(undersat_oil_tables_[is][0].size() >= 2);
+                ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[is][0],
+                                              undersat_oil_tables_[is][item],
+                                              press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[is+1][0],
+                                              undersat_oil_tables_[is+1][item],
+                                              press);
+
+                double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
+                return val;
+            }
+
+
+            }    else {
+            if (Rval < r ) {  // Saturated case
+                return linearInterpolation(saturated_oil_table_[0],
+                                                 saturated_oil_table_[item],
+                                                 press);
+            } else {  // Undersaturated case
+                // Interpolate between table sections
+                int is = tableIndex(saturated_oil_table_[3], r);
+                double w = (r - saturated_oil_table_[3][is]) /
+                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
+                ASSERT(undersat_oil_tables_[is][0].size() >= 2);
+                ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[is][0],
+                                              undersat_oil_tables_[is][item],
+                                              press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[is+1][0],
+                                              undersat_oil_tables_[is+1][item],
+                                              press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        }
+    }
+
 } // namespace Opm

opm/core/props/pvt/SinglePvtLiveOil.hpp

@@ -27,8 +27,10 @@
 namespace Opm
 {
     /// Class for miscible live oil (with dissolved gas in liquid phase).
-    /// For all the virtual methods, the following apply: p and z
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
+    /// or pressure (p) and gas resolution factor (r).
+    /// For all the virtual methods, the following apply: p, r and z
+    /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
     class SinglePvtLiveOil : public SinglePvtInterface
@@ -45,6 +47,14 @@ namespace Opm
                         const double* z,
                         double* output_mu) const;
 
+        /// Viscosity and its derivatives as a function of p and r.
+        virtual void mu(const int n,
+                        const double* p,
+                        const double* r,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const;
+
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
                        const double* p,
@@ -58,6 +68,20 @@ namespace Opm
                           double* output_B,
                           double* output_dBdp) const;
 
+        /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+        virtual void b(const int n,
+                       const double* p,
+                       const double* r,
+                       double* output_b,
+                       double* output_dbdp,
+                       double* output_dbdr) const;
+
+        /// Gas resolution and its derivatives at bublepoint as a function of p.
+        virtual void rbub(const int n,
+                          const double* p,
+                          double* output_rbub,
+                          double* output_drbubdp) const;
+
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
                        const double* p,
@@ -83,6 +107,11 @@ namespace Opm
                             const int item,
                             const bool deriv = false) const;
 
+        double miscible_oil(const double press,
+                            const double r,
+                            const int item,
+                            const int deriv = 0) const;
+
         // PVT properties of live oil (with dissolved gas)
         std::vector<std::vector<double> > saturated_oil_table_;
         std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;

opm/core/props/rock/RockCompressibility.cpp

@@ -42,12 +42,16 @@ namespace Opm
     {
         if (deck.hasField("ROCKTAB")) {
             const table_t& rt = deck.getROCKTAB().rocktab_;
-            int n = rt[0][0].size();
+            if (rt.size() != 1) {
+                THROW("Can only handle a single region in ROCKTAB.");
+            }
+            const int n = rt[0][0].size();
             p_.resize(n);
             poromult_.resize(n);
             for (int i = 0; i < n; ++i) {
                 p_[i] = rt[0][0][i];
                 poromult_[i] = rt[0][1][i];
+                transmult_[i] = rt[0][2][i];
             }
         } else if (deck.hasField("ROCK")) {
             const ROCK& r = deck.getROCK();
@@ -70,11 +74,39 @@ namespace Opm
             const double cpnorm = rock_comp_*(pressure - pref_);
             return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
         } else {
-            // return Opm::linearInterpolation(p_, poromult_, pressure);
             return Opm::linearInterpolation(p_, poromult_, pressure);
         }
     }
 
+    double RockCompressibility::poroMultDeriv(double pressure) const
+    {
+        if (p_.empty()) {
+            // Approximating poro multiplier with a quadratic curve,
+            // we must use its derivative.
+            return rock_comp_ + 2 * rock_comp_ * rock_comp_ * (pressure - pref_);
+        } else {
+            return Opm::linearInterpolationDerivative(p_, poromult_, pressure);
+        }
+    }
+
+    double RockCompressibility::transMult(double pressure) const
+    {
+        if (p_.empty()) {
+            return 1.0;
+        } else {
+            return Opm::linearInterpolation(p_, transmult_, pressure);
+        }
+    }
+
+    double RockCompressibility::transMultDeriv(double pressure) const
+    {
+        if (p_.empty()) {
+            return 0.0;
+        } else {
+            return Opm::linearInterpolationDerivative(p_, transmult_, pressure);
+        }
+    }
+
     double RockCompressibility::rockComp(double pressure) const
     {
         if (p_.empty()) {

opm/core/props/rock/RockCompressibility.hpp

@@ -47,12 +47,22 @@ namespace Opm
         /// Porosity multiplier.
         double poroMult(double pressure) const;
 
+        /// Derivative of porosity multiplier with respect to pressure.
+        double poroMultDeriv(double pressure) const;
+
+        /// Transmissibility multiplier.
+        double transMult(double pressure) const;
+
+        /// Derivative of transmissibility multiplier with respect to pressure.
+        double transMultDeriv(double pressure) const;
+
         /// Rock compressibility = (d poro / d p)*(1 / poro).
         double rockComp(double pressure) const;
 
     private:
         std::vector<double> p_;
         std::vector<double> poromult_;
+        std::vector<double> transmult_;
         double pref_;
         double rock_comp_;
     };

opm/core/props/satfunc/SatFuncGwseg.cpp

@@ -136,14 +136,14 @@ namespace Opm
             const double sg = s[Vapour];
             // xw and xg are the fractions occupied by water and gas zones.
             const double eps = 1e-6;
-            const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
-            const double xg = 1 - xw;
             const double ssw = sg + sw;
-            const double ssg = sw - swco + sg;
+            const double ssg = std::max(sg + sw - swco, eps);
             const double krw = krw_(ssw);
             const double krg = krg_(ssg);
             const double krow = krow_(ssw);
             const double krog = krog_(ssg);
+            const double xw = (sw - swco) / ssg;
+            const double xg = 1 - xw;
             kr[Aqua]   = xw*krw;
             kr[Vapour] = xg*krg;
             kr[Liquid] = xw*krow + xg*krog;
@@ -341,14 +341,14 @@ namespace Opm
             const double sg = s[Vapour];
             // xw and xg are the fractions occupied by water and gas zones.
             const double eps = 1e-6;
-            const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
-            const double xg = 1 - xw;
             const double ssw = sg + sw;
-            const double ssg = sw - swco + sg;
+            const double ssg = std::max(sg + sw - swco, eps);
             const double krw = krw_(ssw);
             const double krg = krg_(ssg);
             const double krow = krow_(ssw);
             const double krog = krog_(ssg);
+            const double xw = (sw - swco) / ssg;
+            const double xg = 1 - xw;
             kr[Aqua]   = xw*krw;
             kr[Vapour] = xg*krg;
             kr[Liquid] = xw*krow + xg*krog;

opm/core/simulator/BlackoilState.hpp

@@ -22,6 +22,7 @@
 
 #include <opm/core/grid.h>
 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
 #include <vector>
 
 namespace Opm
@@ -47,6 +48,7 @@ namespace Opm
                 // but use available phase information instead.
                 sat_[num_phases*cell + 1] = 1.0;
             }
+            gor_.resize(g.number_of_cells, 0.0);
         }
 
         enum ExtremalSat { MinSat, MaxSat };
@@ -59,7 +61,11 @@ namespace Opm
                          const Opm::BlackoilPropertiesInterface& props,
                          ExtremalSat es)
         {
+            if (cells.empty()) {
+                return;
+            }
             const int n = cells.size();
+            ASSERT(n > 0);
             std::vector<double> smin(num_phases_*n);
             std::vector<double> smax(num_phases_*n);
             props.satRange(n, &cells[0], &smin[0], &smax[0]);
@@ -81,12 +87,14 @@ namespace Opm
         std::vector<double>& faceflux    () { return flux_  ; }
         std::vector<double>& surfacevol  () { return surfvol_; }
         std::vector<double>& saturation  () { return sat_   ; }
+        std::vector<double>& gasoilratio () { return gor_   ; }
 
         const std::vector<double>& pressure    () const { return press_ ; }
         const std::vector<double>& facepressure() const { return fpress_; }
         const std::vector<double>& faceflux    () const { return flux_  ; }
         const std::vector<double>& surfacevol  () const { return surfvol_; }
         const std::vector<double>& saturation  () const { return sat_   ; }
+        const std::vector<double>& gasoilratio () const { return gor_   ; }
 
     private:
         int num_phases_;
@@ -95,6 +103,7 @@ namespace Opm
         std::vector<double> flux_  ;
         std::vector<double> surfvol_;
         std::vector<double> sat_   ;
+        std::vector<double> gor_   ;
     };
 
 } // namespace Opm

opm/core/simulator/WellState.hpp

@@ -31,27 +31,46 @@ namespace Opm
     {
     public:
         /// Allocate and initialize if wells is non-null.
+        /// Also tries to give useful initial values to the bhp() and
+        /// wellRates() fields, depending on controls.  The
+        /// perfRates() field is filled with zero, and perfPress()
+        /// with -1e100.
         template <class State>
         void init(const Wells* wells, const State& state)
         {
             if (wells) {
                 const int nw = wells->number_of_wells;
+                const int np = wells->number_of_phases;
                 bhp_.resize(nw);
-                // Initialize bhp to be target pressure
-                // if bhp-controlled well, otherwise set
-                // to pressure in first perforation cell.
+                wellrates_.resize(nw * np, 0.0);
                 for (int w = 0; w < nw; ++w) {
                     const WellControls* ctrl = wells->ctrls[w];
-
+                    // Initialize bhp to be target pressure if
+                    // bhp-controlled well, otherwise set to a little
+                    // above or below (depending on if the well is an
+                    // injector or producer) pressure in first perforation
+                    // cell.
                     if ((ctrl->current < 0) || // SHUT
                         (ctrl->type[ctrl->current] != BHP)) {
-                        const int cell = wells->well_cells[wells->well_connpos[w]];
-                        bhp_[w] = state.pressure()[cell];
-                    }
-                    else {
+                        const int first_cell = wells->well_cells[wells->well_connpos[w]];
+                        const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
+                        bhp_[w] = safety_factor*state.pressure()[first_cell];
+                    } else {
                         bhp_[w] = ctrl->target[ctrl->current];
                     }
+                    // Initialize well rates to match controls if type is SURFACE_RATE
+                    if ((ctrl->current >= 0) && // open well
+                        (ctrl->type[ctrl->current] == SURFACE_RATE)) {
+                        const double rate_target = ctrl->target[ctrl->current];
+                        for (int p = 0; p < np; ++p) {
+                            const double phase_distr = ctrl->distr[np * ctrl->current + p];
+                            wellrates_[np*w + p] = rate_target * phase_distr;
+                        }
+                    }
                 }
+                // The perforation rates and perforation pressures are
+                // not expected to be consistent with bhp_ and wellrates_
+                // after init().
                 perfrates_.resize(wells->well_connpos[nw], 0.0);
                 perfpress_.resize(wells->well_connpos[nw], -1e100);
             }
@@ -61,6 +80,10 @@ namespace Opm
         std::vector<double>& bhp() { return bhp_; }
         const std::vector<double>& bhp() const { return bhp_; }
 
+        /// One rate per well and phase.
+        std::vector<double>& wellRates() { return wellrates_; }
+        const std::vector<double>& wellRates() const { return wellrates_; }
+
         /// One rate per well connection.
         std::vector<double>& perfRates() { return perfrates_; }
         const std::vector<double>& perfRates() const { return perfrates_; }
@@ -71,6 +94,7 @@ namespace Opm
 
     private:
         std::vector<double> bhp_;
+        std::vector<double> wellrates_;
         std::vector<double> perfrates_;
         std::vector<double> perfpress_;
     };

opm/core/simulator/initState_impl.hpp

@@ -599,6 +599,29 @@ namespace Opm
     }
 
 
+    /// Initialize a blackoil state from input deck.
+    template <class Props, class State>
+    void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
+                                   const Props& props,
+                                   const EclipseGridParser& deck,
+                                   const double gravity,
+                                   State& state)
+    {
+        initStateFromDeck(grid, props, deck, gravity, state);
+        initBlackoilSurfvol(grid, props, state);
+        if (deck.hasField("RS")) {
+            const std::vector<double>& rs_deck = deck.getFloatingPointValue("RS");
+            const int num_cells = grid.number_of_cells;
+            for (int c = 0; c < num_cells; ++c) {
+                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                state.gasoilratio()[c] = rs_deck[c_deck];
+            }
+        } else {
+            THROW("Temporarily, we require the RS field.");
+        }
+    }
+
+
 } // namespace Opm
 
 

opm/core/wells/WellsManager.cpp

@@ -505,13 +505,33 @@ namespace Opm
                         int control_pos[5] = { -1, -1, -1, -1, -1 };
                         if (ok && wci_line.surface_flow_max_rate_ >= 0.0) {
                             control_pos[InjectionControl::RATE] = w_->ctrls[wix]->num;
-                            const double distr[3] = { 1.0, 1.0, 1.0 };
+                            double distr[3] = { 0.0, 0.0, 0.0 };
+                            if (wci_line.injector_type_ == "WATER") {
+                                distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                            } else if (wci_line.injector_type_ == "OIL") {
+                                distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                            } else if (wci_line.injector_type_ == "GAS") {
+                                distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                            } else {
+                                THROW("Injector type " << wci_line.injector_type_ << " not supported."
+                                      "WellsManager only supports WATER, OIL and GAS injector types.");
+                            }
                             ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_,
                                                       distr, wix, w_);
                         }
                         if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) {
                             control_pos[InjectionControl::RESV] = w_->ctrls[wix]->num;
-                            const double distr[3] = { 1.0, 1.0, 1.0 };
+                            double distr[3] = { 0.0, 0.0, 0.0 };
+                            if (wci_line.injector_type_ == "WATER") {
+                                distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                            } else if (wci_line.injector_type_ == "OIL") {
+                                distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                            } else if (wci_line.injector_type_ == "GAS") {
+                                distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                            } else {
+                                THROW("Injector type " << wci_line.injector_type_ << " not supported."
+                                      "WellsManager only supports WATER, OIL and GAS injector types.");
+                            }
                             ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_,
                                                       distr, wix, w_);
                         }