Fixes to the special case with almost only water

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -307,7 +307,11 @@ protected:
             else if (pvIdx == Indices::waterSwitchIdx)
                 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
                     delta *= satAlpha;
+
                 else {
+                    // The damping from the saturation change also applies to the rsw/rvw factors
+                    // for consistency
+                    delta *= satAlpha;
                     //Ensure Rvw and Rsw factor does not become negative
                     if (delta > currentValue[ Indices::waterSwitchIdx]) 
                         delta = currentValue[ Indices::waterSwitchIdx];
@@ -321,6 +325,9 @@ protected:
                 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
                     delta *= satAlpha;
                 else {
+                    // The damping from the saturation change also applies to the rs/rv factors
+                    // for consistency
+                    delta *= satAlpha;
                     //Ensure Rv and Rs factor does not become negative
                     if (delta > currentValue[Indices::compositionSwitchIdx])
                         delta = currentValue[Indices::compositionSwitchIdx];
@@ -330,6 +337,7 @@ protected:
                 // solvent saturation updates are also subject to the Appleyard chop
                 delta *= satAlpha;
             else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
+                delta *= satAlpha;
                 // z fraction updates are also subject to the Appleyard chop
                 if (delta > currentValue[Indices::zFractionIdx])
                         delta = currentValue[Indices::zFractionIdx];

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -495,8 +495,6 @@ public:
      */
     bool adaptPrimaryVariables(const Problem& problem, unsigned globalDofIdx, 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
@@ -548,25 +546,30 @@ public:
 
         // special case cells with almost only water
         // use both saturations (if the phase is enabled)
+        // to avoid a singular matrix.
+        // almost is defined as sw >= 1.0-1e-6.
+        static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
+
         // 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)
+        // as sw >= 1.0 -> sg <= 0 (i.e. gas phase disappears)
         if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
-
-            // make sure water saturations does not exceed 1.0
+            // make sure water saturations does not exceed 1.01
+            // we dont stricly force it to 1.0 to avoid convergence issue.
             if constexpr (waterEnabled) {
-                (*this)[Indices::waterSwitchIdx] = 1.0;
+                (*this)[Indices::waterSwitchIdx] = std::min(1.01, sw);
                 assert(primaryVarsMeaningWater() == WaterMeaning::Sw);
             }
+
             // the hydrocarbon gas saturation is set to 0.0
             if constexpr (compositionSwitchEnabled)
                 (*this)[Indices::compositionSwitchIdx] = 0.0;
 
             changed = primaryVarsMeaningGas() != GasMeaning::Sg;
             if(changed) {
-                if constexpr (compositionSwitchEnabled)
+                if constexpr (compositionSwitchEnabled) {
                     setPrimaryVarsMeaningGas(GasMeaning::Sg);
-
+                }
                 // use water pressure?
             }
             return changed;