From 64bdc84a9f339e76faed64ec6fee31a424006386 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Fri, 24 Feb 2017 17:09:40 +0100 Subject: [PATCH] Add option to store/restore hysteresis parameters --- .../EclDefaultMaterial.hpp | 62 ++++ .../EclHysteresisTwoPhaseLawParams.hpp | 21 +- .../EclMaterialLawManager.hpp | 44 +++ .../EclMultiplexerMaterial.hpp | 132 +++++++++ .../EclStone1Material.hpp | 62 ++++ .../EclStone2Material.hpp | 62 ++++ .../EclTwoPhaseMaterial.hpp | 62 ++++ tests/test_eclmateriallawmanager.cpp | 272 +++++++++++++----- 8 files changed, 636 insertions(+), 81 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index 4474d58d3..768ffe701 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -148,6 +148,68 @@ public: Valgrind::CheckDefined(values[waterPhaseIdx]); } + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.oilWaterParams().pcSwMdc(); + krnSwMdc = params.oilWaterParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const double krwSw = 2.0; //Should not be used + params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.gasOilParams().pcSwMdc(); + krnSwMdc = params.gasOilParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const double krwSw = 2.0; //Should not be used + params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc); + } + /*! * \brief Capillary pressure between the gas and the non-wetting * liquid (i.e., oil) phase. diff --git a/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLawParams.hpp b/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLawParams.hpp index b7f642f26..6bf9bc597 100644 --- a/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLawParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLawParams.hpp @@ -60,6 +60,9 @@ public: EclHysteresisTwoPhaseLawParams() { + // These are initialized to two (even though they represent saturations) + // to signify that the values are larger than physically possible, and force + // using the drainage curve before the first saturation update. pcSwMdc_ = 2.0; krnSwMdc_ = 2.0; // krwSwMdc_ = 2.0; @@ -152,7 +155,7 @@ public: { pcSwMdc_ = value; } /*! - * \brief Set the saturation of the wetting phase where the last switch from the main + * \brief Get the saturation of the wetting phase where the last switch from the main * drainage curve to imbibition happend on the capillary pressure curve. */ Scalar pcSwMdc() const @@ -168,7 +171,7 @@ public: // { krwSwMdc_ = value; }; /*! - * \brief Set the saturation of the wetting phase where the last switch from the main + * \brief Get the saturation of the wetting phase where the last switch from the main * drainage curve to imbibition happend on the relperm curve for the * wetting phase. */ @@ -185,7 +188,7 @@ public: { krnSwMdc_ = value; } /*! - * \brief Set the saturation of the wetting phase where the last switch from the main + * \brief Get the saturation of the wetting phase where the last switch from the main * drainage curve to imbibition happend on the relperm curve for the * non-wetting phase. */ @@ -286,9 +289,9 @@ private: Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage); deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_; - Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_); - Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage); - deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_; + // Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_); + // Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage); + // deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_; // assert(std::abs(EffLawT::twoPhaseSatPcnw(imbibitionParams(), pcSwMdc_ + deltaSwImbPc_) // - EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_)) < 1e-8); @@ -311,16 +314,16 @@ private: // largest wettinging phase saturation which is on the main-drainage curve. These are // three different values because the sourounding code can choose to use different // definitions for the saturations for different quantities -// Scalar krwSwMdc_; + //Scalar krwSwMdc_; Scalar krnSwMdc_; Scalar pcSwMdc_; // offsets added to wetting phase saturation uf using the imbibition curves need to // be used to calculate the wetting phase relperm, the non-wetting phase relperm and // the capillary pressure -// Scalar deltaSwImbKrw_; + //Scalar deltaSwImbKrw_; Scalar deltaSwImbKrn_; - Scalar deltaSwImbPc_; + //Scalar deltaSwImbPc_; // trapped non-wetting phase saturation //Scalar Sncrt_; diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index 2df1fab1e..1f40c9072 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp @@ -242,6 +242,50 @@ public: MaterialLaw::updateHysteresis(*threePhaseParams, fluidState); } + void oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + unsigned elemIdx) const + { + if (!enableHysteresis()) { + OPM_THROW(std::runtime_error, "Cannot get hysteresis parameters if hysteresis not enabled."); + } + const auto& params = materialLawParams(elemIdx); + MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params); + } + + void setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + unsigned elemIdx) + { + if (!enableHysteresis()) { + OPM_THROW(std::runtime_error, "Cannot set hysteresis parameters if hysteresis not enabled."); + } + auto& params = materialLawParams(elemIdx); + MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params); + } + + void gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + unsigned elemIdx) const + { + if (!enableHysteresis()) { + OPM_THROW(std::runtime_error, "Cannot get hysteresis parameters if hysteresis not enabled."); + } + const auto& params = materialLawParams(elemIdx); + MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params); + } + + void setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + unsigned elemIdx) + { + if (!enableHysteresis()) { + OPM_THROW(std::runtime_error, "Cannot set hysteresis parameters if hysteresis not enabled."); + } + auto& params = materialLawParams(elemIdx); + MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params); + } + EclEpsScalingPoints& oilWaterScaledEpsPointsDrainage(unsigned elemIdx) { auto& materialParams = *materialLawParams_[elemIdx]; diff --git a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp index ad7f60edf..212fab2af 100644 --- a/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp @@ -159,6 +159,138 @@ public: } } + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + switch (params.approach()) { + case EclStone1Approach: + Stone1Material::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclStone2Approach: + Stone2Material::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclDefaultApproach: + DefaultMaterial::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclTwoPhaseApproach: + TwoPhaseMaterial::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + } + } + + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + switch (params.approach()) { + case EclStone1Approach: + Stone1Material::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclStone2Approach: + Stone2Material::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclDefaultApproach: + DefaultMaterial::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclTwoPhaseApproach: + TwoPhaseMaterial::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + } + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + switch (params.approach()) { + case EclStone1Approach: + Stone1Material::gasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclStone2Approach: + Stone2Material::gasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclDefaultApproach: + DefaultMaterial::gasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclTwoPhaseApproach: + TwoPhaseMaterial::gasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + } + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + switch (params.approach()) { + case EclStone1Approach: + Stone1Material::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclStone2Approach: + Stone2Material::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclDefaultApproach: + DefaultMaterial::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + + case EclTwoPhaseApproach: + TwoPhaseMaterial::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, + params.template getRealParams()); + break; + } + } + /*! * \brief Capillary pressure between the gas and the non-wetting * liquid (i.e., oil) phase. diff --git a/opm/material/fluidmatrixinteractions/EclStone1Material.hpp b/opm/material/fluidmatrixinteractions/EclStone1Material.hpp index 3ce17cb0a..86087dfa6 100644 --- a/opm/material/fluidmatrixinteractions/EclStone1Material.hpp +++ b/opm/material/fluidmatrixinteractions/EclStone1Material.hpp @@ -146,6 +146,68 @@ public: Valgrind::CheckDefined(values[waterPhaseIdx]); } + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.oilWaterParams().pcSwMdc(); + krnSwMdc = params.oilWaterParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const double krwSw = 2.0; //Should not be used + params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.gasOilParams().pcSwMdc(); + krnSwMdc = params.gasOilParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const double krwSw = 2.0; //Should not be used + params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc); + } + /*! * \brief Capillary pressure between the gas and the non-wetting * liquid (i.e., oil) phase. diff --git a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp index 2eabd726b..a46915d05 100644 --- a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp +++ b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp @@ -147,6 +147,68 @@ public: Valgrind::CheckDefined(values[waterPhaseIdx]); } + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.oilWaterParams().pcSwMdc(); + krnSwMdc = params.oilWaterParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const double krwSw = 2.0; //Should not be used + params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.gasOilParams().pcSwMdc(); + krnSwMdc = params.gasOilParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const double krwSw = 2.0; //Should not be used + params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc); + } + /*! * \brief Capillary pressure between the gas and the non-wetting * liquid (i.e., oil) phase. diff --git a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp index 49017cbea..fe4a8a8df 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp @@ -163,6 +163,68 @@ public: } } + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.oilWaterParams().pcSwMdc(); + krnSwMdc = params.oilWaterParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for oil-water + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const Scalar krwSw = 2.0; //Should not be used + params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + const Params& params) + { + pcSwMdc = params.gasOilParams().pcSwMdc(); + krnSwMdc = params.gasOilParams().krnSwMdc(); + + Valgrind::CheckDefined(pcSwMdc); + Valgrind::CheckDefined(krnSwMdc); + } + + /* + * Hysteresis parameters for gas-oil + * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...) + * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...) + * \param params Parameters + */ + static void setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + Params& params) + { + const Scalar krwSw = 2.0; //Should not be used + params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc); + } + /*! * \brief Capillary pressure between the gas and the non-wetting * liquid (i.e., oil) phase. diff --git a/tests/test_eclmateriallawmanager.cpp b/tests/test_eclmateriallawmanager.cpp index b6b829810..4ef0c64e4 100644 --- a/tests/test_eclmateriallawmanager.cpp +++ b/tests/test_eclmateriallawmanager.cpp @@ -68,52 +68,52 @@ static const char* fam1DeckString = "GRID\n" "\n" "DX\n" - " 300*1000 /\n" + " 300*1000 /\n" "DY\n" - " 300*1000 /\n" + " 300*1000 /\n" "DZ\n" - " 100*20 100*30 100*50 /\n" + " 100*20 100*30 100*50 /\n" "\n" "TOPS\n" - " 100*8325 /\n" + " 100*8325 /\n" "\n" "\n" "PROPS\n" "\n" "SWOF\n" - "0.12 0 1 0\n" - "0.18 4.64876033057851E-008 1 0\n" - "0.24 0.000000186 0.997 0\n" - "0.3 4.18388429752066E-007 0.98 0\n" - "0.36 7.43801652892562E-007 0.7 0\n" - "0.42 1.16219008264463E-006 0.35 0\n" - "0.48 1.67355371900826E-006 0.2 0\n" - "0.54 2.27789256198347E-006 0.09 0\n" - "0.6 2.97520661157025E-006 0.021 0\n" - "0.66 3.7654958677686E-006 0.01 0\n" - "0.72 4.64876033057851E-006 0.001 0\n" - "0.78 0.000005625 0.0001 0\n" - "0.84 6.69421487603306E-006 0 0\n" - "0.91 8.05914256198347E-006 0 0\n" - "1 0.984 0 0 /\n" + "0.12 0 1 0\n" + "0.18 4.64876033057851E-008 1 0\n" + "0.24 0.000000186 0.997 0\n" + "0.3 4.18388429752066E-007 0.98 0\n" + "0.36 7.43801652892562E-007 0.7 0\n" + "0.42 1.16219008264463E-006 0.35 0\n" + "0.48 1.67355371900826E-006 0.2 0\n" + "0.54 2.27789256198347E-006 0.09 0\n" + "0.6 2.97520661157025E-006 0.021 0\n" + "0.66 3.7654958677686E-006 0.01 0\n" + "0.72 4.64876033057851E-006 0.001 0\n" + "0.78 0.000005625 0.0001 0\n" + "0.84 6.69421487603306E-006 0 0\n" + "0.91 8.05914256198347E-006 0 0\n" + "1 0.984 0 0 /\n" "\n" "\n" "SGOF\n" - "0 0 1 0\n" - "0.001 0 1 0\n" - "0.02 0 0.997 0\n" - "0.05 0.005 0.980 0\n" - "0.12 0.025 0.700 0\n" - "0.2 0.075 0.350 0\n" - "0.25 0.125 0.200 0\n" - "0.3 0.190 0.090 0\n" - "0.4 0.410 0.021 0\n" - "0.45 0.60 0.010 0\n" - "0.5 0.72 0.001 0\n" - "0.6 0.87 0.0001 0\n" - "0.7 0.94 0.000 0\n" - "0.85 0.98 0.000 0\n" - "0.88 0.984 0.000 0 /\n"; + "0 0 1 0\n" + "0.001 0 1 0\n" + "0.02 0 0.997 0\n" + "0.05 0.005 0.980 0\n" + "0.12 0.025 0.700 0\n" + "0.2 0.075 0.350 0\n" + "0.25 0.125 0.200 0\n" + "0.3 0.190 0.090 0\n" + "0.4 0.410 0.021 0\n" + "0.45 0.60 0.010 0\n" + "0.5 0.72 0.001 0\n" + "0.6 0.87 0.0001 0\n" + "0.7 0.94 0.000 0\n" + "0.85 0.98 0.000 0\n" + "0.88 0.984 0.000 0 /\n"; static const char* fam2DeckString = "RUNSPEC\n" @@ -135,62 +135,62 @@ static const char* fam2DeckString = "GRID\n" "\n" "DX\n" - " 300*1000 /\n" + " 300*1000 /\n" "DY\n" - " 300*1000 /\n" + " 300*1000 /\n" "DZ\n" - " 100*20 100*30 100*50 /\n" + " 100*20 100*30 100*50 /\n" "\n" "TOPS\n" - " 100*8325 /\n" + " 100*8325 /\n" "\n" "\n" "PROPS\n" "\n" "PVTW\n" - " 4017.55 1.038 3.22E-6 0.318 0.0 /\n" + " 4017.55 1.038 3.22E-6 0.318 0.0 /\n" "\n" "\n" "SWFN\n" - "0.12 0 0\n" - "0.18 4.64876033057851E-008 0\n" - "0.24 0.000000186 0\n" - "0.3 4.18388429752066E-007 0\n" - "0.36 7.43801652892562E-007 0\n" - "0.42 1.16219008264463E-006 0\n" - "0.48 1.67355371900826E-006 0\n" - "0.54 2.27789256198347E-006 0\n" - "0.6 2.97520661157025E-006 0\n" - "0.66 3.7654958677686E-006 0\n" - "0.72 4.64876033057851E-006 0\n" - "0.78 0.000005625 0\n" - "0.84 6.69421487603306E-006 0\n" - "0.91 8.05914256198347E-006 0\n" - "1 0.984 0 /\n" + "0.12 0 0\n" + "0.18 4.64876033057851E-008 0\n" + "0.24 0.000000186 0\n" + "0.3 4.18388429752066E-007 0\n" + "0.36 7.43801652892562E-007 0\n" + "0.42 1.16219008264463E-006 0\n" + "0.48 1.67355371900826E-006 0\n" + "0.54 2.27789256198347E-006 0\n" + "0.6 2.97520661157025E-006 0\n" + "0.66 3.7654958677686E-006 0\n" + "0.72 4.64876033057851E-006 0\n" + "0.78 0.000005625 0\n" + "0.84 6.69421487603306E-006 0\n" + "0.91 8.05914256198347E-006 0\n" + "1 0.984 0 /\n" "\n" "\n" "SGFN\n" - "0 0 0\n" - "0.001 0 0\n" - "0.02 0 0\n" - "0.05 0.005 0\n" - "0.12 0.025 0\n" - "0.2 0.075 0\n" - "0.25 0.125 0\n" - "0.3 0.190 0\n" - "0.4 0.410 0\n" - "0.45 0.60 0\n" - "0.5 0.72 0\n" - "0.6 0.87 0\n" - "0.7 0.94 0\n" - "0.85 0.98 0\n" - "0.88 0.984 0 /\n" + "0 0 0\n" + "0.001 0 0\n" + "0.02 0 0\n" + "0.05 0.005 0\n" + "0.12 0.025 0\n" + "0.2 0.075 0\n" + "0.25 0.125 0\n" + "0.3 0.190 0\n" + "0.4 0.410 0\n" + "0.45 0.60 0\n" + "0.5 0.72 0\n" + "0.6 0.87 0\n" + "0.7 0.94 0\n" + "0.85 0.98 0\n" + "0.88 0.984 0 /\n" "\n" "SOF3\n" " 0 0 0 \n" " 0.03 0 0 \n" " 0.09 0 0 \n" - " 0.16 0 0 \n" + " 0.16 0 0 \n" " 0.18 1* 0 \n" " 0.22 0.0001 1* \n" " 0.28 0.001 0.0001 \n" @@ -213,6 +213,82 @@ static const char* fam2DeckString = " 0.88 1 1 / \n" "\n"; +//Taken as a mix of the SPE1 cases above, and Norne to enable hysteresis +static const char* hysterDeckString = + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 10 10 3 /\n" + "\n" + "TABDIMS\n" + "/\n" + "\n" + "OIL\n" + "GAS\n" + "WATER\n" + "\n" + "DISGAS\n" + "\n" + "FIELD\n" + "\n" + "GRID\n" + "\n" + "DX\n" + " 300*1000 /\n" + "DY\n" + " 300*1000 /\n" + "DZ\n" + " 100*20 100*30 100*50 /\n" + "\n" + "TOPS\n" + " 100*8325 /\n" + "\n" + "\n" + "EHYSTR\n" + "0.1 0 0.1 1* KR /\n" + "\n" + "SATOPTS\n" + "HYSTER /\n" + "\n" + "PROPS\n" + "\n" + "SWOF\n" + "0.12 0 1 0\n" + "0.18 4.64876033057851E-008 1 0\n" + "0.24 0.000000186 0.997 0\n" + "0.3 4.18388429752066E-007 0.98 0\n" + "0.36 7.43801652892562E-007 0.7 0\n" + "0.42 1.16219008264463E-006 0.35 0\n" + "0.48 1.67355371900826E-006 0.2 0\n" + "0.54 2.27789256198347E-006 0.09 0\n" + "0.6 2.97520661157025E-006 0.021 0\n" + "0.66 3.7654958677686E-006 0.01 0\n" + "0.72 4.64876033057851E-006 0.001 0\n" + "0.78 0.000005625 0.0001 0\n" + "0.84 6.69421487603306E-006 0 0\n" + "0.91 8.05914256198347E-006 0 0\n" + "1 0.984 0 0 /\n" + "\n" + "\n" + "SGOF\n" + "0 0 1 0\n" + "0.001 0 1 0\n" + "0.02 0 0.997 0\n" + "0.05 0.005 0.980 0\n" + "0.12 0.025 0.700 0\n" + "0.2 0.075 0.350 0\n" + "0.25 0.125 0.200 0\n" + "0.3 0.190 0.090 0\n" + "0.4 0.410 0.021 0\n" + "0.45 0.60 0.010 0\n" + "0.5 0.72 0.001 0\n" + "0.6 0.87 0.0001 0\n" + "0.7 0.94 0.000 0\n" + "0.85 0.98 0.000 0\n" + "0.88 0.984 0.000 0 /\n"; + + + template inline void testAll() { @@ -280,8 +356,24 @@ inline void testAll() OPM_THROW(std::logic_error, "Discrepancy between the deck and the EclMaterialLawManager"); + const auto hysterDeck = parser.parseString(hysterDeckString, parseContext); + const Opm::EclipseState hysterEclState(hysterDeck, parseContext); + + Opm::EclMaterialLawManager hysterMaterialLawManager; + hysterMaterialLawManager.initFromDeck(hysterDeck, hysterEclState, compressedToCartesianIdx); + + if (hysterMaterialLawManager.enableEndPointScaling()) + OPM_THROW(std::logic_error, + "Discrepancy between the deck and the EclMaterialLawManager"); + + if (hysterMaterialLawManager.enableHysteresis() != true) + OPM_THROW(std::logic_error, + "Discrepancy between the deck and the EclMaterialLawManager"); + + + // make sure that the saturation functions for both keyword families are - // identical + // identical, and that setting and getting the hysteresis parameters works for (unsigned elemIdx = 0; elemIdx < n; ++ elemIdx) { for (int i = -10; i < 120; ++ i) { Scalar Sw = Scalar(i)/100; @@ -318,7 +410,43 @@ inline void testAll() "Discrepancy between capillary pressure of family 1 and family 2 keywords"); if (std::abs(krFam1[phaseIdx] - krFam2[phaseIdx]) > 1e-1) OPM_THROW(std::logic_error, - "Discrepancy between capillary pressure of family 1 and family 2 keywords"); + "Discrepancy between relative permeabilities of family 1 and family 2 keywords"); + } + + + // This should ideally test each of the materials (stone1, stone2, default, two-phase), + // but currently only tests default + const Scalar pcSwMdc_in[2] = { 1.0/2.0, 1.0/3.0 }; + const Scalar krnSwMdc_in[2] = { 1.0/5.0, 1.0/7.0 }; + hysterMaterialLawManager.setOilWaterHysteresisParams( + pcSwMdc_in[0], + krnSwMdc_in[0], + elemIdx); + hysterMaterialLawManager.setGasOilHysteresisParams( + pcSwMdc_in[1], + krnSwMdc_in[1], + elemIdx); + + Scalar pcSwMdc_out[2]; + Scalar krnSwMdc_out[2]; + Scalar deltaSwImbKrn_out[2]; + hysterMaterialLawManager.oilWaterHysteresisParams( + pcSwMdc_out[0], + krnSwMdc_out[0], + elemIdx); + hysterMaterialLawManager.gasOilHysteresisParams( + pcSwMdc_out[1], + krnSwMdc_out[1], + elemIdx); + + for (unsigned phasePairIdx = 0; phasePairIdx < 2; ++ phasePairIdx) { + if ((pcSwMdc_in[phasePairIdx] - pcSwMdc_out[phasePairIdx]) != 0.0) + OPM_THROW(std::logic_error, + "Hysteresis parameters did not propagate correctly"); + if ((krnSwMdc_in[phasePairIdx] - krnSwMdc_out[phasePairIdx]) != 0.0) + OPM_THROW(std::logic_error, + "Hysteresis parameters did not propagate correctly"); + } } }