From 7d1b59d69743c6a49317b40301147814ee19192a Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Wed, 28 Sep 2016 16:38:40 +0200
Subject: [PATCH] set the Rv and Rs factors to the saturated values for cells
 which have no gas and no oil

the purpose of this is to get a more defined behaviour when doing the
gravity correction/upstream cell determination in the flux term.

I consider this to be just a kludge, so if anyone has a better idea of
what the composition for the non-existing gas and oil phases is,
please tell me. (note that generic compositional models do not exhibit
this issue because the composition of all fluids is always fully
defined because each component is assumed to dissolve in every phase.)
---
 opm/core/simulator/initStateEquil.hpp      |  4 +--
 opm/core/simulator/initStateEquil_impl.hpp | 32 +++++++++++++++++++++-
 2 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
index cbda500ee..861c86a34 100644
--- a/opm/core/simulator/initStateEquil.hpp
+++ b/opm/core/simulator/initStateEquil.hpp
@@ -25,7 +25,7 @@
 #include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/simulator/EquilibrationHelpers.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/utility/RegionMapping.hpp>
 #include <opm/core/utility/Units.hpp>
@@ -163,7 +163,7 @@ namespace Opm
         phaseSaturations(const Grid&             grid,
                          const Region&           reg,
                          const CellRange&        cells,
-                         BlackoilPropertiesInterface& props,
+                         BlackoilPropertiesFromDeck& props,
                          const std::vector<double> swat_init,
                          std::vector< std::vector<double> >& phase_pressures);
 
diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp
index 20f56c5a4..1c0286dff 100644
--- a/opm/core/simulator/initStateEquil_impl.hpp
+++ b/opm/core/simulator/initStateEquil_impl.hpp
@@ -885,7 +885,7 @@ namespace Opm
      */
     template<class Grid>
     void initStateEquil(const Grid& grid,
-                        BlackoilPropertiesInterface& props,
+                        BlackoilPropertiesFromDeck& props,
                         const Opm::DeckConstPtr deck,
                         const Opm::EclipseStateConstPtr eclipseState,
                         const double gravity,
@@ -901,6 +901,36 @@ namespace Opm
         state.saturation() = Details::convertSats(isc.saturation());
         state.gasoilratio() = isc.rs();
         state.rv() = isc.rv();
+
+        // set the Rs and Rv factors to the saturated ones if both the gas and the oil
+        // saturations are zero. in most cases this does not matter, but when doing the
+        // gravity correction (in particular, when determining the upstream direction) it
+        // does.
+        if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
+            const int posOil = pu.phase_pos[BlackoilPhases::Liquid];
+            const int posGas = pu.phase_pos[BlackoilPhases::Vapour];
+            for (unsigned cellIdx = 0; cellIdx < state.pressure().size(); ++ cellIdx) {
+                const double so = isc.saturation()[posOil][cellIdx];
+                const double sg = isc.saturation()[posGas][cellIdx];
+
+                if (so <= 0.0 && sg <= 0.0) {
+                    int pvtRegionIdx = props.cellPvtRegionIndex()[cellIdx];
+                    const double po = state.pressure()[cellIdx];
+                    const double pg = po; // hack: include capillary pressure!
+                    const double T = 273.15 + 25.0; // does not matter for isothermal simulations!
+                    const double rsSat = props.oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx,
+                                                                                      T,
+                                                                                      pg);
+                    const double rvSat = props.gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx,
+                                                                                       T,
+                                                                                       po);
+
+                    state.gasoilratio()[cellIdx] = rsSat;
+                    state.rv()[cellIdx] = rvSat;
+                }
+            }
+        }
+
         initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
     }