add parameter to set threshold for only water cells. Default is 1.0

Author: Tor Harald Sandve

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_;

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,8 +544,11 @@ 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)
+        // special case for cells with almost only water
+        thresholdWaterFilledCell = std::min(thresholdWaterFilledCell, 1.0 - eps);
+        // 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)