From da06e423d506630d8aa7823b188177993b219bd1 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Wed, 16 Aug 2023 16:19:41 +0200
Subject: [PATCH] add parameter to set threshold for only water cells. Default
 is 1.0

---
 opm/models/blackoil/blackoilnewtonmethod.hh     | 16 ++++++++++++++--
 opm/models/blackoil/blackoilprimaryvariables.hh | 11 +++++------
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/opm/models/blackoil/blackoilnewtonmethod.hh b/opm/models/blackoil/blackoilnewtonmethod.hh
index 17cc57d3a..a849de11a 100644
--- a/opm/models/blackoil/blackoilnewtonmethod.hh
+++ b/opm/models/blackoil/blackoilnewtonmethod.hh
@@ -57,6 +57,8 @@ template<class TypeTag, class MyTypeTag>
 struct TemperatureMin { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct MaximumWaterSaturation { using type = UndefinedProperty; };
+template<class TypeTag, class MyTypeTag>
+struct WaterOnlyThreshold { using type = UndefinedProperty; };
 template<class TypeTag>
 struct DpMaxRel<TypeTag, TTag::NewtonMethod>
 {
@@ -101,6 +103,12 @@ struct MaximumWaterSaturation<TypeTag, TTag::NewtonMethod>
     using type = GetPropType<TypeTag, Scalar>;
     static constexpr type value = 1.0;
 };
+template<class TypeTag>
+struct WaterOnlyThreshold<TypeTag, TTag::NewtonMethod>
+{
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1.0;
+};
 } // namespace Opm::Properties
 
 namespace Opm {
@@ -139,6 +147,7 @@ public:
         tempMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, TemperatureMax);
         tempMin_ = EWOMS_GET_PARAM(TypeTag, Scalar, TemperatureMin);
         waterSaturationMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaximumWaterSaturation);
+        waterOnlyThreshold_ = EWOMS_GET_PARAM(TypeTag, Scalar, WaterOnlyThreshold);
     }
 
     /*!
@@ -168,6 +177,7 @@ public:
         EWOMS_REGISTER_PARAM(TypeTag, Scalar, TemperatureMax, "Maximum absolute temperature");
         EWOMS_REGISTER_PARAM(TypeTag, Scalar, TemperatureMin, "Minimum absolute temperature");
         EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaximumWaterSaturation, "Maximum water saturation");
+        EWOMS_REGISTER_PARAM(TypeTag, Scalar, WaterOnlyThreshold, "Cells with water saturation above or equal is considered one-phase water only");
     }
 
     /*!
@@ -458,9 +468,9 @@ protected:
         // use a threshold value after a switch to make it harder to switch back
         // immediately.
         if (wasSwitched_[globalDofIdx])
-            wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, waterSaturationMax_, priVarOscilationThreshold_);
+            wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, waterSaturationMax_, waterOnlyThreshold_, priVarOscilationThreshold_);
         else
-            wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, waterSaturationMax_);
+            wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, waterSaturationMax_, waterOnlyThreshold_);
 
         if (wasSwitched_[globalDofIdx])
             ++ numPriVarsSwitched_;
@@ -476,6 +486,8 @@ private:
 
     Scalar priVarOscilationThreshold_;
     Scalar waterSaturationMax_;
+    Scalar waterOnlyThreshold_;
+
     Scalar dpMaxRel_;
     Scalar dsMax_;
     bool projectSaturations_;
diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh
index f0750276d..b7925e2fb 100644
--- a/opm/models/blackoil/blackoilprimaryvariables.hh
+++ b/opm/models/blackoil/blackoilprimaryvariables.hh
@@ -493,10 +493,8 @@ public:
      *
      * \return true Iff the interpretation of one of the switching variables was changed
      */
-    bool adaptPrimaryVariables(const Problem& problem, unsigned globalDofIdx, Scalar swMaximum, Scalar eps = 0.0)
+    bool adaptPrimaryVariables(const Problem& problem, unsigned globalDofIdx, Scalar swMaximum, Scalar thresholdWaterFilledCell, Scalar eps = 0.0)
     {
-        static const Scalar thresholdWaterFilledCell = 1.0 - eps;
-
         // this function accesses quite a few black-oil specific low-level functions
         // directly for better performance (instead of going the canonical way through
         // the IntensiveQuantities). The reason is that most intensive quantities are not
@@ -546,9 +544,10 @@ public:
         // keep track if any meaning has changed
         bool changed = false;
 
-        // special case cells with almost only water
-        // use both saturations (if the phase is enabled)
-        // if dissolved gas in water is enabled we shouldn't enter
+        // Special case for cells with almost only water
+        // for these cells both saturations (if the phase is enabled) is used
+        // to avoid singular systems.
+        // If dissolved gas in water is enabled we shouldn't enter
         // here but instead switch to Rsw as primary variable
         // as sw >= 1.0 -> gas <= 0 (i.e. gas phase disappears)
         if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {