diff --git a/opm/material/fluidmatrixinteractions/EclEpsConfig.hpp b/opm/material/fluidmatrixinteractions/EclEpsConfig.hpp index f9a404a88..b40ee8dc7 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsConfig.hpp +++ b/opm/material/fluidmatrixinteractions/EclEpsConfig.hpp @@ -33,6 +33,8 @@ #include #include #include +#include +#include #endif #include @@ -64,6 +66,7 @@ public: // by default, do not scale anything enableSatScaling_ = false; enablePcScaling_ = false; + enableLeverettScaling_ = false; enableKrwScaling_ = false; enableKrnScaling_ = false; enableThreePointKrSatScaling_ = false; @@ -131,6 +134,58 @@ public: bool enablePcScaling() const { return enablePcScaling_; } + /*! + * \brief Specify whether the Leverett capillary pressure scaling is enabled. + * + * If setting this to true, Leverett capillary pressure scaling will be used instead + * of the normal capillary pressure scaling and the value of enablePcScaling() will + * not matter anymore. + */ + void setEnableLeverettScaling(bool yesno) + { enableLeverettScaling_ = yesno; } + + /*! + * \brief Specify the conversion factor between deck and SI pressure units + * + * i.e. for METRIC decks this should be set to 10^5, for FIELD decks to 6894.7573, + * and so on. This factor is required to be set if Leverett pressure scaling is + * enabled because the ECL Leverett function implementation "misuses" the pressure + * column of the underlying tables (SWOF, SGOF, SWFN, etc.) for the (dimensionless) + * J-Function. Since opm-parser automatically applies the conversion factor for deck + * to SI pressure values to the column, opm-material needs this value in order to + * restore the value of the J-function from the "pressure" value. + * + * Note that this is a BIG *FAT* *_HACK_* + */ + void setDeckPressureConversionFactor(double pToSiFactor) + { pToSiFactor_ = pToSiFactor; } + + /*! + * \brief Returns the conversion factor between deck and SI pressure units + * + * i.e. for METRIC decks this should be set to 10^5, for FIELD decks to 6894.7573, + * and so on. This factor is required to be set if Leverett pressure scaling is + * enabled because the ECL Leverett function implementation "misuses" the pressure + * column of the underlying tables (SWOF, SGOF, SWFN, etc.) for the (dimensionless) + * J-Function. Since opm-parser automatically applies the conversion factor for deck + * to SI pressure values to the column, opm-material needs this value in order to + * restore the value of the J-function from the "pressure" value. + * + * Note that this is a BIG *FAT* *_HACK_* + */ + double deckPressureConversionFactor() const + { return pToSiFactor_; } + + /*! + * \brief Returns whether the Leverett capillary pressure scaling is enabled. + * + * If this returns true, Leverett capillary pressure scaling will be used instead of + * the normal capillary pressure scaling and the value of enablePcScaling() does not + * matter anymore. + */ + bool enableLeverettScaling() const + { return enableLeverettScaling_; } + #if HAVE_OPM_PARSER /*! * \brief Reads all relevant material parameters form a cell of a parsed ECL deck. @@ -148,6 +203,7 @@ public: enableSatScaling_ = false; enableThreePointKrSatScaling_ = false; enablePcScaling_ = false; + enableLeverettScaling_ = false; enableKrwScaling_ = false; enableKrnScaling_ = false; return; @@ -177,16 +233,59 @@ public: else enableThreePointKrSatScaling_ = false; + // determine the conversion factor for deck pressures to Pascals + pToSiFactor_ = deck.getActiveUnitSystem().getDimension("Pressure").getSIScaling(); + auto& props = eclState.get3DProperties(); // check if we are supposed to scale the Y axis of the capillary pressure - if (twoPhaseSystemType == EclOilWaterSystem) + if (twoPhaseSystemType == EclOilWaterSystem) { enablePcScaling_ = props.hasDeckDoubleGridProperty("PCW") || props.hasDeckDoubleGridProperty("SWATINIT"); + // check if Leverett capillary pressure scaling is requested + if (deck.hasKeyword("JFUNC")) { + const auto& jfuncKeyword = deck.getKeyword("JFUNC"); + std::string flagString = + jfuncKeyword.getRecord(0).getItem("FLAG").getTrimmedString(0); + std::transform(flagString.begin(), + flagString.end(), + flagString.begin(), + ::toupper); + + if (flagString == "BOTH" || flagString == "WATER") + enableLeverettScaling_ = true; + } + + if (enablePcScaling_ && enableLeverettScaling_) + OPM_THROW(std::runtime_error, + "Capillary pressure scaling and the Leverett scaling function are " + "mutually exclusive: The deck contains the PCW property and the " + "JFUNC keyword applies to the water phase."); + } else { assert(twoPhaseSystemType == EclGasOilSystem); enablePcScaling_ = props.hasDeckDoubleGridProperty("PCG"); + + // check if Leverett capillary pressure scaling is requested + if (deck.hasKeyword("JFUNC")) { + const auto& jfuncKeyword = deck.getKeyword("JFUNC"); + std::string flagString = + jfuncKeyword.getRecord(0).getItem("FLAG").getTrimmedString(0); + std::transform(flagString.begin(), + flagString.end(), + flagString.begin(), + ::toupper); + + if (flagString == "BOTH" || flagString == "GAS") + enableLeverettScaling_ = true; + } + + if (enablePcScaling_ && enableLeverettScaling_) + OPM_THROW(std::runtime_error, + "Capillary pressure scaling and the Leverett scaling function are " + "mutually exclusive: The deck contains the PCG property and the " + "JFUNC keyword applies to the gas phase."); } // check if we are supposed to scale the Y axis of the wetting phase relperm @@ -208,6 +307,9 @@ public: #endif private: + // conversion factor for deck pressure units to SI quantities + double pToSiFactor_; + // enable scaling of the input saturations (i.e., rescale the x-Axis) bool enableSatScaling_; @@ -221,6 +323,7 @@ private: // enable the scaling of the capillary pressure and relative permeability outputs // (i.e., rescale the y-Axis) bool enablePcScaling_; + bool enableLeverettScaling_; bool enableKrwScaling_; bool enableKrnScaling_; }; diff --git a/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp b/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp index 21ca822d5..06464e2c2 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp +++ b/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp @@ -43,6 +43,8 @@ #include #endif +#include + #include #include #include @@ -87,6 +89,18 @@ public: retrieveGridPropertyData_(&krw, eclState, kwPrefix+"KRW"); retrieveGridPropertyData_(&kro, eclState, kwPrefix+"KRO"); retrieveGridPropertyData_(&krg, eclState, kwPrefix+"KRG"); + + // _may_ be needed to calculate the Leverett capillary pressure scaling factor + const auto& ecl3dProps = eclState.get3DProperties(); + poro = &ecl3dProps.getDoubleGridProperty("PORO").getData(); + + if (ecl3dProps.hasDeckDoubleGridProperty("PERMX")) { + permx = &ecl3dProps.getDoubleGridProperty("PERMX").getData(); + permy = permx; + } + + if (ecl3dProps.hasDeckDoubleGridProperty("PERMY")) + permy = &ecl3dProps.getDoubleGridProperty("PERMY").getData(); } #endif @@ -105,6 +119,9 @@ public: const DoubleData* krw; const DoubleData* kro; const DoubleData* krg; + const DoubleData* poro; + const DoubleData* permx; + const DoubleData* permy; private: #if HAVE_OPM_PARSER @@ -153,6 +170,11 @@ struct EclEpsScalingPointsInfo Scalar maxPcow; // maximum capillary pressure of the oil-water system Scalar maxPcgo; // maximum capillary pressure of the gas-oil system + // the Leverett capillary pressure scaling factors. (those only make sense for the + // scaled points, for the unscaled ones they are 1.0.) + Scalar pcowLeverettFactor; + Scalar pcgoLeverettFactor; + // maximum relative permabilities Scalar maxKrw; // maximum relative permability of water Scalar maxKrow; // maximum relative permability of oil in the oil-water system @@ -175,6 +197,8 @@ struct EclEpsScalingPointsInfo << " Sogu: " << Sogu << "\n" << " maxPcow: " << maxPcow << "\n" << " maxPcgo: " << maxPcgo << "\n" + << " pcowLeverettFactor: " << pcowLeverettFactor << "\n" + << " pcgoLeverettFactor: " << pcgoLeverettFactor << "\n" << " maxKrw: " << maxKrw << "\n" << " maxKrg: " << maxKrg << "\n" << " maxKrow: " << maxKrow << "\n" @@ -257,15 +281,20 @@ struct EclEpsScalingPointsInfo throw std::domain_error("No valid saturation keyword family specified"); } + // there are no "unscaled" Leverett factors, so we just set them to 1.0 + pcowLeverettFactor = 1.0; + pcgoLeverettFactor = 1.0; } -#endif /*! * \brief Extract the values of the scaled scaling parameters. * * I.e., the values which are "seen" by the physical model. */ - void extractScaled(const EclEpsGridProperties& epsProperties, unsigned cartesianCellIdx) + void extractScaled(const Opm::Deck& deck, + const Opm::EclipseState& /*eclState*/, + const EclEpsGridProperties& epsProperties, + unsigned cartesianCellIdx) { // overwrite the unscaled values with the values for the cell if it is // explicitly specified by the corresponding keyword. @@ -287,7 +316,75 @@ struct EclEpsScalingPointsInfo // quite likely that's wrong! extractGridPropertyValue_(maxKrow, epsProperties.kro, cartesianCellIdx); extractGridPropertyValue_(maxKrog, epsProperties.kro, cartesianCellIdx); + + // compute the Leverett capillary pressure scaling factors if applicable. note + // that this needs to be done using non-SI units to make it correspond to the + // documentation. + pcowLeverettFactor = 1.0; + pcgoLeverettFactor = 1.0; + if (deck.hasKeyword("JFUNC")) { + const auto& jfuncKeyword = deck.getKeyword("JFUNC"); + const auto& jfuncRecord = jfuncKeyword.getRecord(0); + const std::string& jfuncFlag = jfuncRecord.getItem("FLAG").template get(0); + const std::string& jfuncDir = + jfuncRecord.getItem("DIRECTION").template get(0); + + Scalar perm; + if (jfuncDir == "X") + perm = + (*epsProperties.permx)[cartesianCellIdx]; + else if (jfuncDir == "Y") + perm = + (*epsProperties.permy)[cartesianCellIdx]; + else if (jfuncDir == "XY") + // TODO: verify that this really is the arithmetic mean. (the + // documentation just says that the "average" should be used, IMO the + // harmonic mean would be more appropriate because that's what's usually + // applied when calculating the fluxes.) + perm = + Opm::arithmeticMean((*epsProperties.permx)[cartesianCellIdx], + (*epsProperties.permy)[cartesianCellIdx]); + else + OPM_THROW(std::runtime_error, "Illegal direction indicator for the JFUNC " + "keyword ('"<(0); + Scalar beta = jfuncRecord.getItem("BETA_FACTOR").template get(0); + + // the part of the Leverett capillary pressure which does not depend on + // surface tension. + Scalar commonFactor = std::pow(poro, alpha)/std::pow(perm, beta); + + // multiply the documented constant by 10^5 because we want the pressures + // in [Pa], not in [bar] + const Scalar Uconst = 0.318316 * 1e5; + + // compute the oil-water Leverett factor. + if (jfuncFlag == "WATER" || jfuncFlag == "BOTH") { + // note that we use the surface tension in terms of [dyn/cm] + Scalar gamma = + jfuncRecord.getItem("OW_SURFACE_TENSION").getSIDouble(0) + * (/*dyn/N=*/1e5 * /*m/cm=*/1e-2); + + pcowLeverettFactor = commonFactor*gamma*Uconst; + } + + // compute the gas-oil Leverett factor. + if (jfuncFlag == "GAS" || jfuncFlag == "BOTH") { + // note that we use the surface tension in terms of [dyn/cm] + Scalar gamma = + jfuncRecord.getItem("GO_SURFACE_TENSION").getSIDouble(0) + * (/*dyn/N=*/1e5 * /*m/cm=*/1e-2); + + pcgoLeverettFactor = commonFactor*gamma*Uconst; + } + } } +#endif private: #if HAVE_OPM_PARSER @@ -538,7 +635,10 @@ public: saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl; } - maxPcnw_ = epsInfo.maxPcow; + if (config.enableLeverettScaling()) + maxPcnwOrLeverettFactor_ = epsInfo.pcowLeverettFactor; + else + maxPcnwOrLeverettFactor_ = epsInfo.maxPcow; maxKrw_ = epsInfo.maxKrw; maxKrn_ = epsInfo.maxKrow; } @@ -574,7 +674,11 @@ public: saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu; } - maxPcnw_ = epsInfo.maxPcgo; + if (config.enableLeverettScaling()) + maxPcnwOrLeverettFactor_ = epsInfo.pcgoLeverettFactor; + else + maxPcnwOrLeverettFactor_ = epsInfo.maxPcgo; + maxKrw_ = epsInfo.maxKrog; maxKrn_ = epsInfo.maxKrg; } @@ -620,13 +724,25 @@ public: * \brief Sets the maximum capillary pressure */ void setMaxPcnw(Scalar value) - { maxPcnw_ = value; } + { maxPcnwOrLeverettFactor_ = value; } /*! * \brief Returns the maximum capillary pressure */ Scalar maxPcnw() const - { return maxPcnw_; } + { return maxPcnwOrLeverettFactor_; } + + /*! + * \brief Sets the Leverett scaling factor for capillary pressure + */ + void setLeverettFactor(Scalar value) + { maxPcnwOrLeverettFactor_ = value; } + + /*! + * \brief Returns the Leverett scaling factor for capillary pressure + */ + Scalar leverettFactor() const + { return maxPcnwOrLeverettFactor_; } /*! * \brief Sets the maximum wetting phase relative permeability @@ -661,7 +777,7 @@ public: private: // The the points used for the "y-axis" scaling of capillary pressure - Scalar maxPcnw_; + Scalar maxPcnwOrLeverettFactor_; // The the points used for the "y-axis" scaling of wetting phase relative permability Scalar maxKrw_; diff --git a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp index 599268c24..a00f72e93 100644 --- a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp +++ b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp @@ -468,21 +468,40 @@ private: template static Evaluation unscaledToScaledPcnw_(const Params ¶ms, const Evaluation& unscaledPcnw) { - if (!params.config().enablePcScaling()) - return unscaledPcnw; + if (params.config().enableLeverettScaling()) { + // HACK: we need to divide the "unscaled capillary pressure" by the + // conversion factor for deck to SI pressures because this quantity is + // actually *not* a pressure but the (dimensionless) value of the + // J-function. Since opm-parser currently always treads this column as a + // pressure, it multiplies it with the conversion factor and we have to + // revert this here. + Scalar alpha = params.scaledPoints().leverettFactor(); + return unscaledPcnw*(alpha/params.config().deckPressureConversionFactor()); + } + else if (params.config().enablePcScaling()) { + Scalar alpha = params.scaledPoints().maxPcnw()/params.unscaledPoints().maxPcnw(); + return unscaledPcnw*alpha; + } - Scalar alpha = params.scaledPoints().maxPcnw()/params.unscaledPoints().maxPcnw(); - return unscaledPcnw*alpha; + return unscaledPcnw; } template static Evaluation scaledToUnscaledPcnw_(const Params ¶ms, const Evaluation& scaledPcnw) { - if (!params.config().enablePcScaling()) - return scaledPcnw; + if (params.config().enableLeverettScaling()) { + // HACK: we need to multiply the "scaled capillary pressure" by the conversion + // factor for deck to SI pressures because this quantity is actually *not* a + // pressure but the (dimensionless) value of the J-function. + Scalar alpha = params.scaledPoints().leverettFactor(); + return scaledPcnw*(params.config().deckPressureConversionFactor()/alpha); + } + else if (params.config().enablePcScaling()) { + Scalar alpha = params.unscaledPoints().maxPcnw()/params.scaledPoints().maxPcnw(); + return scaledPcnw/alpha; + } - Scalar alpha = params.unscaledPoints().maxPcnw()/params.scaledPoints().maxPcnw(); - return scaledPcnw/alpha; + return scaledPcnw; } /*! diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index a6bf2284a..9ece64c20 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp @@ -274,7 +274,6 @@ private: oilWaterEclEpsConfig_-> initFromDeck(deck, eclState, Opm::EclOilWaterSystem); enableEndPointScaling_ = deck.hasKeyword("ENDSCALE"); - } void readGlobalHysteresisOptions_(const Opm::Deck& deck) @@ -379,12 +378,16 @@ private: readGasOilScaledPoints_(gasOilScaledInfoVector, gasOilScaledPointsVector, gasOilConfig, + deck, + eclState, epsGridProperties, elemIdx, cartElemIdx); readOilWaterScaledPoints_(oilWaterScaledEpsInfoDrainage_, oilWaterScaledEpsPointsDrainage, oilWaterConfig, + deck, + eclState, epsGridProperties, elemIdx, cartElemIdx); @@ -393,12 +396,16 @@ private: readGasOilScaledPoints_(gasOilScaledImbInfoVector, gasOilScaledImbPointsVector, gasOilConfig, + deck, + eclState, epsImbGridProperties, elemIdx, cartElemIdx); readOilWaterScaledPoints_(oilWaterScaledImbInfoVector, oilWaterScaledImbPointsVector, oilWaterConfig, + deck, + eclState, epsImbGridProperties, elemIdx, cartElemIdx); @@ -771,6 +778,8 @@ private: void readGasOilScaledPoints_(InfoContainer& destInfo, PointsContainer& destPoints, std::shared_ptr config, + const Opm::Deck& deck, + const Opm::EclipseState& eclState, const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, unsigned cartElemIdx) @@ -778,7 +787,7 @@ private: unsigned satnumIdx = static_cast((*epsGridProperties.satnum)[cartElemIdx]) - 1; // ECL uses Fortran indices! destInfo[elemIdx] = std::make_shared >(unscaledEpsInfo_[satnumIdx]); - destInfo[elemIdx]->extractScaled(epsGridProperties, cartElemIdx); + destInfo[elemIdx]->extractScaled(deck, eclState, epsGridProperties, cartElemIdx); destPoints[elemIdx] = std::make_shared >(); destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclGasOilSystem); @@ -788,6 +797,8 @@ private: void readOilWaterScaledPoints_(InfoContainer& destInfo, PointsContainer& destPoints, std::shared_ptr config, + const Opm::Deck& deck, + const Opm::EclipseState& eclState, const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, unsigned cartElemIdx) @@ -795,7 +806,7 @@ private: unsigned satnumIdx = static_cast((*epsGridProperties.satnum)[cartElemIdx]) - 1; // ECL uses Fortran indices! destInfo[elemIdx] = std::make_shared >(unscaledEpsInfo_[satnumIdx]); - destInfo[elemIdx]->extractScaled(epsGridProperties, cartElemIdx); + destInfo[elemIdx]->extractScaled(deck, eclState, epsGridProperties, cartElemIdx); destPoints[elemIdx] = std::make_shared >(); destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclOilWaterSystem);