From afba6c8e8eb917eceeec0d742df2393ae90b5fb3 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Wed, 25 Sep 2019 12:28:58 +0200
Subject: [PATCH] Fix solvent restart

---
 ebos/ecloutputblackoilmodule.hh | 19 +++++++++++--------
 ebos/eclproblem.hh              | 20 ++++++++++++++++----
 ebos/eclwriter.hh               |  3 +++
 3 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh
index 98c6e6758..0fe5bf301 100644
--- a/ebos/ecloutputblackoilmodule.hh
+++ b/ebos/ecloutputblackoilmodule.hh
@@ -1041,6 +1041,17 @@ public:
             so -= sol.data("SGAS")[globalDofIndex];
         }
 
+        if (sSol_.size() > 0) {
+            // keep the SSOL option for backward compatibility
+            // should be removed after 10.2018 release
+            if (sol.has("SSOL"))
+                sSol_[elemIdx] = sol.data("SSOL")[globalDofIndex];
+            else if (sol.has("SSOLVENT"))
+                sSol_[elemIdx] = sol.data("SSOLVENT")[globalDofIndex];
+
+            so -= sSol_[elemIdx];
+        }
+
         assert(saturation_[oilPhaseIdx].size() > 0);
         saturation_[oilPhaseIdx][elemIdx] = so;
 
@@ -1052,14 +1063,6 @@ public:
             rs_[elemIdx] = sol.data("RS")[globalDofIndex];
         if (rv_.size() > 0 && sol.has("RV"))
             rv_[elemIdx] = sol.data("RV")[globalDofIndex];
-        if (sSol_.size() > 0) {
-            // keep the SSOL option for backward compatibility
-            // should be removed after 10.2018 release
-            if (sol.has("SSOL"))
-                sSol_[elemIdx] = sol.data("SSOL")[globalDofIndex];
-            else if (sol.has("SSOLVENT"))
-                sSol_[elemIdx] = sol.data("SSOLVENT")[globalDofIndex];
-        }
         if (cPolymer_.size() > 0 && sol.has("POLYMER"))
             cPolymer_[elemIdx] = sol.data("POLYMER")[globalDofIndex];
         if (cFoam_.size() > 0 && sol.has("FOAM"))
diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index d1ee4e5b6..154f7caf7 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -2492,12 +2492,14 @@ private:
             eclWriter_->eclOutputModule().initHysteresisParams(simulator, elemIdx);
             eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
 
-            processRestartSaturations_(elemFluidState);
+            if (enableSolvent)
+                 solventSaturation_[elemIdx] = eclWriter_->eclOutputModule().getSolventSaturation(elemIdx);
+
+            processRestartSaturations_(elemFluidState, solventSaturation_[elemIdx]);
 
             lastRs_[elemIdx] = elemFluidState.Rs();
             lastRv_[elemIdx] = elemFluidState.Rv();
-            if (enableSolvent)
-                 solventSaturation_[elemIdx] = eclWriter_->eclOutputModule().getSolventSaturation(elemIdx);
+
             if (enablePolymer)
                  polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
             // if we need to restart for polymer molecular weight simulation, we need to add related here
@@ -2548,7 +2550,7 @@ private:
         eclWriter_->endRestart();
     }
 
-    void processRestartSaturations_(InitialFluidState& elemFluidState)
+    void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
     {
         // each phase needs to be above certain value to be claimed to be existing
         // this is used to recover some RESTART running with the defaulted single-precision format
@@ -2561,6 +2563,13 @@ private:
 
                 sumSaturation += elemFluidState.saturation(phaseIdx);
             }
+
+        }
+        if (enableSolvent) {
+            if (solventSaturation < smallSaturationTolerance)
+                solventSaturation = 0.0;
+
+           sumSaturation += solventSaturation;
         }
 
         assert(sumSaturation > 0.0);
@@ -2571,6 +2580,9 @@ private:
                 elemFluidState.setSaturation(phaseIdx, saturation);
             }
         }
+        if (enableSolvent) {
+            solventSaturation = solventSaturation / sumSaturation;
+        }
     }
 
     void readExplicitInitialCondition_()
diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh
index 0b173f1dd..c9f737dcb 100644
--- a/ebos/eclwriter.hh
+++ b/ebos/eclwriter.hh
@@ -157,6 +157,8 @@ class EclWriter
     typedef std::vector<Scalar> ScalarBuffer;
 
     enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+
 
 public:
     static void registerParameters()
@@ -392,6 +394,7 @@ public:
             {"SWAT", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))},
             {"SGAS", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))},
             {"TEMP" , Opm::UnitSystem::measure::temperature, enableEnergy},
+            {"SSOLVENT" , Opm::UnitSystem::measure::identity, enableSolvent},
             {"RS", Opm::UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
             {"RV", Opm::UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
             {"SOMAX", Opm::UnitSystem::measure::identity, simulator_.problem().vapparsActive()},