Whether the fluid is saturated or not is explicitly passed to the pvts

The criteria for whether the fluid is saturated or not is moved from the
within the pvt calculations to the solver, and passed to the pvt
calculations as a array of boolean values.

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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* /*isSat*/,
                              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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* isSat,
                                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(), isSat,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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* /*isSat*/,
                             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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* /*isSat*/,
                               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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                 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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                   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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                  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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* isSat,
                                      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(),isSat,
                                                 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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* isSat,
                                        const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -300,7 +304,7 @@ namespace Opm
         V dmudp(n);
         V dmudr(n);
 
-        props_[phase_usage_.phase_pos[Oil]]->mu(n, po.value().data(), rs.value().data(),
+        props_[phase_usage_.phase_pos[Oil]]->mu(n, po.value().data(), rs.value().data(), isSat,
                                                 mu.data(), dmudp.data(), dmudr.data());
 
         ADB::M dmudp_diag = spdiag(dmudp);
@@ -391,6 +395,7 @@ namespace Opm
     /// \return            Array of n formation volume factor values.
     V BlackoilPropsAdFromDeck::bOil(const V& po,
                                     const V& rs,
+                                    const bool* isSat,
                                     const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -403,7 +408,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(),isSat,
                                                b.data(), dbdp.data(), dbdr.data());
 
         return b;
@@ -466,10 +471,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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
     /// \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 bool* isSat,
                                       const Cells& cells) const
     {
         if (!phase_usage_.phase_used[Oil]) {
@@ -482,7 +489,7 @@ namespace Opm
         V dbdp(n);
         V dbdr(n);
 
-        props_[phase_usage_.phase_pos[Oil]]->b(n, po.value().data(), rs.value().data(),
+        props_[phase_usage_.phase_pos[Oil]]->b(n, po.value().data(), rs.value().data(),isSat,
                                                b.data(), dbdp.data(), dbdr.data());
 
         ADB::M dbdp_diag = spdiag(dbdp);

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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                 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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                   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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                  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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                 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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                   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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                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]  isSat  Array of n booleans telling whether the fluid is saturated or not.
         /// \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 bool* isSat,
                  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;
 
+        bool isSat[rs.size()];
+        getSaturatedCells(state,&isSat[0]);
+
         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, &isSat[0], cells_);
                 rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
                 // DUMP(rq_[pos].b);
                 // DUMP(rq_[pos].accum[aix]);
@@ -594,6 +595,8 @@ namespace {
             // DUMP(residual_.mass_balance[phase]);
         }
 
+        bool isSat[grid_.number_of_cells];
+        getSaturatedCells(state,isSat);
         // -------- Extra (optional) sg or rs equation, and rs contributions to the mass balance equations --------
 
         // Add the extra (flux) terms to the gas mass balance equations
@@ -613,8 +616,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);
         }
 
@@ -654,7 +658,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, &isSat[0],cells_);
                 cell_rho_total += state.saturation[pos] * cell_rho;
             }
         }
@@ -664,7 +668,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, &isSat[0],cells_);
                 const V fraction = compi.col(pos);
                 inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells);
             }
@@ -885,6 +889,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 +1007,14 @@ namespace {
                                                  const SolutionState&    state )
     {
         const int phase = canph_[ actph ];
-        const ADB mu    = fluidViscosity(phase, state.pressure, state.rs, cells_);
+        bool isSat[grid_.number_of_cells];
+        getSaturatedCells(state,&isSat[0]);
+        const ADB mu    = fluidViscosity(phase, state.pressure, state.rs, &isSat[0],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, &isSat[0],cells_);
 
         ADB& head = rq_[ actph ].head;
 
@@ -1060,13 +1067,14 @@ namespace {
     FullyImplicitBlackoilSolver::fluidViscosity(const int               phase,
                                                 const ADB&              p    ,
                                                 const ADB&              rs   ,
+                                                const bool*             isSat,
                                                 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, isSat,cells);
         }
         case Gas:
             return fluid_.muGas(p, cells);
@@ -1083,13 +1091,14 @@ namespace {
     FullyImplicitBlackoilSolver::fluidReciprocFVF(const int               phase,
                                                   const ADB&              p    ,
                                                   const ADB&              rs   ,
+                                                  const bool*             isSat,
                                                   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, isSat, cells);
         }
         case Gas:
             return fluid_.bGas(p, cells);
@@ -1106,10 +1115,11 @@ namespace {
     FullyImplicitBlackoilSolver::fluidDensity(const int               phase,
                                               const ADB&              p    ,
                                               const ADB&              rs   ,
+                                              const bool*             isSat,
                                               const std::vector<int>& cells) const
     {
         const double* rhos = fluid_.surfaceDensity();
-        ADB b = fluidReciprocFVF(phase, p, rs, cells);
+        ADB b = fluidReciprocFVF(phase, p, rs, isSat, 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 +1205,20 @@ namespace {
     }
 
 
+    void
+    FullyImplicitBlackoilSolver::getSaturatedCells(const SolutionState& state, bool* isSat) const
+    {
+        const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
+        const V sg = state.saturation[pg].value();
+        for (int c=0; c < sg.size(); ++ c) {
+            if (sg[c]>0){
+                isSat[c] = true;
+            }
+            else{
+                isSat[c] = false;
+            }
+        }
+    }
+
+
 } // namespace Opm

opm/autodiff/FullyImplicitBlackoilSolver.hpp

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

opm/autodiff/ImpesTPFAAD.cpp

@@ -533,7 +533,8 @@ namespace {
             return fluid_.muWat(p, cells);
         case Oil: {
             V dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.muOil(p, dummy_rs, cells);
+            bool dummy_isSat[p.size()];
+            return fluid_.muOil(p, dummy_rs, dummy_isSat, cells);
         }
         case Gas:
             return fluid_.muGas(p, cells);
@@ -553,7 +554,8 @@ namespace {
             return fluid_.muWat(p, cells);
         case Oil: {
             ADB dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.muOil(p, dummy_rs, cells);
+            bool dummy_isSat[p.size()];
+            return fluid_.muOil(p, dummy_rs, dummy_isSat, cells);
         }
         case Gas:
             return fluid_.muGas(p, cells);
@@ -573,7 +575,8 @@ namespace {
             return fluid_.bWat(p, cells);
         case Oil: {
             V dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.bOil(p, dummy_rs, cells);
+            bool dummy_isSat[p.size()];
+            return fluid_.bOil(p, dummy_rs, dummy_isSat,cells);
         }
         case Gas:
             return fluid_.bGas(p, cells);
@@ -593,7 +596,8 @@ namespace {
             return fluid_.bWat(p, cells);
         case Oil: {
             ADB dummy_rs = V::Zero(p.size(), 1) * p;
-            return fluid_.bOil(p, dummy_rs, cells);
+            bool dummy_isSat[p.size()];
+            return fluid_.bOil(p, dummy_rs, dummy_isSat,cells);
         }
         case Gas:
             return fluid_.bGas(p, cells);