diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index aaeb4bdab7c..540dcbad2cf 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -235,6 +235,7 @@ if(ENABLE_ECL_INPUT) src/opm/material/fluidmatrixinteractions/EclEpsConfig.cpp src/opm/material/fluidmatrixinteractions/EclEpsGridProperties.cpp src/opm/material/fluidmatrixinteractions/EclHysteresisConfig.cpp + src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp ) diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp index e0583673315..b936ee140e2 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp @@ -31,41 +31,34 @@ #ifndef OPM_ECL_MATERIAL_LAW_MANAGER_HPP #define OPM_ECL_MATERIAL_LAW_MANAGER_HPP +#include + #include -#include #include #include -#include -#include -#include -#include #include #include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include #include #include -#include +#include #include namespace Opm { +class EclipseState; +class EclEpsConfig; +template class EclEpsScalingPoints; +template struct EclEpsScalingPointsInfo; +class EclHysteresisConfig; +enum class EclTwoPhaseSystemType; +class Runspec; +class SgfnTable; +class SgofTable; +class SlgofTable; +class Sof2Table; +class Sof3Table; + /*! * \ingroup fluidmatrixinteractions * @@ -117,6 +110,9 @@ class EclMaterialLawManager using MaterialLaw = EclMultiplexerMaterial; using MaterialLawParams = typename MaterialLaw::Params; + EclMaterialLawManager(); + ~EclMaterialLawManager(); + private: // internal typedefs using GasOilEffectiveParamVector = std::vector>; @@ -133,298 +129,9 @@ class EclMaterialLawManager using MaterialLawParamsVector = std::vector>; public: - EclMaterialLawManager() - {} - - void initFromState(const EclipseState& eclState) - { - // get the number of saturation regions and the number of cells in the deck - const auto& runspec = eclState.runspec(); - const size_t numSatRegions = runspec.tabdims().getNumSatTables(); - - const auto& ph = runspec.phases(); - this->hasGas = ph.active(Phase::GAS); - this->hasOil = ph.active(Phase::OIL); - this->hasWater = ph.active(Phase::WATER); - - readGlobalEpsOptions_(eclState); - readGlobalHysteresisOptions_(eclState); - readGlobalThreePhaseOptions_(runspec); - - // Read the end point scaling configuration (once per run). - gasOilConfig = std::make_shared(); - oilWaterConfig = std::make_shared(); - gasWaterConfig = std::make_shared(); - gasOilConfig->initFromState(eclState, EclTwoPhaseSystemType::GasOil); - oilWaterConfig->initFromState(eclState, EclTwoPhaseSystemType::OilWater); - gasWaterConfig->initFromState(eclState, EclTwoPhaseSystemType::GasWater); - - - const auto& tables = eclState.getTableManager(); - - { - const auto& stone1exTables = tables.getStone1exTable(); - - if (! stone1exTables.empty()) { - stoneEtas.clear(); - stoneEtas.reserve(numSatRegions); - - for (const auto& table : stone1exTables) { - stoneEtas.push_back(table.eta); - } - } - } - - this->unscaledEpsInfo_.resize(numSatRegions); - - if (this->hasGas + this->hasOil + this->hasWater == 1) { - // Single-phase simulation. Special case. Nothing to do here. - return; - } - - // Multiphase simulation. Common case. - const auto tolcrit = runspec.saturationFunctionControls() - .minimumRelpermMobilityThreshold(); - - const auto rtep = satfunc::getRawTableEndpoints(tables, ph, tolcrit); - const auto rfunc = satfunc::getRawFunctionValues(tables, ph, rtep); - - for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { - this->unscaledEpsInfo_[satRegionIdx] - .extractUnscaled(rtep, rfunc, satRegionIdx); - } - } - - void initParamsForElements(const EclipseState& eclState, size_t numCompressedElems) - { - // get the number of saturation regions - const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables(); - - // setup the saturation region specific parameters - gasOilUnscaledPointsVector_.resize(numSatRegions); - oilWaterUnscaledPointsVector_.resize(numSatRegions); - gasWaterUnscaledPointsVector_.resize(numSatRegions); - - gasOilEffectiveParamVector_.resize(numSatRegions); - oilWaterEffectiveParamVector_.resize(numSatRegions); - gasWaterEffectiveParamVector_.resize(numSatRegions); - for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { - // unscaled points for end-point scaling - readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, eclState, satRegionIdx); - readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, eclState, satRegionIdx); - readGasWaterUnscaledPoints_(gasWaterUnscaledPointsVector_, gasWaterConfig, eclState, satRegionIdx); - - // the parameters for the effective two-phase matererial laws - readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, eclState, satRegionIdx); - readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, eclState, satRegionIdx); - readGasWaterEffectiveParameters_(gasWaterEffectiveParamVector_, eclState, satRegionIdx); - } - - // copy the SATNUM grid property. in some cases this is not necessary, but it - // should not require much memory anyway... - satnumRegionArray_.resize(numCompressedElems); - if (eclState.fieldProps().has_int("SATNUM")) { - const auto& satnumRawData = eclState.fieldProps().get_int("SATNUM"); - for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { - satnumRegionArray_[elemIdx] = satnumRawData[elemIdx] - 1; - } - } - else { - std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0); - } - auto copy_krnum = [&eclState, numCompressedElems](std::vector& dest, const std::string keyword) { - if (eclState.fieldProps().has_int(keyword)) { - dest.resize(numCompressedElems); - const auto& satnumRawData = eclState.fieldProps().get_int(keyword); - for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { - dest[elemIdx] = satnumRawData[elemIdx] - 1; - } - } - }; - copy_krnum(krnumXArray_, "KRNUMX"); - copy_krnum(krnumYArray_, "KRNUMY"); - copy_krnum(krnumZArray_, "KRNUMZ"); - - // create the information for the imbibition region (IMBNUM). By default this is - // the same as the saturation region (SATNUM) - imbnumRegionArray_ = satnumRegionArray_; - if (eclState.fieldProps().has_int("IMBNUM")) { - const auto& imbnumRawData = eclState.fieldProps().get_int("IMBNUM"); - for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { - imbnumRegionArray_[elemIdx] = imbnumRawData[elemIdx] - 1; - } - } - - assert(numCompressedElems == satnumRegionArray_.size()); - assert(!enableHysteresis() || numCompressedElems == imbnumRegionArray_.size()); - - // read the scaled end point scaling parameters which are specific for each - // element - oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems); - - std::unique_ptr epsImbGridProperties; - - if (enableHysteresis()) { - epsImbGridProperties = std::make_unique(eclState, true); - } - - EclEpsGridProperties epsGridProperties(eclState, false); - materialLawParams_.resize(numCompressedElems); - - for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { - unsigned satRegionIdx = static_cast(satnumRegionArray_[elemIdx]); - auto gasOilParams = std::make_shared(); - auto oilWaterParams = std::make_shared(); - auto gasWaterParams = std::make_shared(); - gasOilParams->setConfig(hysteresisConfig_); - oilWaterParams->setConfig(hysteresisConfig_); - gasWaterParams->setConfig(hysteresisConfig_); - - auto [gasOilScaledInfo, gasOilScaledPoint] = - readScaledPoints_(*gasOilConfig, - eclState, - epsGridProperties, - elemIdx, - EclTwoPhaseSystemType::GasOil); - - auto [owinfo, oilWaterScaledEpsPointDrainage] = - readScaledPoints_(*oilWaterConfig, - eclState, - epsGridProperties, - elemIdx, - EclTwoPhaseSystemType::OilWater); - oilWaterScaledEpsInfoDrainage_[elemIdx] = owinfo; - - auto [gasWaterScaledInfo, gasWaterScaledPoint] = - readScaledPoints_(*gasWaterConfig, - eclState, - epsGridProperties, - elemIdx, - EclTwoPhaseSystemType::GasWater); - - if (hasGas && hasOil) { - GasOilEpsTwoPhaseParams gasOilDrainParams; - gasOilDrainParams.setConfig(gasOilConfig); - gasOilDrainParams.setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - gasOilDrainParams.setScaledPoints(gasOilScaledPoint); - gasOilDrainParams.setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); - gasOilDrainParams.finalize(); - - gasOilParams->setDrainageParams(gasOilDrainParams, - gasOilScaledInfo, - EclTwoPhaseSystemType::GasOil); - } - - if (hasOil && hasWater) { - OilWaterEpsTwoPhaseParams oilWaterDrainParams; - oilWaterDrainParams.setConfig(oilWaterConfig); - oilWaterDrainParams.setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - oilWaterDrainParams.setScaledPoints(oilWaterScaledEpsPointDrainage); - oilWaterDrainParams.setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); - oilWaterDrainParams.finalize(); - - oilWaterParams->setDrainageParams(oilWaterDrainParams, - owinfo, - EclTwoPhaseSystemType::OilWater); - } - - if (hasGas && hasWater && !hasOil) { - GasWaterEpsTwoPhaseParams gasWaterDrainParams; - gasWaterDrainParams.setConfig(gasWaterConfig); - gasWaterDrainParams.setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]); - gasWaterDrainParams.setScaledPoints(gasWaterScaledPoint); - gasWaterDrainParams.setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]); - gasWaterDrainParams.finalize(); - - gasWaterParams->setDrainageParams(gasWaterDrainParams, - gasWaterScaledInfo, - EclTwoPhaseSystemType::GasWater); - } - - if (enableHysteresis()) { - auto [gasOilScaledImbInfo, gasOilScaledImbPoint] = - readScaledPoints_(*gasOilConfig, - eclState, - *epsImbGridProperties, - elemIdx, - EclTwoPhaseSystemType::GasOil); - - auto [oilWaterScaledImbInfo, oilWaterScaledImbPoint] = - readScaledPoints_(*oilWaterConfig, - eclState, - *epsImbGridProperties, - elemIdx, - EclTwoPhaseSystemType::OilWater); - - auto [gasWaterScaledImbInfo, gasWaterScaledImbPoint] = - readScaledPoints_(*gasWaterConfig, - eclState, - *epsImbGridProperties, - elemIdx, - EclTwoPhaseSystemType::GasWater); - - unsigned imbRegionIdx = imbnumRegionArray_[elemIdx]; - if (hasGas && hasOil) { - GasOilEpsTwoPhaseParams gasOilImbParamsHyst; - gasOilImbParamsHyst.setConfig(gasOilConfig); - gasOilImbParamsHyst.setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]); - gasOilImbParamsHyst.setScaledPoints(gasOilScaledImbPoint); - gasOilImbParamsHyst.setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]); - gasOilImbParamsHyst.finalize(); - - gasOilParams->setImbibitionParams(gasOilImbParamsHyst, - gasOilScaledImbInfo, - EclTwoPhaseSystemType::GasOil); - } - - if (hasOil && hasWater) { - OilWaterEpsTwoPhaseParams oilWaterImbParamsHyst; - oilWaterImbParamsHyst.setConfig(oilWaterConfig); - oilWaterImbParamsHyst.setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]); - oilWaterImbParamsHyst.setScaledPoints(oilWaterScaledImbPoint); - oilWaterImbParamsHyst.setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]); - oilWaterImbParamsHyst.finalize(); - - oilWaterParams->setImbibitionParams(oilWaterImbParamsHyst, - oilWaterScaledImbInfo, - EclTwoPhaseSystemType::OilWater); - } - - if (hasGas && hasWater && !hasOil) { - GasWaterEpsTwoPhaseParams gasWaterImbParamsHyst; - gasWaterImbParamsHyst.setConfig(gasWaterConfig); - gasWaterImbParamsHyst.setUnscaledPoints(gasWaterUnscaledPointsVector_[imbRegionIdx]); - gasWaterImbParamsHyst.setScaledPoints(gasWaterScaledImbPoint); - gasWaterImbParamsHyst.setEffectiveLawParams(gasWaterEffectiveParamVector_[imbRegionIdx]); - gasWaterImbParamsHyst.finalize(); - - gasWaterParams->setImbibitionParams(gasWaterImbParamsHyst, - gasWaterScaledImbInfo, - EclTwoPhaseSystemType::GasWater); - } - } - - if (hasGas && hasOil) - gasOilParams->finalize(); - - if (hasOil && hasWater) - oilWaterParams->finalize(); - - if (hasGas && hasWater && !hasOil) - gasWaterParams->finalize(); - - initThreePhaseParams_(eclState, - materialLawParams_[elemIdx], - satRegionIdx, - oilWaterScaledEpsInfoDrainage_[elemIdx], - oilWaterParams, - gasOilParams, - gasWaterParams); - - materialLawParams_[elemIdx].finalize(); - } - } + void initFromState(const EclipseState& eclState); + void initParamsForElements(const EclipseState& eclState, size_t numCompressedElems); /*! * \brief Modify the initial condition according to the SWATINIT keyword. @@ -436,53 +143,7 @@ class EclMaterialLawManager */ Scalar applySwatinit(unsigned elemIdx, Scalar pcow, - Scalar Sw) - { - auto& elemScaledEpsInfo = oilWaterScaledEpsInfoDrainage_[elemIdx]; - - // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74 - - if (pcow < 0.0) - Sw = elemScaledEpsInfo.Swu; - else { - - if (Sw <= elemScaledEpsInfo.Swl) - Sw = elemScaledEpsInfo.Swl; - - // specify a fluid state which only stores the saturations - using FluidState = SimpleModularFluidState don't care */ - /*storePressure=*/false, - /*storeTemperature=*/false, - /*storeComposition=*/false, - /*storeFugacity=*/false, - /*storeSaturation=*/true, - /*storeDensity=*/false, - /*storeViscosity=*/false, - /*storeEnthalpy=*/false>; - FluidState fs; - fs.setSaturation(waterPhaseIdx, Sw); - fs.setSaturation(gasPhaseIdx, 0); - fs.setSaturation(oilPhaseIdx, 0); - std::array pc = { 0 }; - MaterialLaw::capillaryPressures(pc, materialLawParams(elemIdx), fs); - - Scalar pcowAtSw = pc[oilPhaseIdx] - pc[waterPhaseIdx]; - constexpr const Scalar pcowAtSwThreshold = 1.0; //Pascal - // avoid divison by very small number - if (std::abs(pcowAtSw) > pcowAtSwThreshold) { - elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw; - auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx); - elemEclEpsScalingPoints.init(elemScaledEpsInfo, - *oilWaterEclEpsConfig_, - EclTwoPhaseSystemType::OilWater); - } - } - - return Sw; - } + Scalar Sw); bool enableEndPointScaling() const { return enableEndPointScaling_; } @@ -510,126 +171,21 @@ class EclMaterialLawManager * wells with its own saturation table idx. In order to reset the saturation table idx * in the materialLawparams_ call the method with the cells satRegionIdx */ - const MaterialLawParams& connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const - { - MaterialLawParams& mlp = const_cast(materialLawParams_[elemIdx]); - - if (enableHysteresis()) - OpmLog::warning("Warning: Using non-default satnum regions for connection is not tested in combination with hysteresis"); - // Currently we don't support COMPIMP. I.e. use the same table lookup for the hysteresis curves. - // unsigned impRegionIdx = satRegionIdx; - - // change the sat table it points to. - switch (mlp.approach()) { - case EclMultiplexerApproach::Stone1: { - auto& realParams = mlp.template getRealParams(); - - realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); -// if (enableHysteresis()) { -// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); -// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); -// } - } - break; - - case EclMultiplexerApproach::Stone2: { - auto& realParams = mlp.template getRealParams(); - realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); -// if (enableHysteresis()) { -// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); -// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); -// } - } - break; - - case EclMultiplexerApproach::Default: { - auto& realParams = mlp.template getRealParams(); - realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); -// if (enableHysteresis()) { -// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); -// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); -// } - } - break; - - case EclMultiplexerApproach::TwoPhase: { - auto& realParams = mlp.template getRealParams(); - realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); -// if (enableHysteresis()) { -// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); -// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); -// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); -// } - } - break; - - default: - throw std::logic_error("Enum value for material approach unknown!"); - } - - return mlp; - } + const MaterialLawParams& connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const; int satnumRegionIdx(unsigned elemIdx) const { return satnumRegionArray_[elemIdx]; } - int getKrnumSatIdx(unsigned elemIdx, FaceDir::DirEnum facedir) const { - using Dir = FaceDir::DirEnum; - const std::vector* array = nullptr; - switch(facedir) { - case Dir::XPlus: - array = &krnumXArray_; - break; - case Dir::YPlus: - array = &krnumYArray_; - break; - case Dir::ZPlus: - array = &krnumZArray_; - break; - default: - throw std::runtime_error("Unknown face direction"); - } - if (array->size() > 0) { - return (*array)[elemIdx]; - } - else { - return satnumRegionArray_[elemIdx]; - } - } - bool hasDirectionalRelperms() const { - if (krnumXArray_.size() > 0 || krnumYArray_.size() > 0 || krnumZArray_.size() > 0) { - return true; - } - return false; - } - int imbnumRegionIdx(unsigned elemIdx) const - { return imbnumRegionArray_[elemIdx]; } + int getKrnumSatIdx(unsigned elemIdx, FaceDir::DirEnum facedir) const; - std::shared_ptr& materialLawParamsPointerReferenceHack(unsigned elemIdx) + bool hasDirectionalRelperms() const { - assert(0 <= elemIdx && elemIdx < materialLawParams_.size()); - return materialLawParams_[elemIdx]; + return !krnumXArray_.empty() || !krnumYArray_.empty() || !krnumZArray_.empty(); } + int imbnumRegionIdx(unsigned elemIdx) const + { return imbnumRegionArray_[elemIdx]; } + template void updateHysteresis(const FluidState& fluidState, unsigned elemIdx) { @@ -641,519 +197,86 @@ class EclMaterialLawManager void oilWaterHysteresisParams(Scalar& pcSwMdc, Scalar& krnSwMdc, - unsigned elemIdx) const - { - if (!enableHysteresis()) - throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled."); - - const auto& params = materialLawParams(elemIdx); - MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params); - } + unsigned elemIdx) const; void setOilWaterHysteresisParams(const Scalar& pcSwMdc, const Scalar& krnSwMdc, - unsigned elemIdx) - { - if (!enableHysteresis()) - throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled."); - - auto& params = materialLawParams(elemIdx); - MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params); - } + unsigned elemIdx); void gasOilHysteresisParams(Scalar& pcSwMdc, Scalar& krnSwMdc, - unsigned elemIdx) const - { - if (!enableHysteresis()) - throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled."); - - const auto& params = materialLawParams(elemIdx); - MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params); - } + unsigned elemIdx) const; void setGasOilHysteresisParams(const Scalar& pcSwMdc, const Scalar& krnSwMdc, - unsigned elemIdx) - { - if (!enableHysteresis()) - throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled."); + unsigned elemIdx); - auto& params = materialLawParams(elemIdx); - MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params); - } - - EclEpsScalingPoints& oilWaterScaledEpsPointsDrainage(unsigned elemIdx) - { - auto& materialParams = materialLawParams_[elemIdx]; - switch (materialParams.approach()) { - case EclMultiplexerApproach::Stone1: { - auto& realParams = materialParams.template getRealParams(); - return realParams.oilWaterParams().drainageParams().scaledPoints(); - } - - case EclMultiplexerApproach::Stone2: { - auto& realParams = materialParams.template getRealParams(); - return realParams.oilWaterParams().drainageParams().scaledPoints(); - } - - case EclMultiplexerApproach::Default: { - auto& realParams = materialParams.template getRealParams(); - return realParams.oilWaterParams().drainageParams().scaledPoints(); - } - - case EclMultiplexerApproach::TwoPhase: { - auto& realParams = materialParams.template getRealParams(); - return realParams.oilWaterParams().drainageParams().scaledPoints(); - } - default: - throw std::logic_error("Enum value for material approach unknown!"); - } - } + EclEpsScalingPoints& oilWaterScaledEpsPointsDrainage(unsigned elemIdx); const EclEpsScalingPointsInfo& oilWaterScaledEpsInfoDrainage(size_t elemIdx) const { return oilWaterScaledEpsInfoDrainage_[elemIdx]; } private: - void readGlobalEpsOptions_(const EclipseState& eclState) - { - oilWaterEclEpsConfig_ = std::make_shared(); - oilWaterEclEpsConfig_->initFromState(eclState, EclTwoPhaseSystemType::OilWater); + void readGlobalEpsOptions_(const EclipseState& eclState); - enableEndPointScaling_ = eclState.getTableManager().hasTables("ENKRVD"); - } + void readGlobalHysteresisOptions_(const EclipseState& state); - void readGlobalHysteresisOptions_(const EclipseState& state) - { - hysteresisConfig_ = std::make_shared(); - hysteresisConfig_->initFromState(state.runspec()); - } - - void readGlobalThreePhaseOptions_(const Runspec& runspec) - { - bool gasEnabled = runspec.phases().active(Phase::GAS); - bool oilEnabled = runspec.phases().active(Phase::OIL); - bool waterEnabled = runspec.phases().active(Phase::WATER); - - int numEnabled = - (gasEnabled?1:0) - + (oilEnabled?1:0) - + (waterEnabled?1:0); - - if (numEnabled == 0) { - throw std::runtime_error("At least one fluid phase must be enabled. (Is: "+std::to_string(numEnabled)+")"); - } else if (numEnabled == 1) { - threePhaseApproach_ = EclMultiplexerApproach::OnePhase; - } else if ( numEnabled == 2) { - threePhaseApproach_ = EclMultiplexerApproach::TwoPhase; - if (!gasEnabled) - twoPhaseApproach_ = EclTwoPhaseApproach::OilWater; - else if (!oilEnabled) - twoPhaseApproach_ = EclTwoPhaseApproach::GasWater; - else if (!waterEnabled) - twoPhaseApproach_ = EclTwoPhaseApproach::GasOil; - } - else { - assert(numEnabled == 3); - - threePhaseApproach_ = EclMultiplexerApproach::Default; - const auto& satctrls = runspec.saturationFunctionControls(); - if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone2) - threePhaseApproach_ = EclMultiplexerApproach::Stone2; - else if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone1) - threePhaseApproach_ = EclMultiplexerApproach::Stone1; - } - } + void readGlobalThreePhaseOptions_(const Runspec& runspec); template void readGasOilEffectiveParameters_(Container& dest, const EclipseState& eclState, - unsigned satRegionIdx) - { - if (!hasGas || !hasOil) - // we don't read anything if either the gas or the oil phase is not active - return; - - dest[satRegionIdx] = std::make_shared(); - - auto& effParams = *dest[satRegionIdx]; - - // the situation for the gas phase is complicated that all saturations are - // shifted by the connate water saturation. - const Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl; - const auto tolcrit = eclState.runspec().saturationFunctionControls() - .minimumRelpermMobilityThreshold(); - - const auto& tableManager = eclState.getTableManager(); - - switch (eclState.runspec().saturationFunctionControls().family()) { - case SatFuncControls::KeywordFamily::Family_I: - { - const TableContainer& sgofTables = tableManager.getSgofTables(); - const TableContainer& slgofTables = tableManager.getSlgofTables(); - if (!sgofTables.empty()) - readGasOilEffectiveParametersSgof_(effParams, Swco, tolcrit, - sgofTables.getTable(satRegionIdx)); - else if (!slgofTables.empty()) - readGasOilEffectiveParametersSlgof_(effParams, Swco, tolcrit, - slgofTables.getTable(satRegionIdx)); - else if ( !tableManager.getSgofletTable().empty() ) { - const auto& letSgofTab = tableManager.getSgofletTable()[satRegionIdx]; - const std::vector dum; // dummy arg to comform with existing interface - - effParams.setApproach(SatCurveMultiplexerApproach::LET); - auto& realParams = effParams.template getRealParams(); - - // S=(So-Sogcr)/(1-Sogcr-Sgcr-Swco), krog = Krt*S^L/[S^L+E*(1.0-S)^T] - const Scalar s_min_w = letSgofTab.s2_critical; - const Scalar s_max_w = 1.0-letSgofTab.s1_critical-Swco; - const std::vector& letCoeffsOil = {s_min_w, s_max_w, - static_cast(letSgofTab.l2_relperm), - static_cast(letSgofTab.e2_relperm), - static_cast(letSgofTab.t2_relperm), - static_cast(letSgofTab.krt2_relperm)}; - realParams.setKrwSamples(letCoeffsOil, dum); - - // S=(1-So-Sgcr-Swco)/(1-Sogcr-Sgcr-Swco), krg = Krt*S^L/[S^L+E*(1.0-S)^T] - const Scalar s_min_nw = letSgofTab.s1_critical+Swco; - const Scalar s_max_nw = 1.0-letSgofTab.s2_critical; - const std::vector& letCoeffsGas = {s_min_nw, s_max_nw, - static_cast(letSgofTab.l1_relperm), - static_cast(letSgofTab.e1_relperm), - static_cast(letSgofTab.t1_relperm), - static_cast(letSgofTab.krt1_relperm)}; - realParams.setKrnSamples(letCoeffsGas, dum); - - // S=(So-Sorg)/(1-Sorg-Sgl-Swco), Pc = Pct + (pcir_pc-Pct)*(1-S)^L/[(1-S)^L+E*S^T] - const std::vector& letCoeffsPc = {static_cast(letSgofTab.s2_residual), - static_cast(letSgofTab.s1_residual+Swco), - static_cast(letSgofTab.l_pc), - static_cast(letSgofTab.e_pc), - static_cast(letSgofTab.t_pc), - static_cast(letSgofTab.pcir_pc), - static_cast(letSgofTab.pct_pc)}; - realParams.setPcnwSamples(letCoeffsPc, dum); - - realParams.finalize(); - } - break; - } - - case SatFuncControls::KeywordFamily::Family_II: - { - const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable( satRegionIdx ); - if (!hasWater) { - // oil and gas case - const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable( satRegionIdx ); - readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable); - } - else { - const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable( satRegionIdx ); - readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable); - } - break; - } - - case SatFuncControls::KeywordFamily::Undefined: - throw std::domain_error("No valid saturation keyword family specified"); - } - } + unsigned satRegionIdx); void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, - const SgofTable& sgofTable) - { - // convert the saturations of the SGOF keyword from gas to oil saturations - std::vector SoSamples(sgofTable.numRows()); - for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) { - SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx); - } - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - realParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KROG"))); - realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KRG"))); - realParams.setPcnwSamples(SoSamples, sgofTable.getColumn("PCOG").vectorCopy()); - realParams.finalize(); - } + const SgofTable& sgofTable); void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, - const SlgofTable& slgofTable) - { - // convert the saturations of the SLGOF keyword from "liquid" to oil saturations - std::vector SoSamples(slgofTable.numRows()); - for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) { - SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco; - } - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - realParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KROG"))); - realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KRG"))); - realParams.setPcnwSamples(SoSamples, slgofTable.getColumn("PCOG").vectorCopy()); - realParams.finalize(); - } + const SlgofTable& slgofTable); void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, const Sof3Table& sof3Table, - const SgfnTable& sgfnTable) - { - // convert the saturations of the SGFN keyword from gas to oil saturations - std::vector SoSamples(sgfnTable.numRows()); - std::vector SoColumn = sof3Table.getColumn("SO").vectorCopy(); - for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) { - SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); - } - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROG"))); - realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); - realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy()); - realParams.finalize(); - } + const SgfnTable& sgfnTable); void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams, const Scalar Swco, const double tolcrit, const Sof2Table& sof2Table, - const SgfnTable& sgfnTable) - { - // convert the saturations of the SGFN keyword from gas to oil saturations - std::vector SoSamples(sgfnTable.numRows()); - std::vector SoColumn = sof2Table.getColumn("SO").vectorCopy(); - for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) { - SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); - } - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO"))); - realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); - realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy()); - realParams.finalize(); - } + const SgfnTable& sgfnTable); template void readOilWaterEffectiveParameters_(Container& dest, const EclipseState& eclState, - unsigned satRegionIdx) - { - if (!hasOil || !hasWater) - // we don't read anything if either the water or the oil phase is not active - return; - - dest[satRegionIdx] = std::make_shared(); - - const auto tolcrit = eclState.runspec().saturationFunctionControls() - .minimumRelpermMobilityThreshold(); - - const auto& tableManager = eclState.getTableManager(); - auto& effParams = *dest[satRegionIdx]; - - switch (eclState.runspec().saturationFunctionControls().family()) { - case SatFuncControls::KeywordFamily::Family_I: - { - if (tableManager.hasTables("SWOF")) { - const auto& swofTable = tableManager.getSwofTables().getTable(satRegionIdx); - const std::vector SwColumn = swofTable.getColumn("SW").vectorCopy(); - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KRW"))); - realParams.setKrnSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KROW"))); - realParams.setPcnwSamples(SwColumn, swofTable.getColumn("PCOW").vectorCopy()); - realParams.finalize(); - } - else if ( !tableManager.getSwofletTable().empty() ) { - const auto& letTab = tableManager.getSwofletTable()[satRegionIdx]; - const std::vector dum; // dummy arg to conform with existing interface - - effParams.setApproach(SatCurveMultiplexerApproach::LET); - auto& realParams = effParams.template getRealParams(); - - // S=(Sw-Swcr)/(1-Sowcr-Swcr), krw = Krt*S^L/[S^L+E*(1.0-S)^T] - const Scalar s_min_w = letTab.s1_critical; - const Scalar s_max_w = 1.0-letTab.s2_critical; - const std::vector& letCoeffsWat = {s_min_w, s_max_w, - static_cast(letTab.l1_relperm), - static_cast(letTab.e1_relperm), - static_cast(letTab.t1_relperm), - static_cast(letTab.krt1_relperm)}; - realParams.setKrwSamples(letCoeffsWat, dum); - - // S=(So-Sowcr)/(1-Sowcr-Swcr), krow = Krt*S^L/[S^L+E*(1.0-S)^T] - const Scalar s_min_nw = letTab.s2_critical; - const Scalar s_max_nw = 1.0-letTab.s1_critical; - const std::vector& letCoeffsOil = {s_min_nw, s_max_nw, - static_cast(letTab.l2_relperm), - static_cast(letTab.e2_relperm), - static_cast(letTab.t2_relperm), - static_cast(letTab.krt2_relperm)}; - realParams.setKrnSamples(letCoeffsOil, dum); - - // S=(Sw-Swco)/(1-Swco-Sorw), Pc = Pct + (Pcir-Pct)*(1-S)^L/[(1-S)^L+E*S^T] - const std::vector& letCoeffsPc = {static_cast(letTab.s1_residual), - static_cast(letTab.s2_residual), - static_cast(letTab.l_pc), - static_cast(letTab.e_pc), - static_cast(letTab.t_pc), - static_cast(letTab.pcir_pc), - static_cast(letTab.pct_pc)}; - realParams.setPcnwSamples(letCoeffsPc, dum); - - realParams.finalize(); - } - break; - } - - case SatFuncControls::KeywordFamily::Family_II: - { - const auto& swfnTable = tableManager.getSwfnTables().getTable(satRegionIdx); - const std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); - realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy()); - - if (!hasGas) { - const auto& sof2Table = tableManager.getSof2Tables().getTable(satRegionIdx); - // convert the saturations of the SOF2 keyword from oil to water saturations - std::vector SwSamples(sof2Table.numRows()); - for (size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx) - SwSamples[sampleIdx] = 1 - sof2Table.get("SO", sampleIdx); - - realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO"))); - } else { - const auto& sof3Table = tableManager.getSof3Tables().getTable(satRegionIdx); - // convert the saturations of the SOF3 keyword from oil to water saturations - std::vector SwSamples(sof3Table.numRows()); - for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx) - SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx); - - realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROW"))); - } - realParams.finalize(); - break; - } - - case SatFuncControls::KeywordFamily::Undefined: - throw std::domain_error("No valid saturation keyword family specified"); - } - } + unsigned satRegionIdx); template void readGasWaterEffectiveParameters_(Container& dest, const EclipseState& eclState, - unsigned satRegionIdx) - { - if (!hasGas || !hasWater || hasOil) - // we don't read anything if either the gas or the water phase is not active or if oil is present - return; - - dest[satRegionIdx] = std::make_shared(); - - auto& effParams = *dest[satRegionIdx]; - - const auto tolcrit = eclState.runspec().saturationFunctionControls() - .minimumRelpermMobilityThreshold(); - - const auto& tableManager = eclState.getTableManager(); - - switch (eclState.runspec().saturationFunctionControls().family()) { - case SatFuncControls::KeywordFamily::Family_I: - { - throw std::domain_error("Saturation keyword family I is not applicable for a gas-water system"); - } - - case SatFuncControls::KeywordFamily::Family_II: - { - //Todo: allow also for Sgwfn table input as alternative to Sgfn and Swfn table input - const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable( satRegionIdx ); - const SwfnTable& swfnTable = tableManager.getSwfnTables().getTable( satRegionIdx ); - - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); - - std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); - - realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); - std::vector SwSamples(sgfnTable.numRows()); - for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) - SwSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx); - realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); - //Capillary pressure is read from SWFN. - //For gas-water system the capillary pressure column values are set to 0 in SGFN - realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy()); - realParams.finalize(); - - break; - } - - case SatFuncControls::KeywordFamily::Undefined: - throw std::domain_error("No valid saturation keyword family specified"); - } - } + unsigned satRegionIdx); template void readGasOilUnscaledPoints_(Container& dest, std::shared_ptr config, const EclipseState& /* eclState */, - unsigned satRegionIdx) - { - if (!hasGas || !hasOil) - // we don't read anything if either the gas or the oil phase is not active - return; - - dest[satRegionIdx] = std::make_shared >(); - dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], - *config, - EclTwoPhaseSystemType::GasOil); - } + unsigned satRegionIdx); template void readOilWaterUnscaledPoints_(Container& dest, std::shared_ptr config, const EclipseState& /* eclState */, - unsigned satRegionIdx) - { - if (!hasOil || !hasWater) - // we don't read anything if either the water or the oil phase is not active - return; - - dest[satRegionIdx] = std::make_shared >(); - dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], - *config, - EclTwoPhaseSystemType::OilWater); - } + unsigned satRegionIdx); template void readGasWaterUnscaledPoints_(Container& dest, std::shared_ptr config, const EclipseState& /* eclState */, - unsigned satRegionIdx) - { - if (hasOil) - // we don't read anything if oil phase is active - return; - - dest[satRegionIdx] = std::make_shared >(); - dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], - *config, - EclTwoPhaseSystemType::GasWater); - } + unsigned satRegionIdx); std::tuple, EclEpsScalingPoints> @@ -1161,18 +284,7 @@ class EclMaterialLawManager const EclipseState& eclState, const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, - EclTwoPhaseSystemType type) - { - unsigned satRegionIdx = epsGridProperties.satRegion( elemIdx ); - - EclEpsScalingPointsInfo destInfo(unscaledEpsInfo_[satRegionIdx]); - destInfo.extractScaled(eclState, epsGridProperties, elemIdx); - - EclEpsScalingPoints destPoint; - destPoint.init(destInfo, config, type); - - return {destInfo, destPoint}; - } + EclTwoPhaseSystemType type); void initThreePhaseParams_(const EclipseState& /* eclState */, MaterialLawParams& materialParams, @@ -1180,74 +292,7 @@ class EclMaterialLawManager const EclEpsScalingPointsInfo& epsInfo, std::shared_ptr oilWaterParams, std::shared_ptr gasOilParams, - std::shared_ptr gasWaterParams) - { - materialParams.setApproach(threePhaseApproach_); - - switch (materialParams.approach()) { - case EclMultiplexerApproach::Stone1: { - auto& realParams = materialParams.template getRealParams(); - realParams.setGasOilParams(gasOilParams); - realParams.setOilWaterParams(oilWaterParams); - realParams.setSwl(epsInfo.Swl); - - if (!stoneEtas.empty()) { - realParams.setEta(stoneEtas[satRegionIdx]); - } - else - realParams.setEta(1.0); - realParams.finalize(); - break; - } - - case EclMultiplexerApproach::Stone2: { - auto& realParams = materialParams.template getRealParams(); - realParams.setGasOilParams(gasOilParams); - realParams.setOilWaterParams(oilWaterParams); - realParams.setSwl(epsInfo.Swl); - realParams.finalize(); - break; - } - - case EclMultiplexerApproach::Default: { - auto& realParams = materialParams.template getRealParams(); - realParams.setGasOilParams(gasOilParams); - realParams.setOilWaterParams(oilWaterParams); - realParams.setSwl(epsInfo.Swl); - realParams.finalize(); - break; - } - - case EclMultiplexerApproach::TwoPhase: { - auto& realParams = materialParams.template getRealParams(); - realParams.setGasOilParams(gasOilParams); - realParams.setOilWaterParams(oilWaterParams); - realParams.setGasWaterParams(gasWaterParams); - realParams.setApproach(twoPhaseApproach_); - realParams.finalize(); - break; - } - - case EclMultiplexerApproach::OnePhase: { - // Nothing to do, no parameters. - break; - } - } - } - - // Relative permeability values not strictly greater than 'tolcrit' treated as zero. - std::vector normalizeKrValues_(const double tolcrit, - const TableColumn& krValues) const - { - auto kr = krValues.vectorCopy(); - std::transform(kr.begin(), kr.end(), kr.begin(), - [tolcrit](const double kri) - { - return (kri > tolcrit) ? kri : 0.0; - }); - - return kr; - } + std::shared_ptr gasWaterParams); bool enableEndPointScaling_; std::shared_ptr hysteresisConfig_; diff --git a/src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp b/src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp new file mode 100644 index 00000000000..fbb8402ee5f --- /dev/null +++ b/src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp @@ -0,0 +1,1162 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +namespace { + +// Relative permeability values not strictly greater than 'tolcrit' treated as zero. +std::vector normalizeKrValues_(const double tolcrit, + const Opm::TableColumn& krValues) +{ + auto kr = krValues.vectorCopy(); + std::transform(kr.begin(), kr.end(), kr.begin(), + [tolcrit](const double kri) + { + return (kri > tolcrit) ? kri : 0.0; + }); + + return kr; +} + +} + +namespace Opm { + +template +EclMaterialLawManager::EclMaterialLawManager() = default; + +template +EclMaterialLawManager::~EclMaterialLawManager() = default; + +template +void EclMaterialLawManager:: +initFromState(const EclipseState& eclState) +{ + // get the number of saturation regions and the number of cells in the deck + const auto& runspec = eclState.runspec(); + const size_t numSatRegions = runspec.tabdims().getNumSatTables(); + + const auto& ph = runspec.phases(); + this->hasGas = ph.active(Phase::GAS); + this->hasOil = ph.active(Phase::OIL); + this->hasWater = ph.active(Phase::WATER); + + readGlobalEpsOptions_(eclState); + readGlobalHysteresisOptions_(eclState); + readGlobalThreePhaseOptions_(runspec); + + // Read the end point scaling configuration (once per run). + gasOilConfig = std::make_shared(); + oilWaterConfig = std::make_shared(); + gasWaterConfig = std::make_shared(); + gasOilConfig->initFromState(eclState, EclTwoPhaseSystemType::GasOil); + oilWaterConfig->initFromState(eclState, EclTwoPhaseSystemType::OilWater); + gasWaterConfig->initFromState(eclState, EclTwoPhaseSystemType::GasWater); + + + const auto& tables = eclState.getTableManager(); + + { + const auto& stone1exTables = tables.getStone1exTable(); + + if (! stone1exTables.empty()) { + stoneEtas.clear(); + stoneEtas.reserve(numSatRegions); + + for (const auto& table : stone1exTables) { + stoneEtas.push_back(table.eta); + } + } + } + + this->unscaledEpsInfo_.resize(numSatRegions); + + if (this->hasGas + this->hasOil + this->hasWater == 1) { + // Single-phase simulation. Special case. Nothing to do here. + return; + } + + // Multiphase simulation. Common case. + const auto tolcrit = runspec.saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto rtep = satfunc::getRawTableEndpoints(tables, ph, tolcrit); + const auto rfunc = satfunc::getRawFunctionValues(tables, ph, rtep); + + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { + this->unscaledEpsInfo_[satRegionIdx] + .extractUnscaled(rtep, rfunc, satRegionIdx); + } +} + +template +void EclMaterialLawManager:: +initParamsForElements(const EclipseState& eclState, size_t numCompressedElems) +{ + // get the number of saturation regions + const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables(); + + // setup the saturation region specific parameters + gasOilUnscaledPointsVector_.resize(numSatRegions); + oilWaterUnscaledPointsVector_.resize(numSatRegions); + gasWaterUnscaledPointsVector_.resize(numSatRegions); + + gasOilEffectiveParamVector_.resize(numSatRegions); + oilWaterEffectiveParamVector_.resize(numSatRegions); + gasWaterEffectiveParamVector_.resize(numSatRegions); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { + // unscaled points for end-point scaling + readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, eclState, satRegionIdx); + readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, eclState, satRegionIdx); + readGasWaterUnscaledPoints_(gasWaterUnscaledPointsVector_, gasWaterConfig, eclState, satRegionIdx); + + // the parameters for the effective two-phase matererial laws + readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, eclState, satRegionIdx); + readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, eclState, satRegionIdx); + readGasWaterEffectiveParameters_(gasWaterEffectiveParamVector_, eclState, satRegionIdx); + } + + // copy the SATNUM grid property. in some cases this is not necessary, but it + // should not require much memory anyway... + satnumRegionArray_.resize(numCompressedElems); + if (eclState.fieldProps().has_int("SATNUM")) { + const auto& satnumRawData = eclState.fieldProps().get_int("SATNUM"); + for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { + satnumRegionArray_[elemIdx] = satnumRawData[elemIdx] - 1; + } + } + else { + std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0); + } + auto copy_krnum = [&eclState, numCompressedElems](std::vector& dest, const std::string keyword) { + if (eclState.fieldProps().has_int(keyword)) { + dest.resize(numCompressedElems); + const auto& satnumRawData = eclState.fieldProps().get_int(keyword); + for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { + dest[elemIdx] = satnumRawData[elemIdx] - 1; + } + } + }; + copy_krnum(krnumXArray_, "KRNUMX"); + copy_krnum(krnumYArray_, "KRNUMY"); + copy_krnum(krnumZArray_, "KRNUMZ"); + + // create the information for the imbibition region (IMBNUM). By default this is + // the same as the saturation region (SATNUM) + imbnumRegionArray_ = satnumRegionArray_; + if (eclState.fieldProps().has_int("IMBNUM")) { + const auto& imbnumRawData = eclState.fieldProps().get_int("IMBNUM"); + for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { + imbnumRegionArray_[elemIdx] = imbnumRawData[elemIdx] - 1; + } + } + + assert(numCompressedElems == satnumRegionArray_.size()); + assert(!enableHysteresis() || numCompressedElems == imbnumRegionArray_.size()); + + // read the scaled end point scaling parameters which are specific for each + // element + oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems); + + std::unique_ptr epsImbGridProperties; + + if (enableHysteresis()) { + epsImbGridProperties = std::make_unique(eclState, true); + } + + EclEpsGridProperties epsGridProperties(eclState, false); + materialLawParams_.resize(numCompressedElems); + + for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) { + unsigned satRegionIdx = static_cast(satnumRegionArray_[elemIdx]); + auto gasOilParams = std::make_shared(); + auto oilWaterParams = std::make_shared(); + auto gasWaterParams = std::make_shared(); + gasOilParams->setConfig(hysteresisConfig_); + oilWaterParams->setConfig(hysteresisConfig_); + gasWaterParams->setConfig(hysteresisConfig_); + + auto [gasOilScaledInfo, gasOilScaledPoint] = + readScaledPoints_(*gasOilConfig, + eclState, + epsGridProperties, + elemIdx, + EclTwoPhaseSystemType::GasOil); + + auto [owinfo, oilWaterScaledEpsPointDrainage] = + readScaledPoints_(*oilWaterConfig, + eclState, + epsGridProperties, + elemIdx, + EclTwoPhaseSystemType::OilWater); + oilWaterScaledEpsInfoDrainage_[elemIdx] = owinfo; + + auto [gasWaterScaledInfo, gasWaterScaledPoint] = + readScaledPoints_(*gasWaterConfig, + eclState, + epsGridProperties, + elemIdx, + EclTwoPhaseSystemType::GasWater); + + if (hasGas && hasOil) { + GasOilEpsTwoPhaseParams gasOilDrainParams; + gasOilDrainParams.setConfig(gasOilConfig); + gasOilDrainParams.setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); + gasOilDrainParams.setScaledPoints(gasOilScaledPoint); + gasOilDrainParams.setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); + gasOilDrainParams.finalize(); + + gasOilParams->setDrainageParams(gasOilDrainParams, + gasOilScaledInfo, + EclTwoPhaseSystemType::GasOil); + } + + if (hasOil && hasWater) { + OilWaterEpsTwoPhaseParams oilWaterDrainParams; + oilWaterDrainParams.setConfig(oilWaterConfig); + oilWaterDrainParams.setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); + oilWaterDrainParams.setScaledPoints(oilWaterScaledEpsPointDrainage); + oilWaterDrainParams.setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); + oilWaterDrainParams.finalize(); + + oilWaterParams->setDrainageParams(oilWaterDrainParams, + owinfo, + EclTwoPhaseSystemType::OilWater); + } + + if (hasGas && hasWater && !hasOil) { + GasWaterEpsTwoPhaseParams gasWaterDrainParams; + gasWaterDrainParams.setConfig(gasWaterConfig); + gasWaterDrainParams.setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]); + gasWaterDrainParams.setScaledPoints(gasWaterScaledPoint); + gasWaterDrainParams.setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]); + gasWaterDrainParams.finalize(); + + gasWaterParams->setDrainageParams(gasWaterDrainParams, + gasWaterScaledInfo, + EclTwoPhaseSystemType::GasWater); + } + + if (enableHysteresis()) { + auto [gasOilScaledImbInfo, gasOilScaledImbPoint] = + readScaledPoints_(*gasOilConfig, + eclState, + *epsImbGridProperties, + elemIdx, + EclTwoPhaseSystemType::GasOil); + + auto [oilWaterScaledImbInfo, oilWaterScaledImbPoint] = + readScaledPoints_(*oilWaterConfig, + eclState, + *epsImbGridProperties, + elemIdx, + EclTwoPhaseSystemType::OilWater); + + auto [gasWaterScaledImbInfo, gasWaterScaledImbPoint] = + readScaledPoints_(*gasWaterConfig, + eclState, + *epsImbGridProperties, + elemIdx, + EclTwoPhaseSystemType::GasWater); + + unsigned imbRegionIdx = imbnumRegionArray_[elemIdx]; + if (hasGas && hasOil) { + GasOilEpsTwoPhaseParams gasOilImbParamsHyst; + gasOilImbParamsHyst.setConfig(gasOilConfig); + gasOilImbParamsHyst.setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]); + gasOilImbParamsHyst.setScaledPoints(gasOilScaledImbPoint); + gasOilImbParamsHyst.setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]); + gasOilImbParamsHyst.finalize(); + + gasOilParams->setImbibitionParams(gasOilImbParamsHyst, + gasOilScaledImbInfo, + EclTwoPhaseSystemType::GasOil); + } + + if (hasOil && hasWater) { + OilWaterEpsTwoPhaseParams oilWaterImbParamsHyst; + oilWaterImbParamsHyst.setConfig(oilWaterConfig); + oilWaterImbParamsHyst.setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]); + oilWaterImbParamsHyst.setScaledPoints(oilWaterScaledImbPoint); + oilWaterImbParamsHyst.setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]); + oilWaterImbParamsHyst.finalize(); + + oilWaterParams->setImbibitionParams(oilWaterImbParamsHyst, + oilWaterScaledImbInfo, + EclTwoPhaseSystemType::OilWater); + } + + if (hasGas && hasWater && !hasOil) { + GasWaterEpsTwoPhaseParams gasWaterImbParamsHyst; + gasWaterImbParamsHyst.setConfig(gasWaterConfig); + gasWaterImbParamsHyst.setUnscaledPoints(gasWaterUnscaledPointsVector_[imbRegionIdx]); + gasWaterImbParamsHyst.setScaledPoints(gasWaterScaledImbPoint); + gasWaterImbParamsHyst.setEffectiveLawParams(gasWaterEffectiveParamVector_[imbRegionIdx]); + gasWaterImbParamsHyst.finalize(); + + gasWaterParams->setImbibitionParams(gasWaterImbParamsHyst, + gasWaterScaledImbInfo, + EclTwoPhaseSystemType::GasWater); + } + } + + if (hasGas && hasOil) + gasOilParams->finalize(); + + if (hasOil && hasWater) + oilWaterParams->finalize(); + + if (hasGas && hasWater && !hasOil) + gasWaterParams->finalize(); + + initThreePhaseParams_(eclState, + materialLawParams_[elemIdx], + satRegionIdx, + oilWaterScaledEpsInfoDrainage_[elemIdx], + oilWaterParams, + gasOilParams, + gasWaterParams); + + materialLawParams_[elemIdx].finalize(); + } +} + +template +typename TraitsT::Scalar EclMaterialLawManager:: +applySwatinit(unsigned elemIdx, + Scalar pcow, + Scalar Sw) +{ + auto& elemScaledEpsInfo = oilWaterScaledEpsInfoDrainage_[elemIdx]; + + // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74 + + if (pcow < 0.0) + Sw = elemScaledEpsInfo.Swu; + else { + + if (Sw <= elemScaledEpsInfo.Swl) + Sw = elemScaledEpsInfo.Swl; + + // specify a fluid state which only stores the saturations + using FluidState = SimpleModularFluidState don't care */ + /*storePressure=*/false, + /*storeTemperature=*/false, + /*storeComposition=*/false, + /*storeFugacity=*/false, + /*storeSaturation=*/true, + /*storeDensity=*/false, + /*storeViscosity=*/false, + /*storeEnthalpy=*/false>; + FluidState fs; + fs.setSaturation(waterPhaseIdx, Sw); + fs.setSaturation(gasPhaseIdx, 0); + fs.setSaturation(oilPhaseIdx, 0); + std::array pc = { 0 }; + MaterialLaw::capillaryPressures(pc, materialLawParams(elemIdx), fs); + + Scalar pcowAtSw = pc[oilPhaseIdx] - pc[waterPhaseIdx]; + constexpr const Scalar pcowAtSwThreshold = 1.0; //Pascal + // avoid divison by very small number + if (std::abs(pcowAtSw) > pcowAtSwThreshold) { + elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw; + auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx); + elemEclEpsScalingPoints.init(elemScaledEpsInfo, + *oilWaterEclEpsConfig_, + EclTwoPhaseSystemType::OilWater); + } + } + + return Sw; +} + +template +const typename EclMaterialLawManager::MaterialLawParams& +EclMaterialLawManager:: +connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const +{ + MaterialLawParams& mlp = const_cast(materialLawParams_[elemIdx]); + + if (enableHysteresis()) + OpmLog::warning("Warning: Using non-default satnum regions for connection is not tested in combination with hysteresis"); + // Currently we don't support COMPIMP. I.e. use the same table lookup for the hysteresis curves. + // unsigned impRegionIdx = satRegionIdx; + + // change the sat table it points to. + switch (mlp.approach()) { + case EclMultiplexerApproach::Stone1: { + auto& realParams = mlp.template getRealParams(); + + realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); + realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); +// if (enableHysteresis()) { +// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); +// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); +// } + } + break; + + case EclMultiplexerApproach::Stone2: { + auto& realParams = mlp.template getRealParams(); + realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); + realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); +// if (enableHysteresis()) { +// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); +// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); +// } + } + break; + + case EclMultiplexerApproach::Default: { + auto& realParams = mlp.template getRealParams(); + realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); + realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); +// if (enableHysteresis()) { +// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); +// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); +// } + } + break; + + case EclMultiplexerApproach::TwoPhase: { + auto& realParams = mlp.template getRealParams(); + realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); + realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); +// if (enableHysteresis()) { +// realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]); +// realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]); +// realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]); +// } + } + break; + + default: + throw std::logic_error("Enum value for material approach unknown!"); + } + + return mlp; +} + +template +int EclMaterialLawManager:: +getKrnumSatIdx(unsigned elemIdx, FaceDir::DirEnum facedir) const +{ + using Dir = FaceDir::DirEnum; + const std::vector* array = nullptr; + switch(facedir) { + case Dir::XPlus: + array = &krnumXArray_; + break; + case Dir::YPlus: + array = &krnumYArray_; + break; + case Dir::ZPlus: + array = &krnumZArray_; + break; + default: + throw std::runtime_error("Unknown face direction"); + } + if (array->size() > 0) { + return (*array)[elemIdx]; + } + else { + return satnumRegionArray_[elemIdx]; + } +} + +template +void EclMaterialLawManager:: +oilWaterHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + unsigned elemIdx) const +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled."); + + const auto& params = materialLawParams(elemIdx); + MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params); +} + +template +void EclMaterialLawManager:: +setOilWaterHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + unsigned elemIdx) +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled."); + + auto& params = materialLawParams(elemIdx); + MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params); +} + +template +void EclMaterialLawManager:: +gasOilHysteresisParams(Scalar& pcSwMdc, + Scalar& krnSwMdc, + unsigned elemIdx) const +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled."); + + const auto& params = materialLawParams(elemIdx); + MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params); +} + +template +void EclMaterialLawManager:: +setGasOilHysteresisParams(const Scalar& pcSwMdc, + const Scalar& krnSwMdc, + unsigned elemIdx) +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled."); + + auto& params = materialLawParams(elemIdx); + MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params); +} + +template +EclEpsScalingPoints& +EclMaterialLawManager:: +oilWaterScaledEpsPointsDrainage(unsigned elemIdx) +{ + auto& materialParams = materialLawParams_[elemIdx]; + switch (materialParams.approach()) { + case EclMultiplexerApproach::Stone1: { + auto& realParams = materialParams.template getRealParams(); + return realParams.oilWaterParams().drainageParams().scaledPoints(); + } + + case EclMultiplexerApproach::Stone2: { + auto& realParams = materialParams.template getRealParams(); + return realParams.oilWaterParams().drainageParams().scaledPoints(); + } + + case EclMultiplexerApproach::Default: { + auto& realParams = materialParams.template getRealParams(); + return realParams.oilWaterParams().drainageParams().scaledPoints(); + } + + case EclMultiplexerApproach::TwoPhase: { + auto& realParams = materialParams.template getRealParams(); + return realParams.oilWaterParams().drainageParams().scaledPoints(); + } + default: + throw std::logic_error("Enum value for material approach unknown!"); + } +} + +template +void EclMaterialLawManager:: +readGlobalEpsOptions_(const EclipseState& eclState) +{ + oilWaterEclEpsConfig_ = std::make_shared(); + oilWaterEclEpsConfig_->initFromState(eclState, EclTwoPhaseSystemType::OilWater); + + enableEndPointScaling_ = eclState.getTableManager().hasTables("ENKRVD"); +} + +template +void EclMaterialLawManager:: +readGlobalHysteresisOptions_(const EclipseState& state) +{ + hysteresisConfig_ = std::make_shared(); + hysteresisConfig_->initFromState(state.runspec()); +} + +template +void EclMaterialLawManager:: +readGlobalThreePhaseOptions_(const Runspec& runspec) +{ + bool gasEnabled = runspec.phases().active(Phase::GAS); + bool oilEnabled = runspec.phases().active(Phase::OIL); + bool waterEnabled = runspec.phases().active(Phase::WATER); + + int numEnabled = + (gasEnabled?1:0) + + (oilEnabled?1:0) + + (waterEnabled?1:0); + + if (numEnabled == 0) { + throw std::runtime_error("At least one fluid phase must be enabled. (Is: "+std::to_string(numEnabled)+")"); + } else if (numEnabled == 1) { + threePhaseApproach_ = EclMultiplexerApproach::OnePhase; + } else if ( numEnabled == 2) { + threePhaseApproach_ = EclMultiplexerApproach::TwoPhase; + if (!gasEnabled) + twoPhaseApproach_ = EclTwoPhaseApproach::OilWater; + else if (!oilEnabled) + twoPhaseApproach_ = EclTwoPhaseApproach::GasWater; + else if (!waterEnabled) + twoPhaseApproach_ = EclTwoPhaseApproach::GasOil; + } + else { + assert(numEnabled == 3); + + threePhaseApproach_ = EclMultiplexerApproach::Default; + const auto& satctrls = runspec.saturationFunctionControls(); + if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone2) + threePhaseApproach_ = EclMultiplexerApproach::Stone2; + else if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone1) + threePhaseApproach_ = EclMultiplexerApproach::Stone1; + } +} + +template +template +void EclMaterialLawManager:: +readGasOilEffectiveParameters_(Container& dest, + const EclipseState& eclState, + unsigned satRegionIdx) +{ + if (!hasGas || !hasOil) + // we don't read anything if either the gas or the oil phase is not active + return; + + dest[satRegionIdx] = std::make_shared(); + + auto& effParams = *dest[satRegionIdx]; + + // the situation for the gas phase is complicated that all saturations are + // shifted by the connate water saturation. + const Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl; + const auto tolcrit = eclState.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto& tableManager = eclState.getTableManager(); + + switch (eclState.runspec().saturationFunctionControls().family()) { + case SatFuncControls::KeywordFamily::Family_I: + { + const TableContainer& sgofTables = tableManager.getSgofTables(); + const TableContainer& slgofTables = tableManager.getSlgofTables(); + if (!sgofTables.empty()) + readGasOilEffectiveParametersSgof_(effParams, Swco, tolcrit, + sgofTables.getTable(satRegionIdx)); + else if (!slgofTables.empty()) + readGasOilEffectiveParametersSlgof_(effParams, Swco, tolcrit, + slgofTables.getTable(satRegionIdx)); + else if ( !tableManager.getSgofletTable().empty() ) { + const auto& letSgofTab = tableManager.getSgofletTable()[satRegionIdx]; + const std::vector dum; // dummy arg to comform with existing interface + + effParams.setApproach(SatCurveMultiplexerApproach::LET); + auto& realParams = effParams.template getRealParams(); + + // S=(So-Sogcr)/(1-Sogcr-Sgcr-Swco), krog = Krt*S^L/[S^L+E*(1.0-S)^T] + const Scalar s_min_w = letSgofTab.s2_critical; + const Scalar s_max_w = 1.0-letSgofTab.s1_critical-Swco; + const std::vector& letCoeffsOil = {s_min_w, s_max_w, + static_cast(letSgofTab.l2_relperm), + static_cast(letSgofTab.e2_relperm), + static_cast(letSgofTab.t2_relperm), + static_cast(letSgofTab.krt2_relperm)}; + realParams.setKrwSamples(letCoeffsOil, dum); + + // S=(1-So-Sgcr-Swco)/(1-Sogcr-Sgcr-Swco), krg = Krt*S^L/[S^L+E*(1.0-S)^T] + const Scalar s_min_nw = letSgofTab.s1_critical+Swco; + const Scalar s_max_nw = 1.0-letSgofTab.s2_critical; + const std::vector& letCoeffsGas = {s_min_nw, s_max_nw, + static_cast(letSgofTab.l1_relperm), + static_cast(letSgofTab.e1_relperm), + static_cast(letSgofTab.t1_relperm), + static_cast(letSgofTab.krt1_relperm)}; + realParams.setKrnSamples(letCoeffsGas, dum); + + // S=(So-Sorg)/(1-Sorg-Sgl-Swco), Pc = Pct + (pcir_pc-Pct)*(1-S)^L/[(1-S)^L+E*S^T] + const std::vector& letCoeffsPc = {static_cast(letSgofTab.s2_residual), + static_cast(letSgofTab.s1_residual+Swco), + static_cast(letSgofTab.l_pc), + static_cast(letSgofTab.e_pc), + static_cast(letSgofTab.t_pc), + static_cast(letSgofTab.pcir_pc), + static_cast(letSgofTab.pct_pc)}; + realParams.setPcnwSamples(letCoeffsPc, dum); + + realParams.finalize(); + } + break; + } + + case SatFuncControls::KeywordFamily::Family_II: + { + const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable( satRegionIdx ); + if (!hasWater) { + // oil and gas case + const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable( satRegionIdx ); + readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable); + } + else { + const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable( satRegionIdx ); + readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable); + } + break; + } + + case SatFuncControls::KeywordFamily::Undefined: + throw std::domain_error("No valid saturation keyword family specified"); + } +} + +template +void EclMaterialLawManager:: +readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const SgofTable& sgofTable) +{ + // convert the saturations of the SGOF keyword from gas to oil saturations + std::vector SoSamples(sgofTable.numRows()); + for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) { + SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx); + } + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + realParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KROG"))); + realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KRG"))); + realParams.setPcnwSamples(SoSamples, sgofTable.getColumn("PCOG").vectorCopy()); + realParams.finalize(); +} + +template +void EclMaterialLawManager:: +readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const SlgofTable& slgofTable) +{ + // convert the saturations of the SLGOF keyword from "liquid" to oil saturations + std::vector SoSamples(slgofTable.numRows()); + for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) { + SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco; + } + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + realParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KROG"))); + realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KRG"))); + realParams.setPcnwSamples(SoSamples, slgofTable.getColumn("PCOG").vectorCopy()); + realParams.finalize(); +} + +template +void EclMaterialLawManager:: +readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const Sof3Table& sof3Table, + const SgfnTable& sgfnTable) +{ + // convert the saturations of the SGFN keyword from gas to oil saturations + std::vector SoSamples(sgfnTable.numRows()); + std::vector SoColumn = sof3Table.getColumn("SO").vectorCopy(); + for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) { + SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); + } + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROG"))); + realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); + realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy()); + realParams.finalize(); +} + +template +void EclMaterialLawManager:: +readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const Sof2Table& sof2Table, + const SgfnTable& sgfnTable) +{ + // convert the saturations of the SGFN keyword from gas to oil saturations + std::vector SoSamples(sgfnTable.numRows()); + std::vector SoColumn = sof2Table.getColumn("SO").vectorCopy(); + for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) { + SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); + } + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO"))); + realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); + realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy()); + realParams.finalize(); +} + +template +template +void EclMaterialLawManager:: +readOilWaterEffectiveParameters_(Container& dest, + const EclipseState& eclState, + unsigned satRegionIdx) +{ + if (!hasOil || !hasWater) + // we don't read anything if either the water or the oil phase is not active + return; + + dest[satRegionIdx] = std::make_shared(); + + const auto tolcrit = eclState.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto& tableManager = eclState.getTableManager(); + auto& effParams = *dest[satRegionIdx]; + + switch (eclState.runspec().saturationFunctionControls().family()) { + case SatFuncControls::KeywordFamily::Family_I: + { + if (tableManager.hasTables("SWOF")) { + const auto& swofTable = tableManager.getSwofTables().getTable(satRegionIdx); + const std::vector SwColumn = swofTable.getColumn("SW").vectorCopy(); + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KRW"))); + realParams.setKrnSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KROW"))); + realParams.setPcnwSamples(SwColumn, swofTable.getColumn("PCOW").vectorCopy()); + realParams.finalize(); + } + else if ( !tableManager.getSwofletTable().empty() ) { + const auto& letTab = tableManager.getSwofletTable()[satRegionIdx]; + const std::vector dum; // dummy arg to conform with existing interface + + effParams.setApproach(SatCurveMultiplexerApproach::LET); + auto& realParams = effParams.template getRealParams(); + + // S=(Sw-Swcr)/(1-Sowcr-Swcr), krw = Krt*S^L/[S^L+E*(1.0-S)^T] + const Scalar s_min_w = letTab.s1_critical; + const Scalar s_max_w = 1.0-letTab.s2_critical; + const std::vector& letCoeffsWat = {s_min_w, s_max_w, + static_cast(letTab.l1_relperm), + static_cast(letTab.e1_relperm), + static_cast(letTab.t1_relperm), + static_cast(letTab.krt1_relperm)}; + realParams.setKrwSamples(letCoeffsWat, dum); + + // S=(So-Sowcr)/(1-Sowcr-Swcr), krow = Krt*S^L/[S^L+E*(1.0-S)^T] + const Scalar s_min_nw = letTab.s2_critical; + const Scalar s_max_nw = 1.0-letTab.s1_critical; + const std::vector& letCoeffsOil = {s_min_nw, s_max_nw, + static_cast(letTab.l2_relperm), + static_cast(letTab.e2_relperm), + static_cast(letTab.t2_relperm), + static_cast(letTab.krt2_relperm)}; + realParams.setKrnSamples(letCoeffsOil, dum); + + // S=(Sw-Swco)/(1-Swco-Sorw), Pc = Pct + (Pcir-Pct)*(1-S)^L/[(1-S)^L+E*S^T] + const std::vector& letCoeffsPc = {static_cast(letTab.s1_residual), + static_cast(letTab.s2_residual), + static_cast(letTab.l_pc), + static_cast(letTab.e_pc), + static_cast(letTab.t_pc), + static_cast(letTab.pcir_pc), + static_cast(letTab.pct_pc)}; + realParams.setPcnwSamples(letCoeffsPc, dum); + + realParams.finalize(); + } + break; + } + + case SatFuncControls::KeywordFamily::Family_II: + { + const auto& swfnTable = tableManager.getSwfnTables().getTable(satRegionIdx); + const std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); + realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy()); + + if (!hasGas) { + const auto& sof2Table = tableManager.getSof2Tables().getTable(satRegionIdx); + // convert the saturations of the SOF2 keyword from oil to water saturations + std::vector SwSamples(sof2Table.numRows()); + for (size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx) + SwSamples[sampleIdx] = 1 - sof2Table.get("SO", sampleIdx); + + realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO"))); + } else { + const auto& sof3Table = tableManager.getSof3Tables().getTable(satRegionIdx); + // convert the saturations of the SOF3 keyword from oil to water saturations + std::vector SwSamples(sof3Table.numRows()); + for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx) + SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx); + + realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROW"))); + } + realParams.finalize(); + break; + } + + case SatFuncControls::KeywordFamily::Undefined: + throw std::domain_error("No valid saturation keyword family specified"); + } +} + +template +template +void EclMaterialLawManager:: +readGasWaterEffectiveParameters_(Container& dest, + const EclipseState& eclState, + unsigned satRegionIdx) +{ + if (!hasGas || !hasWater || hasOil) + // we don't read anything if either the gas or the water phase is not active or if oil is present + return; + + dest[satRegionIdx] = std::make_shared(); + + auto& effParams = *dest[satRegionIdx]; + + const auto tolcrit = eclState.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto& tableManager = eclState.getTableManager(); + + switch (eclState.runspec().saturationFunctionControls().family()) { + case SatFuncControls::KeywordFamily::Family_I: + { + throw std::domain_error("Saturation keyword family I is not applicable for a gas-water system"); + } + + case SatFuncControls::KeywordFamily::Family_II: + { + //Todo: allow also for Sgwfn table input as alternative to Sgfn and Swfn table input + const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable( satRegionIdx ); + const SwfnTable& swfnTable = tableManager.getSwfnTables().getTable( satRegionIdx ); + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); + + realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); + std::vector SwSamples(sgfnTable.numRows()); + for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) + SwSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx); + realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); + //Capillary pressure is read from SWFN. + //For gas-water system the capillary pressure column values are set to 0 in SGFN + realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy()); + realParams.finalize(); + + break; + } + + case SatFuncControls::KeywordFamily::Undefined: + throw std::domain_error("No valid saturation keyword family specified"); + } +} + +template +template +void EclMaterialLawManager:: +readGasOilUnscaledPoints_(Container& dest, + std::shared_ptr config, + const EclipseState& /* eclState */, + unsigned satRegionIdx) +{ + if (!hasGas || !hasOil) + // we don't read anything if either the gas or the oil phase is not active + return; + + dest[satRegionIdx] = std::make_shared >(); + dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], + *config, + EclTwoPhaseSystemType::GasOil); +} + +template +template +void EclMaterialLawManager:: +readOilWaterUnscaledPoints_(Container& dest, + std::shared_ptr config, + const EclipseState& /* eclState */, + unsigned satRegionIdx) +{ + if (!hasOil || !hasWater) + // we don't read anything if either the water or the oil phase is not active + return; + + dest[satRegionIdx] = std::make_shared >(); + dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], + *config, + EclTwoPhaseSystemType::OilWater); +} + +template +template +void EclMaterialLawManager:: +readGasWaterUnscaledPoints_(Container& dest, + std::shared_ptr config, + const EclipseState& /* eclState */, + unsigned satRegionIdx) +{ + if (hasOil) + // we don't read anything if oil phase is active + return; + + dest[satRegionIdx] = std::make_shared >(); + dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], + *config, + EclTwoPhaseSystemType::GasWater); +} + +template +std::tuple, + EclEpsScalingPoints> +EclMaterialLawManager:: +readScaledPoints_(const EclEpsConfig& config, + const EclipseState& eclState, + const EclEpsGridProperties& epsGridProperties, + unsigned elemIdx, + EclTwoPhaseSystemType type) +{ + unsigned satRegionIdx = epsGridProperties.satRegion( elemIdx ); + + EclEpsScalingPointsInfo destInfo(unscaledEpsInfo_[satRegionIdx]); + destInfo.extractScaled(eclState, epsGridProperties, elemIdx); + + EclEpsScalingPoints destPoint; + destPoint.init(destInfo, config, type); + + return {destInfo, destPoint}; +} + +template +void EclMaterialLawManager:: +initThreePhaseParams_(const EclipseState& /* eclState */, + MaterialLawParams& materialParams, + unsigned satRegionIdx, + const EclEpsScalingPointsInfo& epsInfo, + std::shared_ptr oilWaterParams, + std::shared_ptr gasOilParams, + std::shared_ptr gasWaterParams) +{ + materialParams.setApproach(threePhaseApproach_); + + switch (materialParams.approach()) { + case EclMultiplexerApproach::Stone1: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + + if (!stoneEtas.empty()) { + realParams.setEta(stoneEtas[satRegionIdx]); + } + else + realParams.setEta(1.0); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::Stone2: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::Default: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::TwoPhase: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setGasWaterParams(gasWaterParams); + realParams.setApproach(twoPhaseApproach_); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::OnePhase: { + // Nothing to do, no parameters. + break; + } + } +} + +template class EclMaterialLawManager>; +template class EclMaterialLawManager>; + +} // namespace Opm