From bbf6d56000e3103ebfe73331182a343cd90f3627 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <thsa@iris.no>
Date: Mon, 8 Sep 2014 11:14:44 +0200
Subject: [PATCH] Calculate oil saturation from changes in water and gas
 saturation First the change in oil saturation is calculated from changes in
 water and oil saturation. Then oil saturation is updated based on this change
 instead of just fixed to 1-sw-sg. With this change the oil saturation is less
 sensitive towards numerical errors that may cause very small oil saturations.
 Witch again may cause the simulator to think that the gas phase is saturation
 with vaporized oil when it is not.

---
 .../FullyImplicitBlackoilSolver_impl.hpp      | 74 ++++++++++---------
 1 file changed, 40 insertions(+), 34 deletions(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 7fb74d86d..256d9a12b 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -1278,7 +1278,7 @@ namespace {
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
         const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
         const double dsmax = dsMax();
-        V so = one;
+        V so;
         V sw;
         V sg;
 
@@ -1306,17 +1306,52 @@ namespace {
                 const int pos = pu.phase_pos[ Water ];
                 const V sw_old = s_old.col(pos);
                 sw = sw_old - step * dsw;
-                so -= sw;
             }
 
             if (active_[Gas]) {
                 const int pos = pu.phase_pos[ Gas ];
                 const V sg_old = s_old.col(pos);
                 sg = sg_old - step * dsg;
-                so -= sg;
+            }
+
+            const int pos = pu.phase_pos[ Oil ];
+            const V so_old = s_old.col(pos);
+            so = so_old - step * dso;
+        }
+
+        auto ixg = sg < 0;
+        for (int c = 0; c < nc; ++c) {
+            if (ixg[c]) {
+                sw[c] = sw[c] / (1-sg[c]);
+                so[c] = so[c] / (1-sg[c]);
+                sg[c] = 0;
             }
         }
 
+
+        auto ixo = so < 0;
+        for (int c = 0; c < nc; ++c) {
+            if (ixo[c]) {
+                sw[c] = sw[c] / (1-so[c]);
+                sg[c] = sg[c] / (1-so[c]);
+                so[c] = 0;
+            }
+        }
+
+        auto ixw = sw < 0;
+        for (int c = 0; c < nc; ++c) {
+            if (ixw[c]) {
+                so[c] = so[c] / (1-sw[c]);
+                sg[c] = sg[c] / (1-sw[c]);
+                sw[c] = 0;
+            }
+        }
+
+        const V sumSat = sw + so + sg;
+        sw = sw / sumSat;
+        so = so / sumSat;
+        sg = sg / sumSat;
+
         const double drsmaxrel = drsMaxRel();
         const double drvmax = 1e9;//% same as in Mrst
         V rs;
@@ -1351,7 +1386,7 @@ namespace {
 
             // keep oil saturated if previous sg is sufficient large:
             const int pos = pu.phase_pos[ Gas ];
-            auto hadGas = (sg < 0 && s_old.col(pos) > epsilon);
+            auto hadGas = (sg <= 0 && s_old.col(pos) > epsilon);
             // Set oil saturated if previous rs is sufficiently large
             const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
             auto gasVaporized =  ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
@@ -1374,7 +1409,7 @@ namespace {
 
             // keep oil saturated if previous so is sufficient large:
             const int pos = pu.phase_pos[ Oil ];
-            auto hadOil = (so < 0 && s_old.col(pos) > epsilon );
+            auto hadOil = (so <= 0 && s_old.col(pos) > epsilon );
             // Set oil saturated if previous rv is sufficiently large
             const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc);
             auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) );
@@ -1387,35 +1422,6 @@ namespace {
 
         }
 
-        auto ixg = sg < 0;
-        for (int c = 0; c < nc; ++c) {
-            if (ixg[c]) {
-                sw[c] = sw[c] / (1-sg[c]);
-                so[c] = so[c] / (1-sg[c]);
-                sg[c] = 0;
-            }
-        }
-
-
-        auto ixo = so < 0;
-        for (int c = 0; c < nc; ++c) {
-            if (ixo[c]) {
-                sw[c] = sw[c] / (1-so[c]);
-                sg[c] = sg[c] / (1-so[c]);
-                so[c] = 0;
-            }
-        }
-
-        auto ixw = sw < 0;
-        for (int c = 0; c < nc; ++c) {
-            if (ixw[c]) {
-                so[c] = so[c] / (1-sw[c]);
-                sg[c] = sg[c] / (1-sw[c]);
-                sw[c] = 0;
-            }
-        }
-
-
         // Update saturations
 
         for (int c = 0; c < nc; ++c) {