Merge pull request #69 from bska/match-reference/properties

Switch PVT Evaluators to Use Explicit External Condition Interface

opm/autodiff/BlackoilPropsAd.cpp

@@ -121,10 +121,12 @@ namespace Opm
     /// Oil viscosity.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n viscosity values.
     V BlackoilPropsAd::muOil(const V& po,
                              const V& rs,
+                             const std::vector<PhasePresence>& /*cond*/,
                              const Cells& cells) const
     {
         if (!pu_.phase_used[Oil]) {
@@ -197,14 +199,16 @@ namespace Opm
     /// Oil viscosity.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n viscosity values.
     ADB BlackoilPropsAd::muOil(const ADB& po,
                                const ADB& rs,
+                               const std::vector<PhasePresence>& cond,
                                const Cells& cells) const
     {
 #if 1
-        return ADB::constant(muOil(po.value(), rs.value(), cells), po.blockPattern());
+        return ADB::constant(muOil(po.value(), rs.value(), cond, cells), po.blockPattern());
 #else
         if (!pu_.phase_used[Oil]) {
             OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
@@ -306,10 +310,12 @@ namespace Opm
     /// Oil formation volume factor.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n formation volume factor values.
     V BlackoilPropsAd::bOil(const V& po,
                             const V& rs,
+                            const std::vector<PhasePresence>& /*cond*/,
                             const Cells& cells) const
     {
         if (!pu_.phase_used[Oil]) {
@@ -382,10 +388,12 @@ namespace Opm
     /// Oil formation volume factor.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n formation volume factor values.
     ADB BlackoilPropsAd::bOil(const ADB& po,
                               const ADB& rs,
+                              const std::vector<PhasePresence>& /*cond*/,
                               const Cells& cells) const
     {
         if (!pu_.phase_used[Oil]) {

opm/autodiff/BlackoilPropsAd.hpp

@@ -109,10 +109,12 @@ namespace Opm
         /// Oil viscosity.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n viscosity values.
         V muOil(const V& po,
                 const V& rs,
+                const std::vector<PhasePresence>& cond,
                 const Cells& cells) const;
 
         /// Gas viscosity.
@@ -132,10 +134,12 @@ namespace Opm
         /// Oil viscosity.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n viscosity values.
         ADB muOil(const ADB& po,
                   const ADB& rs,
+                  const std::vector<PhasePresence>& cond,
                   const Cells& cells) const;
 
         /// Gas viscosity.
@@ -158,10 +162,12 @@ namespace Opm
         /// Oil formation volume factor.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n formation volume factor values.
         V bOil(const V& po,
                const V& rs,
+               const std::vector<PhasePresence>& cond,
                const Cells& cells) const;
 
         /// Gas formation volume factor.
@@ -181,10 +187,12 @@ namespace Opm
         /// Oil formation volume factor.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n formation volume factor values.
         ADB bOil(const ADB& po,
                  const ADB& rs,
+                 const std::vector<PhasePresence>& cond,
                  const Cells& cells) const;
 
         /// Gas formation volume factor.

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -212,10 +212,12 @@ namespace Opm
     /// Oil viscosity.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n viscosity values.
     V BlackoilPropsAdFromDeck::muOil(const V& po,
                                      const V& rs,
+                                     const std::vector<PhasePresence>& cond,
                                      const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -227,7 +229,7 @@ namespace Opm
         V dmudp(n);
         V dmudr(n);
 
-        props_[phase_usage_.phase_pos[Oil]]->mu(n, po.data(), rs.data(),
+        props_[phase_usage_.phase_pos[Oil]]->mu(n, po.data(), rs.data(), &cond[0],
                                                 mu.data(), dmudp.data(), dmudr.data());
         return mu;
     }
@@ -285,10 +287,12 @@ namespace Opm
     /// Oil viscosity.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n viscosity values.
     ADB BlackoilPropsAdFromDeck::muOil(const ADB& po,
                                        const ADB& rs,
+                                       const std::vector<PhasePresence>& cond,
                                        const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -301,7 +305,7 @@ namespace Opm
         V dmudr(n);
 
         props_[phase_usage_.phase_pos[Oil]]->mu(n, po.value().data(), rs.value().data(),
-                                                mu.data(), dmudp.data(), dmudr.data());
+                                                &cond[0], mu.data(), dmudp.data(), dmudr.data());
 
         ADB::M dmudp_diag = spdiag(dmudp);
         ADB::M dmudr_diag = spdiag(dmudr);
@@ -387,10 +391,12 @@ namespace Opm
     /// Oil formation volume factor.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n formation volume factor values.
     V BlackoilPropsAdFromDeck::bOil(const V& po,
                                     const V& rs,
+                                    const std::vector<PhasePresence>& cond,
                                     const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -403,7 +409,7 @@ namespace Opm
         V dbdp(n);
         V dbdr(n);
 
-        props_[phase_usage_.phase_pos[Oil]]->b(n, po.data(), rs.data(),
+        props_[phase_usage_.phase_pos[Oil]]->b(n, po.data(), rs.data(), &cond[0],
                                                b.data(), dbdp.data(), dbdr.data());
 
         return b;
@@ -466,10 +472,12 @@ namespace Opm
     /// Oil formation volume factor.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  rs     Array of n gas solution factor values.
+    /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
     /// \return            Array of n formation volume factor values.
     ADB BlackoilPropsAdFromDeck::bOil(const ADB& po,
                                       const ADB& rs,
+                                      const std::vector<PhasePresence>& cond,
                                       const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -483,7 +491,7 @@ namespace Opm
         V dbdr(n);
 
         props_[phase_usage_.phase_pos[Oil]]->b(n, po.value().data(), rs.value().data(),
-                                               b.data(), dbdp.data(), dbdr.data());
+                                               &cond[0], b.data(), dbdp.data(), dbdr.data());
 
         ADB::M dbdp_diag = spdiag(dbdp);
         ADB::M dbdr_diag = spdiag(dbdr);

opm/autodiff/BlackoilPropsAdFromDeck.hpp

@@ -110,10 +110,12 @@ namespace Opm
         /// Oil viscosity.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n viscosity values.
         V muOil(const V& po,
                 const V& rs,
+                const std::vector<PhasePresence>& cond,
                 const Cells& cells) const;
 
         /// Gas viscosity.
@@ -133,10 +135,12 @@ namespace Opm
         /// Oil viscosity.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n viscosity values.
         ADB muOil(const ADB& po,
                   const ADB& rs,
+                  const std::vector<PhasePresence>& cond,
                   const Cells& cells) const;
 
         /// Gas viscosity.
@@ -159,10 +163,12 @@ namespace Opm
         /// Oil formation volume factor.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n formation volume factor values.
         V bOil(const V& po,
                const V& rs,
+               const std::vector<PhasePresence>& cond,
                const Cells& cells) const;
 
         /// Gas formation volume factor.
@@ -182,10 +188,12 @@ namespace Opm
         /// Oil formation volume factor.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n formation volume factor values.
         ADB bOil(const ADB& po,
                  const ADB& rs,
+                 const std::vector<PhasePresence>& cond,
                  const Cells& cells) const;
 
         /// Gas formation volume factor.

opm/autodiff/BlackoilPropsAdInterface.hpp

@@ -100,11 +100,13 @@ namespace Opm
         /// Oil viscosity.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n viscosity values.
         virtual
         V muOil(const V& po,
                 const V& rs,
+                const std::vector<PhasePresence>& cond,
                 const Cells& cells) const = 0;
 
         /// Gas viscosity.
@@ -126,11 +128,13 @@ namespace Opm
         /// Oil viscosity.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n viscosity values.
         virtual
         ADB muOil(const ADB& po,
                   const ADB& rs,
+                  const std::vector<PhasePresence>& cond,
                   const Cells& cells) const = 0;
 
         /// Gas viscosity.
@@ -155,11 +159,13 @@ namespace Opm
         /// Oil formation volume factor.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n formation volume factor values.
         virtual
         V bOil(const V& po,
                const V& rs,
+               const std::vector<PhasePresence>& cond,
                const Cells& cells) const = 0;
 
         /// Gas formation volume factor.
@@ -181,11 +187,13 @@ namespace Opm
         /// Oil formation volume factor.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  rs     Array of n gas solution factor values.
+        /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
         /// \return            Array of n formation volume factor values.
         virtual
         ADB bOil(const ADB& po,
                  const ADB& rs,
+                 const std::vector<PhasePresence>& cond,
                  const Cells& cells) const = 0;
 
         /// Gas formation volume factor.

opm/autodiff/FullyImplicitBlackoilSolver.cpp

@@ -71,8 +71,6 @@ namespace {
         return all_cells;
     }
 
-
-
     template <class GeoProps>
     AutoDiffBlock<double>::M
     gravityOperator(const UnstructuredGrid& grid,
@@ -530,13 +528,16 @@ namespace {
         const std::vector<ADB>& sat   = state.saturation;
         const ADB&              rs    = state.rs;
 
+        std::vector<PhasePresence> cond;
+        classifyCondition(state, cond);
+
         const ADB pv_mult = poroMult(press);
 
         const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
         for (int phase = 0; phase < maxnp; ++phase) {
             if (active_[ phase ]) {
                 const int pos = pu.phase_pos[ phase ];
-                rq_[pos].b = fluidReciprocFVF(phase, press, rs, cells_);
+                rq_[pos].b = fluidReciprocFVF(phase, press, rs, cond, cells_);
                 rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
                 // DUMP(rq_[pos].b);
                 // DUMP(rq_[pos].accum[aix]);
@@ -613,8 +614,9 @@ namespace {
             const ADB sg_eq = state.saturation[pg];
             const ADB rs_max = fluidRsMax(state.pressure, cells_);
             const ADB rs_eq = state.rs - rs_max;
-            Selector<double> use_rs_eq(rs_eq.value());
+            // Consider the fluid to be saturated if sg >= 1e-14 (a small number)
-            residual_.rs_or_sg_eq = use_rs_eq.select(rs_eq, sg_eq);
+            Selector<double> use_sat_eq(sg_eq.value()-1e-14);
+            residual_.rs_or_sg_eq = use_sat_eq.select(rs_eq, sg_eq);
             // DUMP(residual_.rs_or_sg_eq);
         }
 
@@ -650,11 +652,15 @@ namespace {
                 assert(g[dd] == 0.0);
             }
         }
+
+        std::vector<PhasePresence> cond;
+        classifyCondition(state, cond);
+
         ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern());
         for (int phase = 0; phase < 3; ++phase) {
             if (active_[phase]) {
                 const int pos = pu.phase_pos[phase];
-                const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, cells_);
+                const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, cond, cells_);
                 cell_rho_total += state.saturation[pos] * cell_rho;
             }
         }
@@ -664,7 +670,7 @@ namespace {
         for (int phase = 0; phase < 3; ++phase) {
             if (active_[phase]) {
                 const int pos = pu.phase_pos[phase];
-                const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, cells_);
+                const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, cond, cells_);
                 const V fraction = compi.col(pos);
                 inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells);
             }
@@ -885,6 +891,7 @@ namespace {
                 // rs(sg > 0) = rs_sat(sg > 0);
                 // rs(rs > rs_sat*rs_adjust) = rs_sat(rs > rs_sat*rs_adjust);
                 for (int c = 0; c < nc; ++c) {
+
                     if (zerosg[c]) {
                         sg[c] = 0.0;
                     }
@@ -1002,12 +1009,16 @@ namespace {
                                                  const SolutionState&    state )
     {
         const int phase = canph_[ actph ];
-        const ADB mu    = fluidViscosity(phase, state.pressure, state.rs, cells_);
+
+        std::vector<PhasePresence> cond;
+        classifyCondition(state, cond);
+
+        const ADB mu    = fluidViscosity(phase, state.pressure, state.rs, cond, cells_);
         const ADB tr_mult = transMult(state.pressure);
 
         rq_[ actph ].mob = tr_mult * kr[ phase ] / mu;
 
-        const ADB rho   = fluidDensity(phase, state.pressure, state.rs, cells_);
+        const ADB rho   = fluidDensity(phase, state.pressure, state.rs, cond, cells_);
 
         ADB& head = rq_[ actph ].head;
 
@@ -1060,13 +1071,14 @@ namespace {
     FullyImplicitBlackoilSolver::fluidViscosity(const int               phase,
                                                 const ADB&              p    ,
                                                 const ADB&              rs   ,
+                                                const std::vector<PhasePresence>& cond,
                                                 const std::vector<int>& cells) const
     {
         switch (phase) {
         case Water:
             return fluid_.muWat(p, cells);
         case Oil: {
-            return fluid_.muOil(p, rs, cells);
+            return fluid_.muOil(p, rs, cond, cells);
         }
         case Gas:
             return fluid_.muGas(p, cells);
@@ -1083,13 +1095,14 @@ namespace {
     FullyImplicitBlackoilSolver::fluidReciprocFVF(const int               phase,
                                                   const ADB&              p    ,
                                                   const ADB&              rs   ,
+                                                  const std::vector<PhasePresence>& cond,
                                                   const std::vector<int>& cells) const
     {
         switch (phase) {
         case Water:
             return fluid_.bWat(p, cells);
         case Oil: {
-            return fluid_.bOil(p, rs, cells);
+            return fluid_.bOil(p, rs, cond, cells);
         }
         case Gas:
             return fluid_.bGas(p, cells);
@@ -1106,10 +1119,11 @@ namespace {
     FullyImplicitBlackoilSolver::fluidDensity(const int               phase,
                                               const ADB&              p    ,
                                               const ADB&              rs   ,
+                                              const std::vector<PhasePresence>& cond,
                                               const std::vector<int>& cells) const
     {
         const double* rhos = fluid_.surfaceDensity();
-        ADB b = fluidReciprocFVF(phase, p, rs, cells);
+        ADB b = fluidReciprocFVF(phase, p, rs, cond, cells);
         ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
         if (phase == Oil && active_[Gas]) {
             // It is correct to index into rhos with canonical phase indices.
@@ -1195,4 +1209,45 @@ namespace {
     }
 
 
+    void
+    FullyImplicitBlackoilSolver::
+    classifyCondition(const SolutionState&        state,
+                      std::vector<PhasePresence>& cond ) const
+    {
+        const PhaseUsage& pu = fluid_.phaseUsage();
+
+        if (active_[ Gas ]) {
+            // Oil/Gas or Water/Oil/Gas system
+            const int po = pu.phase_pos[ Oil ];
+            const int pg = pu.phase_pos[ Gas ];
+
+            const V&  so = state.saturation[ po ].value();
+            const V&  sg = state.saturation[ pg ].value();
+
+            cond.resize(sg.size());
+
+            for (V::Index c = 0, e = sg.size(); c != e; ++c) {
+                if (so[c] > 0)        { cond[c].setFreeOil  (); }
+                if (sg[c] > 0)        { cond[c].setFreeGas  (); }
+                if (active_[ Water ]) { cond[c].setFreeWater(); }
+            }
+        }
+        else {
+            // Water/Oil system
+            assert (active_[ Water ]);
+
+            const int po = pu.phase_pos[ Oil ];
+            const V&  so = state.saturation[ po ].value();
+
+            cond.resize(so.size());
+
+            for (V::Index c = 0, e = so.size(); c != e; ++c) {
+                cond[c].setFreeWater();
+
+                if (so[c] > 0) { cond[c].setFreeOil(); }
+            }
+        }
+    }
+
+
 } // namespace Opm

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -189,18 +189,21 @@ namespace Opm {
         fluidViscosity(const int               phase,
                        const ADB&              p    ,
                        const ADB&              rs   ,
+                       const std::vector<PhasePresence>& cond,
                        const std::vector<int>& cells) const;
 
         ADB
         fluidReciprocFVF(const int               phase,
                          const ADB&              p    ,
                          const ADB&              rs   ,
+                         const std::vector<PhasePresence>& cond,
                          const std::vector<int>& cells) const;
 
         ADB
         fluidDensity(const int               phase,
                      const ADB&              p    ,
                      const ADB&              rs   ,
+                     const std::vector<PhasePresence>& cond,
                      const std::vector<int>& cells) const;
 
         V
@@ -216,6 +219,10 @@ namespace Opm {
 
         ADB
         transMult(const ADB& p) const;
+
+        void
+        classifyCondition(const SolutionState&        state,
+                          std::vector<PhasePresence>& cond ) const;
     };
 } // namespace Opm
 

opm/autodiff/ImpesTPFAAD.cpp

@@ -533,7 +533,9 @@ namespace {
             return fluid_.muWat(p, cells);
         case Oil: {
             V dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.muOil(p, dummy_rs, cells);
+            std::vector<PhasePresence> cond(dummy_rs.size());
+
+            return fluid_.muOil(p, dummy_rs, cond, cells);
         }
         case Gas:
             return fluid_.muGas(p, cells);
@@ -553,7 +555,9 @@ namespace {
             return fluid_.muWat(p, cells);
         case Oil: {
             ADB dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.muOil(p, dummy_rs, cells);
+            std::vector<PhasePresence> cond(dummy_rs.size());
+
+            return fluid_.muOil(p, dummy_rs, cond, cells);
         }
         case Gas:
             return fluid_.muGas(p, cells);
@@ -573,7 +577,9 @@ namespace {
             return fluid_.bWat(p, cells);
         case Oil: {
             V dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.bOil(p, dummy_rs, cells);
+            std::vector<PhasePresence> cond(dummy_rs.size());
+
+            return fluid_.bOil(p, dummy_rs, cond, cells);
         }
         case Gas:
             return fluid_.bGas(p, cells);
@@ -593,7 +599,9 @@ namespace {
             return fluid_.bWat(p, cells);
         case Oil: {
             ADB dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.bOil(p, dummy_rs, cells);
+            std::vector<PhasePresence> cond(dummy_rs.size());
+
+            return fluid_.bOil(p, dummy_rs, cond, cells);
         }
         case Gas:
             return fluid_.bGas(p, cells);