Fix Family I SOGCR Endpoint Hack

The maximum attainable oil saturation in the G/O system is 1-SWL.
Update call sites accordingly.

opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp

@@ -186,8 +186,12 @@ public:
                                        Scalar& krnSwMdc,
                                        const Params& params)
     {
-        pcSwMdc = params.gasOilParams().pcSwMdc();
+        const auto Swco = params.Swl();
-        krnSwMdc = params.gasOilParams().krnSwMdc();
+
+        // Pretend oil saturation is 'Swco' larger than it really is in
+        // order to infer correct maximum Sg values in output layer.
+        pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
+        krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
 
         Valgrind::CheckDefined(pcSwMdc);
         Valgrind::CheckDefined(krnSwMdc);
@@ -203,8 +207,10 @@ public:
                                           const Scalar& krnSwMdc,
                                           Params& params)
     {
+        // Maximum attainable oil saturation is 1-SWL
+        const auto Swco = params.Swl();
         const double krwSw = 2.0; //Should not be used
-        params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc);
+        params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
     }
 
     /*!
@@ -220,7 +226,8 @@ public:
     static Evaluation pcgn(const Params& params,
                            const FluidState& fs)
     {
-        const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
+        // Maximum attainable oil saturation is 1-SWL.
+        const auto Sw = 1.0 - params.Swl() - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
     }
 
@@ -316,7 +323,8 @@ public:
     static Evaluation krg(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation Sw = 1.0 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        // Maximum attainable oil saturation is 1-SWL.
+        const Evaluation Sw = 1.0 - params.Swl() - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
     }
 
@@ -379,12 +387,13 @@ public:
     template <class FluidState>
     static void updateHysteresis(Params& params, const FluidState& fluidState)
     {
-        Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
+        const Scalar Swco = params.Swl();
-        Scalar So = Opm::scalarValue(fluidState.saturation(oilPhaseIdx));
+
-        Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
+        const Scalar Sw = clampSaturation(fluidState, waterPhaseIdx);
+        const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
+        const Scalar Sg = clampSaturation(fluidState, gasPhaseIdx);
 
         if (params.inconsistentHysteresisUpdate()) {
-            Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
             // NOTE: the saturations which are passed to update the hysteresis curves are
             // inconsistent with the ones used to calculate the relative permabilities. We do
             // it like this anyway because (a) the saturation functions of opm-core do it
@@ -394,21 +403,34 @@ public:
             //
             // Though be aware that from a physical perspective this is definitively
             // incorrect!
-            params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So);
+            params.oilWaterParams().update(/*pcSw=*/  1.0 - So,
-            params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
+                                           /*krwSw=*/ 1.0 - So,
+                                           /*krnSw=*/ 1.0 - So);
+
+            params.gasOilParams().update(/*pcSw=*/  1.0 - Swco - Sg,
+                                         /*krwSw=*/ 1.0 - Swco - Sg,
+                                         /*krnSw=*/ 1.0 - Swco - Sg);
         }
         else {
-            const Scalar Swco = params.Swl();
-            Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw));
-            Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
-
             const Scalar Sw_ow = Sg + std::max(Swco, Sw);
-            const Scalar So_go = 1 + Sw_ow;
+            const Scalar So_go = 1.0 - Sw_ow;
 
-            params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
+            params.oilWaterParams().update(/*pcSw=*/  Sw,
-            params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
+                                           /*krwSw=*/ 1 - Sg,
+                                           /*krnSw=*/ Sw_ow);
+
+            params.gasOilParams().update(/*pcSw=*/  1.0 - Swco - Sg,
+                                         /*krwSw=*/ So_go,
+                                         /*krnSw=*/ 1.0 - Swco - Sg);
         }
     }
+
+    template <class FluidState>
+    static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
+    {
+        const auto sat = Opm::scalarValue(fluidState.saturation(phaseIndex));
+        return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
+    }
 };
 } // namespace Opm
 

opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp

@@ -150,10 +150,9 @@ struct EclEpsScalingPointsInfo
      * I.e., the values which are used for the nested Fluid-Matrix interactions and which
      * are produced by them.
      */
-    void extractUnscaled(const Opm::satfunc::RawTableEndPoints&    rtep,
+    void extractUnscaled(const Opm::satfunc::RawTableEndPoints& rtep,
-                         const Opm::satfunc::RawFunctionValues&    rfunc,
+                         const Opm::satfunc::RawFunctionValues& rfunc,
-                         const Opm::SatFuncControls::KeywordFamily family,
+                         const std::vector<double>::size_type   satRegionIdx)
-                         const std::vector<double>::size_type      satRegionIdx)
     {
         this->Swl = rtep.connate.water[satRegionIdx];
         this->Sgl = rtep.connate.gas  [satRegionIdx];
@@ -162,10 +161,6 @@ struct EclEpsScalingPointsInfo
         this->Sgcr  = rtep.critical.gas         [satRegionIdx];
         this->Sowcr = rtep.critical.oil_in_water[satRegionIdx];
         this->Sogcr = rtep.critical.oil_in_gas  [satRegionIdx];
-        if (family == SatFuncControls::KeywordFamily::Family_I) {
-            // Hack.  Papers over unknown problems elsewhere.
-            this->Sogcr += this->Swl;
-        }
 
         this->Swu = rtep.maximum.water[satRegionIdx];
         this->Sgu = rtep.maximum.gas  [satRegionIdx];
@@ -352,21 +347,23 @@ public:
             assert(epsSystemType == EclGasOilSystem);
 
             // saturation scaling for capillary pressure
-            saturationPcPoints_[0] = 1.0 - epsInfo.Sgu;
+            saturationPcPoints_[0] = 1.0 - epsInfo.Swl - epsInfo.Sgu;
-            saturationPcPoints_[2] = saturationPcPoints_[1] = 1.0 - epsInfo.Sgl;
+            saturationPcPoints_[2] = saturationPcPoints_[1] = 1.0 - epsInfo.Swl - epsInfo.Sgl;
 
             // krw saturation scaling endpoints
             saturationKrwPoints_[0] = epsInfo.Sogcr;
-            saturationKrwPoints_[1] = 1 - epsInfo.Sgcr - epsInfo.Swl;
+            saturationKrwPoints_[1] = 1.0 - epsInfo.Sgcr - epsInfo.Swl;
-            saturationKrwPoints_[2] = 1 - epsInfo.Swl - epsInfo.Sgl;
+            saturationKrwPoints_[2] = 1.0 - epsInfo.Swl - epsInfo.Sgl;
 
             // krn saturation scaling endpoints (with the non-wetting phase being gas).
-            // because opm-material specifies non-wetting phase relperms in terms of the
+            //
-            // wetting phase saturations, the code here uses 1 minus the values specified
+            // As opm-material specifies non-wetting phase relative
-            // by the Eclipse TD and the order of the scaling points is reversed
+            // permeabilities in terms of the wetting phase saturations, the
-            saturationKrnPoints_[2] = 1.0 - epsInfo.Sgcr;
+            // code here uses (1-SWL) minus the values specified by the
-            saturationKrnPoints_[1] = epsInfo.Sogcr + epsInfo.Swl;
+            // ECLIPSE TD and the order of the scaling points is reversed.
-            saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
+            saturationKrnPoints_[2] = 1.0 - epsInfo.Swl - epsInfo.Sgcr;
+            saturationKrnPoints_[1] = epsInfo.Sogcr;
+            saturationKrnPoints_[0] = 1.0 - epsInfo.Swl - epsInfo.Sgu;
 
             if (config.enableLeverettScaling())
                 maxPcnwOrLeverettFactor_ = epsInfo.pcgoLeverettFactor;

opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp

@@ -163,13 +163,12 @@ public:
         const auto tolcrit = runspec.saturationFunctionControls()
             .minimumRelpermMobilityThreshold();
 
-        const auto family   = runspec.saturationFunctionControls().family();
         const auto rtepPtr  = satfunc::getRawTableEndpoints(tables, ph, tolcrit);
         const auto rfuncPtr = satfunc::getRawFunctionValues(tables, ph, *rtepPtr);
 
         for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
             this->unscaledEpsInfo_[satRegionIdx]
-                .extractUnscaled(*rtepPtr, *rfuncPtr, family, satRegionIdx);
+                .extractUnscaled(*rtepPtr, *rfuncPtr, satRegionIdx);
         }
     }
 
@@ -753,13 +752,11 @@ private:
     {
         // convert the saturations of the SGOF keyword from gas to oil saturations
         std::vector<double> SoSamples(sgofTable.numRows());
-        std::vector<double> SoKroSamples(sgofTable.numRows());
         for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) {
-            SoSamples[sampleIdx] = 1 - sgofTable.get("SG", sampleIdx);
+            SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx);
-            SoKroSamples[sampleIdx] = SoSamples[sampleIdx] - Swco;
         }
 
-        effParams.setKrwSamples(SoKroSamples, sgofTable.getColumn("KROG").vectorCopy());
+        effParams.setKrwSamples(SoSamples, sgofTable.getColumn("KROG").vectorCopy());
         effParams.setKrnSamples(SoSamples, sgofTable.getColumn("KRG").vectorCopy());
         effParams.setPcnwSamples(SoSamples, sgofTable.getColumn("PCOG").vectorCopy());
         effParams.finalize();
@@ -771,20 +768,18 @@ private:
     {
         // convert the saturations of the SLGOF keyword from "liquid" to oil saturations
         std::vector<double> SoSamples(slgofTable.numRows());
-        std::vector<double> SoKroSamples(slgofTable.numRows());
         for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) {
-            SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx);
+            SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco;
-            SoKroSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco;
         }
 
-        effParams.setKrwSamples(SoKroSamples, slgofTable.getColumn("KROG").vectorCopy());
+        effParams.setKrwSamples(SoSamples, slgofTable.getColumn("KROG").vectorCopy());
         effParams.setKrnSamples(SoSamples, slgofTable.getColumn("KRG").vectorCopy());
         effParams.setPcnwSamples(SoSamples, slgofTable.getColumn("PCOG").vectorCopy());
         effParams.finalize();
     }
 
     void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
-                                               Scalar /* Swco */,
+                                               Scalar Swco,
                                                const Opm::Sof3Table& sof3Table,
                                                const Opm::SgfnTable& sgfnTable)
     {
@@ -792,7 +787,7 @@ private:
         std::vector<double> SoSamples(sgfnTable.numRows());
         std::vector<double> SoColumn = sof3Table.getColumn("SO").vectorCopy();
         for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
-            SoSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx);
+            SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
         }
 
         effParams.setKrwSamples(SoColumn, sof3Table.getColumn("KROG").vectorCopy());
@@ -802,7 +797,7 @@ private:
     }
 
     void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
-                                               Scalar /* Swco */,
+                                               Scalar Swco,
                                                const Opm::Sof2Table& sof2Table,
                                                const Opm::SgfnTable& sgfnTable)
     {
@@ -810,7 +805,7 @@ private:
         std::vector<double> SoSamples(sgfnTable.numRows());
         std::vector<double> SoColumn = sof2Table.getColumn("SO").vectorCopy();
         for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
-            SoSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx);
+            SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
         }
 
         effParams.setKrwSamples(SoColumn, sof2Table.getColumn("KRO").vectorCopy());

opm/material/fluidmatrixinteractions/EclStone1Material.hpp

@@ -186,8 +186,12 @@ public:
                                        Scalar& krnSwMdc,
                                        const Params& params)
     {
-        pcSwMdc = params.gasOilParams().pcSwMdc();
+        const auto Swco = params.Swl();
-        krnSwMdc = params.gasOilParams().krnSwMdc();
+
+        // Pretend oil saturation is 'Swco' larger than it really is in
+        // order to infer correct maximum Sg values in output layer.
+        pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
+        krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
 
         Valgrind::CheckDefined(pcSwMdc);
         Valgrind::CheckDefined(krnSwMdc);
@@ -203,8 +207,10 @@ public:
                                           const Scalar& krnSwMdc,
                                           Params& params)
     {
+        // Maximum attainable oil saturation is 1-SWL
+        const auto Swco = params.Swl();
         const double krwSw = 2.0; //Should not be used
-        params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc);
+        params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
     }
 
     /*!
@@ -220,7 +226,8 @@ public:
     static Evaluation pcgn(const Params& params,
                            const FluidState& fs)
     {
-        const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
+        // Maximum attainable oil saturation is 1-SWL
+        const auto Sw = 1.0 - params.Swl() - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
     }
 
@@ -321,7 +328,8 @@ public:
     static Evaluation krg(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        // Maximum attainable oil saturation is 1-SWL,
+        const Evaluation Sw = 1 - params.Swl() - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
     }
 
@@ -387,11 +395,14 @@ public:
     template <class FluidState>
     static void updateHysteresis(Params& params, const FluidState& fluidState)
     {
+        const Scalar Swco = params.Swl();
         const Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
         const Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
 
         params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
-        params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg);
+        params.gasOilParams().update(/*pcSw=*/  1.0 - Swco - Sg,
+                                     /*krwSw=*/ 1.0 - Swco - Sg,
+                                     /*krnSw=*/ 1.0 - Swco - Sg);
     }
 };
 } // namespace Opm

opm/material/fluidmatrixinteractions/EclStone2Material.hpp

@@ -187,8 +187,12 @@ public:
                                        Scalar& krnSwMdc,
                                        const Params& params)
     {
-        pcSwMdc = params.gasOilParams().pcSwMdc();
+        const auto Swco = params.Swl();
-        krnSwMdc = params.gasOilParams().krnSwMdc();
+
+        // Pretend oil saturation is 'Swco' larger than it really is in
+        // order to infer correct maximum Sg values in output layer.
+        pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
+        krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
 
         Valgrind::CheckDefined(pcSwMdc);
         Valgrind::CheckDefined(krnSwMdc);
@@ -204,8 +208,10 @@ public:
                                           const Scalar& krnSwMdc,
                                           Params& params)
     {
+        // Maximum attainable oil saturation is 1-SWL
+        const auto Swco = params.Swl();
         const double krwSw = 2.0; //Should not be used
-        params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc);
+        params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
     }
 
     /*!
@@ -221,7 +227,8 @@ public:
     static Evaluation pcgn(const Params& params,
                            const FluidState& fs)
     {
-        const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
+        // Maximum attainable oil saturation is 1-SWL.
+        const auto Sw = 1.0 - params.Swl() - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
     }
 
@@ -322,7 +329,8 @@ public:
     static Evaluation krg(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        // Maximum attainable oil saturation is 1-SWL.
+        const Evaluation Sw = 1 - params.Swl() - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
     }
 
@@ -351,8 +359,8 @@ public:
         const Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco);
         const Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
         const Evaluation krw = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
-        const Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Sg);
+        const Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Swco - Sg);
-        const Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg);
+        const Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Swco - Sg);
 
         return krocw * ((krow/krocw + krw) * (krog/krocw + krg) - krw - krg);
     }
@@ -367,11 +375,14 @@ public:
     template <class FluidState>
     static void updateHysteresis(Params& params, const FluidState& fluidState)
     {
-        Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
+        const Scalar Swco = params.Swl();
-        Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
+        const Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
+        const Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
 
         params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
-        params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg);
+        params.gasOilParams().update(/*pcSw=*/  1.0 - Swco - Sg,
+                                     /*krwSw=*/ 1.0 - Swco - Sg,
+                                     /*krnSw=*/ 1.0 - Swco - Sg);
     }
 };
 } // namespace Opm