From c390e356f3bc19a6e5e73956a3b20e90463a2809 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 11 Jan 2021 18:04:18 +0100 Subject: [PATCH] Fix Family I SOGCR Endpoint Hack The maximum attainable oil saturation in the G/O system is 1-SWL. Update call sites accordingly. --- .../EclDefaultMaterial.hpp | 58 +++++++++++++------ .../EclEpsScalingPoints.hpp | 33 +++++------ .../EclMaterialLawManager.hpp | 23 +++----- .../EclStone1Material.hpp | 23 ++++++-- .../EclStone2Material.hpp | 31 ++++++---- 5 files changed, 102 insertions(+), 66 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index 825978d13..7946eafd2 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -186,8 +186,12 @@ public: Scalar& krnSwMdc, const Params& params) { - pcSwMdc = params.gasOilParams().pcSwMdc(); - krnSwMdc = params.gasOilParams().krnSwMdc(); + const auto Swco = params.Swl(); + + // 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(fs.saturation(gasPhaseIdx)); + // Maximum attainable oil saturation is 1-SWL. + const auto Sw = 1.0 - params.Swl() - Opm::decay(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(fluidState.saturation(gasPhaseIdx)); + // Maximum attainable oil saturation is 1-SWL. + const Evaluation Sw = 1.0 - params.Swl() - Opm::decay(fluidState.saturation(gasPhaseIdx)); return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); } @@ -379,12 +387,13 @@ public: template static void updateHysteresis(Params& params, const FluidState& fluidState) { - Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx)); - Scalar So = Opm::scalarValue(fluidState.saturation(oilPhaseIdx)); - Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx)); + const Scalar Swco = params.Swl(); + + 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.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg); + params.oilWaterParams().update(/*pcSw=*/ 1.0 - So, + /*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.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg); + params.oilWaterParams().update(/*pcSw=*/ Sw, + /*krwSw=*/ 1 - Sg, + /*krnSw=*/ Sw_ow); + + params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg, + /*krwSw=*/ So_go, + /*krnSw=*/ 1.0 - Swco - Sg); } } + + template + 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 diff --git a/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp b/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp index 09a634f8a..3e3e08a24 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp +++ b/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, - const Opm::satfunc::RawFunctionValues& rfunc, - const Opm::SatFuncControls::KeywordFamily family, - const std::vector::size_type satRegionIdx) + void extractUnscaled(const Opm::satfunc::RawTableEndPoints& rtep, + const Opm::satfunc::RawFunctionValues& rfunc, + const std::vector::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_[2] = saturationPcPoints_[1] = 1.0 - epsInfo.Sgl; + saturationPcPoints_[0] = 1.0 - epsInfo.Swl - epsInfo.Sgu; + 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_[2] = 1 - epsInfo.Swl - epsInfo.Sgl; + saturationKrwPoints_[1] = 1.0 - epsInfo.Sgcr - epsInfo.Swl; + 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 - // by the Eclipse TD and the order of the scaling points is reversed - saturationKrnPoints_[2] = 1.0 - epsInfo.Sgcr; - saturationKrnPoints_[1] = epsInfo.Sogcr + epsInfo.Swl; - saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu; + // + // As opm-material specifies non-wetting phase relative + // permeabilities in terms of the wetting phase saturations, the + // code here uses (1-SWL) minus the values specified by the + // ECLIPSE TD and the order of the scaling points is reversed. + 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; diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index faa18d71c..504c44988 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/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 SoSamples(sgofTable.numRows()); - std::vector SoKroSamples(sgofTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) { - SoSamples[sampleIdx] = 1 - sgofTable.get("SG", sampleIdx); - SoKroSamples[sampleIdx] = SoSamples[sampleIdx] - Swco; + SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx); } - 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 SoSamples(slgofTable.numRows()); - std::vector SoKroSamples(slgofTable.numRows()); for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) { - SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx); - SoKroSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco; + SoSamples[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 SoSamples(sgfnTable.numRows()); std::vector 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 SoSamples(sgfnTable.numRows()); std::vector 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()); diff --git a/opm/material/fluidmatrixinteractions/EclStone1Material.hpp b/opm/material/fluidmatrixinteractions/EclStone1Material.hpp index d893b16a3..24e6f4198 100644 --- a/opm/material/fluidmatrixinteractions/EclStone1Material.hpp +++ b/opm/material/fluidmatrixinteractions/EclStone1Material.hpp @@ -186,8 +186,12 @@ public: Scalar& krnSwMdc, const Params& params) { - pcSwMdc = params.gasOilParams().pcSwMdc(); - krnSwMdc = params.gasOilParams().krnSwMdc(); + const auto Swco = params.Swl(); + + // 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(fs.saturation(gasPhaseIdx)); + // Maximum attainable oil saturation is 1-SWL + const auto Sw = 1.0 - params.Swl() - Opm::decay(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(fluidState.saturation(gasPhaseIdx)); + // Maximum attainable oil saturation is 1-SWL, + const Evaluation Sw = 1 - params.Swl() - Opm::decay(fluidState.saturation(gasPhaseIdx)); return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); } @@ -387,11 +395,14 @@ public: template 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 diff --git a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp index 4d0e82742..1b94fedeb 100644 --- a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp +++ b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp @@ -187,8 +187,12 @@ public: Scalar& krnSwMdc, const Params& params) { - pcSwMdc = params.gasOilParams().pcSwMdc(); - krnSwMdc = params.gasOilParams().krnSwMdc(); + const auto Swco = params.Swl(); + + // 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(fs.saturation(gasPhaseIdx)); + // Maximum attainable oil saturation is 1-SWL. + const auto Sw = 1.0 - params.Swl() - Opm::decay(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(fluidState.saturation(gasPhaseIdx)); + // Maximum attainable oil saturation is 1-SWL. + const Evaluation Sw = 1 - params.Swl() - Opm::decay(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 krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg); + const Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Swco - 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 static void updateHysteresis(Params& params, const FluidState& fluidState) { - Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx)); - Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx)); + 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