From 56e829329f3d873830a56b6a79b766065bb5901d Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Tue, 22 Jan 2019 12:50:47 +0100
Subject: [PATCH] each phase needs to be above certain value to be treated to
 be present

it helps to recover some RESTART running based on single precision
format.
---
 ebos/eclproblem.hh | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index 7f708b4b5..984ec59be 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -1868,6 +1868,9 @@ private:
             elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
             eclWriter_->eclOutputModule().initHysteresisParams(this->simulator(), elemIdx);
             eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
+
+            processRestartSaturations_(elemFluidState);
+
             lastRs_[elemIdx] = elemFluidState.Rs();
             lastRv_[elemIdx] = elemFluidState.Rv();
             if (enableSolvent)
@@ -1894,6 +1897,30 @@ private:
 
     }
 
+    void processRestartSaturations_(InitialFluidState& elemFluidState) {
+        // 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
+        const Scalar smallSaturationTolerance = 1.e-6;
+        Scalar sumSaturation = 0.;
+        for (size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (FluidSystem::phaseIsActive(phaseIdx)) {
+                if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
+                    elemFluidState.setSaturation(phaseIdx, 0.);
+
+                sumSaturation += elemFluidState.saturation(phaseIdx);
+            }
+        }
+
+        assert(sumSaturation > 0.);
+
+        for (size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (FluidSystem::phaseIsActive(phaseIdx)) {
+                const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
+                elemFluidState.setSaturation(phaseIdx, saturation);
+            }
+        }
+    }
+
     void readExplicitInitialCondition_()
     {
         const auto& vanguard = this->simulator().vanguard();