Fixes to the special case with almost only water

For the almost only water case the primary variables are p, sg and sw.

1) Use sw>=1-1e-6 as threshold for defining almost only water
2) Add sg to sw if sg<0 to make the oil saturation concistent even though sg is set to zero
3) Chop the water saturation at 1.01 and not 1.0 in order to improve convergence
4) Restrict update of rs/rv and rsw/rvw by the saturation scaling factor from the Appleyard-chopping.

Author: Tor Harald Sandve

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

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,33 @@ 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)
-        if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
-
-            // make sure water saturations does not exceed 1.0
+        // as sw >= 1.0 -> sg <= 0 (i.e. gas phase disappears)
+        if ((sw - std::min(sg, 0.0)) >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
+            // make sure water saturations does not exceed 1.01
+            // we dont stricly force it to 1.0 to avoid convergence issue.
+            // note that for sg < 0 we add the "negative" saturation of the gas to the water
+            // phase. Since So = 1 - sg - sw. This will give the correct oil saturation.
+            // even though the gas saturation is set to zero below.
             if constexpr (waterEnabled) {
-                (*this)[Indices::waterSwitchIdx] = 1.0;
+                (*this)[Indices::waterSwitchIdx] = std::min(1.01, sw - std::min(sg, 0.0));
                 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;