Merge remote-tracking branch 'totto82/newfluid2' into combined

CMakeLists_files.cmake

@@ -147,6 +147,7 @@ list (APPEND TEST_SOURCE_FILES
 	tests/test_column_extract.cpp
 	tests/test_geom2d.cpp
 	tests/test_param.cpp
+  tests/not-unit/test_newfluidinterface.cpp
 	)
 
 # originally generated with the command:

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
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// 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
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// 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
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// 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
+        /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-        /// are expected to be of size n and n*num_phases, respectively.
+        /// 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
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// 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
+    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// are expected to be of size n and n*num_phases, respectively.
+    /// 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_;

tests/not-unit/test_newfluidinterface.cpp

@@ -0,0 +1,255 @@
+#include <opm/core/io/eclipse/EclipseGridParser.hpp>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
+#include <opm/core/props/pvt/SinglePvtDead.hpp>
+#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
+#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
+#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
+#include <opm/core/props/phaseUsageFromDeck.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <iostream>
+#include <iterator>
+#include <vector>
+#include <string>
+using namespace Opm;
+using namespace std;
+
+
+// The function object divides a Factor with an element
+template <class Type>
+class MultValue
+{
+   private:
+      Type Factor;   // The value to multiply by
+   public:
+      // Constructor initializes the value to multiply by
+      MultValue ( const Type& _Val ) : Factor ( _Val ) {
+      }
+
+      // The function call for the element to be multiplied
+      int operator ( ) ( Type& elem ) const
+      {
+         return  Factor / elem;
+      }
+};
+int main () {
+    // read parameters from command-line
+    const string filename = "../../opm-core/tests/not-unit/blackoil/SPE9small.DATA";
+    cout << "Reading deck: " << filename << endl;
+    const EclipseGridParser deck (filename);
+    std::string mu_output = "mu_output";
+    std::string b_output = "b_output";
+    std::string rbub_output = "rbub_output";
+
+    PhaseUsage phase_usage_;
+
+    std::vector<std::tr1::shared_ptr<SinglePvtInterface> > props_;
+
+    phase_usage_ = phaseUsageFromDeck(deck);
+    enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
+
+    int samples = 0;
+
+    std::fstream muos(mu_output.c_str(), std::fstream::out | std::fstream::trunc);
+    if(!(muos.good())){
+        std::cout << "Could not open"<< mu_output << std::endl;
+        exit(3);
+    }
+
+    std::fstream bos(b_output.c_str(), std::fstream::out | std::fstream::trunc);
+    bos << setiosflags(ios::scientific) << setprecision(12);
+    if(!(bos.good())){
+        std::cout << "Could not open"<< b_output << std::endl;
+        exit(3);
+    }
+
+    std::fstream rbubos(rbub_output.c_str(), std::fstream::out | std::fstream::trunc);
+    rbubos << setiosflags(ios::scientific) << setprecision(12);
+    if(!(rbubos.good())){
+        std::cout << "Could not open"<< rbub_output << std::endl;
+        exit(3);
+    }
+
+
+
+    // Set the properties.
+        props_.resize(phase_usage_.num_phases);
+        // Water PVT
+        if (phase_usage_.phase_used[Aqua]) {
+            if (deck.hasField("PVTW")) {
+                props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(deck.getPVTW().pvtw_));
+            } else {
+                // Eclipse 100 default.
+                props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
+            }
+        }
+        // Oil PVT
+        if (phase_usage_.phase_used[Liquid]) {
+            if (deck.hasField("PVDO")) {
+                if (samples > 0) {
+                    props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_, samples));
+                } 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")) {
+                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(deck.getPVCDO().pvcdo_));
+            } else {
+                THROW("Input is missing PVDO or PVTO\n");
+            }
+        }
+        // Gas PVT
+        if (phase_usage_.phase_used[Vapour]) {
+            if (deck.hasField("PVDG")) {
+                if (samples > 0) {
+                    props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_, samples));
+                } 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 {
+                THROW("Input is missing PVDG or PVTG\n");
+            }
+        }
+
+
+        int n = 6;
+        int np = 3; //phase_usage_.num_phases;
+        double p[n];
+        double r[n];
+        double z[np*n];
+
+        double mu[n];
+        double dmudp[n];
+        double dmudr[n];
+        double mu_new[n];
+        double dmudp_diff;
+        double dmudr_diff;
+        double dmudp_diff_u;
+        double dmudr_diff_u;
+
+        //double rf[n];
+
+        double h = 1;
+        double rh = 1;
+
+        p[0] = 10000000;
+        p[1] = p[0] + h;
+        p[2] = 10000000;
+
+        p[3] = 10000000;
+        p[4] = p[0] + h;
+        p[5] = 10000000;
+
+
+        // saturated
+        r[0] = 200;
+        r[1] = 200;
+        r[2] = 200 + rh;
+
+        // undersaturated
+        r[3] = 50;
+        r[4] = 50;
+        r[5] = 50 +rh;
+
+
+        for (int i = 0; i < n; ++i) {
+                z[0+i*np] = 0; z[1+i*np] = 1;
+                z[2+i*np] = r[i];
+
+       }
+
+
+
+        // test mu
+        for (int phase = 1; phase < 2; ++phase) {
+            props_[phase]->mu(n, p, r, mu_new,dmudp,dmudr);
+            props_[phase]->mu(n, p, z, mu);
+            dmudp_diff = (mu_new[1]-mu_new[0])/h;
+            dmudr_diff = (mu_new[2]-mu_new[0])/rh;
+            dmudp_diff_u = (mu_new[4]-mu_new[3])/h;
+            dmudr_diff_u = (mu_new[5]-mu_new[3])/rh;
+
+            std::copy(mu,mu + n, std::ostream_iterator<double>(muos, " "));
+            muos << "\n";
+            std::copy(mu_new,mu_new + n, std::ostream_iterator<double>(muos, " "));
+            muos << "\n";
+            std::copy(dmudp,dmudp + n, std::ostream_iterator<double>(muos, " "));
+            muos << "\n";
+            muos << dmudp_diff << " " << dmudp_diff_u << "\n";
+            std::copy(dmudr,dmudr + n, std::ostream_iterator<double>(muos, " "));
+            muos << "\n";
+            muos << dmudr_diff << " " << dmudr_diff_u << "\n";
+        }
+
+        // test b
+        double b[n];
+        double B[n];
+        double invB[n];
+        double dinvBdp[n];
+        double dBdp[n];
+        double dbdr[n];
+        double dbdp[n];
+        double dbdp_diff;
+        double dbdr_diff;
+        double dbdp_diff_u;
+        double dbdr_diff_u;
+
+        for (int phase = 1; phase < 2; ++phase) {
+            props_[phase]->b(n, p, r, b,dbdp,dbdr);
+            //props_[phase]->B(n, p, z, B);
+            props_[phase]->dBdp(n, p, z, B, dBdp);
+            dbdp_diff = (b[1]-b[0])/h;
+            dbdr_diff = (b[2]-b[0])/rh;
+            dbdp_diff_u = (b[4]-b[3])/h;
+            dbdr_diff_u = (b[5]-b[3])/rh;
+            for (int i = 0; i < n; ++i){
+                invB[i] = 1/B[i];
+                dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
+
+            }
+            std::copy(b,b + n, std::ostream_iterator<double>(bos, " "));
+            bos << "\n";
+            std::copy(invB,invB + n, std::ostream_iterator<double>(bos, " "));
+            bos << "\n";
+            std::copy(dinvBdp,dinvBdp + n, std::ostream_iterator<double>(bos, " "));
+            bos << "\n";
+            std::copy(dbdp,dbdp + n, std::ostream_iterator<double>(bos, " "));
+            bos << "\n";
+            bos << dbdp_diff << " " << dbdp_diff_u << "\n";
+            std::copy(dbdr,dbdr + n, std::ostream_iterator<double>(bos, " "));
+            bos << "\n";
+            bos << dbdr_diff << " " << dbdr_diff_u << "\n";
+
+        }
+
+        // test rbub
+
+        double rbub[n];
+        double drbubdp[n];
+        double drbubdp_diff;
+
+        for (int phase = 1; phase < 2; ++phase) {
+            props_[phase] ->rbub(n,p,rbub,drbubdp);
+
+            drbubdp_diff = (rbub[1]-rbub[0])/h;
+            std::copy(rbub,rbub + n, std::ostream_iterator<double>(rbubos, " "));
+            rbubos << "\n";
+            std::copy(drbubdp,drbubdp + n, std::ostream_iterator<double>(rbubos, " "));
+            rbubos << drbubdp_diff;
+            rbubos << "\n";
+
+
+        }
+
+
+
+
+}
+

tutorials/tutorial1.cpp

@@ -1,5 +1,5 @@
+/// \cond SKIP
 /*!
-\cond SKIP
   Copyright 2012 SINTEF ICT, Applied Mathematics.
 
   This file is part of the Open Porous Media project (OPM).
@@ -16,11 +16,8 @@
 
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-\endcond
 */
-#if HAVE_CONFIG_H
+///\endcond
-#include "config.h"
-#endif // HAVE_CONFIG_H
 
 /// \page tutorial1 A simple cartesian grid
 /// This tutorial explains how to construct a simple Cartesian grid,
@@ -34,6 +31,8 @@
 /// \snippet tutorial1.cpp including headers
 
 /// \internal [including headers] 
+#include "config.h"
+
 #include <opm/core/grid.h>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/core/io/vtk/writeVtkData.hpp>

tutorials/tutorial2.cpp

@@ -1,5 +1,5 @@
+/// \cond SKIP
 /*!
-\cond SKIP
   Copyright 2012 SINTEF ICT, Applied Mathematics.
 
   This file is part of the Open Porous Media project (OPM).
@@ -16,8 +16,9 @@
 
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-\endcond
 */
+/// \endcond
+
 /// \page tutorial2 Flow Solver for a single phase
 /// \details The flow equations consist of the mass conservation equation
 /// \f[\nabla\cdot {\bf u}=q\f] and the Darcy law \f[{\bf u} =- \frac{1}{\mu}K\nabla p.\f] Here,
@@ -27,9 +28,7 @@
 /// We solve the flow equations for a Cartesian grid and we set the source term
 /// \f$q\f$ be zero except at the left-lower and right-upper corner, where it is equal
 /// with opposite sign (inflow equal to outflow).
-#if HAVE_CONFIG_H
 #include "config.h"
-#endif // HAVE_CONFIG_H
 
 #include <opm/core/grid.h>
 #include <opm/core/grid/GridManager.hpp>

tutorials/tutorial3.cpp

@@ -1,5 +1,5 @@
+/// \cond SKIP
 /*!
-\cond SKIP
   Copyright 2012 SINTEF ICT, Applied Mathematics.
 
   This file is part of the Open Porous Media project (OPM).
@@ -16,11 +16,9 @@
 
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-\endcond
 */
-#if HAVE_CONFIG_H
+/// \endcond
 #include "config.h"
-#endif // HAVE_CONFIG_H
 
 #include <iostream>
 #include <iomanip>

tutorials/tutorial4.cpp

@@ -1,5 +1,5 @@
+/// \cond SKIP
 /*!
-\cond SKIP
   Copyright 2012 SINTEF ICT, Applied Mathematics.
 
   This file is part of the Open Porous Media project (OPM).
@@ -16,11 +16,9 @@
 
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-\endcond
 */
-#if HAVE_CONFIG_H
+/// \endcond
 #include "config.h"
-#endif // HAVE_CONFIG_H
 
 #include <iostream>
 #include <iomanip>