diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2d49dfcc291..f4d9b056fca 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -344,7 +344,11 @@ if(ENABLE_ECL_INPUT) opm/material/fluidmatrixinteractions/EclHysteresisConfig.cpp opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp + opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp + opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp + opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp + opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp opm/material/thermal/EclThermalLawManager.cpp ) @@ -997,6 +1001,7 @@ list( APPEND PUBLIC_HEADER_FILES opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp opm/material/fluidmatrixinteractions/TwoPhaseLETCurves.hpp opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp + opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp opm/material/fluidmatrixinteractions/DirectionalMaterialLawParams.hpp opm/material/fluidmatrixinteractions/DirectionalMaterialLawParams.hpp opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp diff --git a/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp b/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp index 061c288db78..50fef84183f 100644 --- a/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp +++ b/opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp @@ -90,6 +90,8 @@ class EclHysteresisTwoPhaseLaw : public EffectiveLawT::Traits //! are dependent on the phase composition static constexpr bool isCompositionDependent = false; + static constexpr bool isHysteresisDependent = true; + /*! * \brief The capillary pressure-saturation curves depending on absolute saturations. * diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp new file mode 100644 index 00000000000..8b9a956db15 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp @@ -0,0 +1,437 @@ +// -*- 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 + +namespace Opm { + +template +EclMaterialLawManagerSimple::EclMaterialLawManagerSimple() = default; + +template +EclMaterialLawManagerSimple::~EclMaterialLawManagerSimple() = default; + +template +void EclMaterialLawManagerSimple:: +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); + } + } + + const auto& ppcwmaxTables = tables.getPpcwmax(); + this->enablePpcwmax_ = !ppcwmaxTables.empty(); + + if (this->enablePpcwmax_) { + maxAllowPc_.clear(); + modifySwl_.clear(); + + maxAllowPc_.reserve(numSatRegions); + modifySwl_.reserve(numSatRegions); + + for (const auto& table : ppcwmaxTables) { + maxAllowPc_.push_back(table.max_cap_pres); + modifySwl_.push_back(table.option); + } + } + } + + 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); + } + + // WAG hysteresis parameters per SATNUM. + // if (eclState.runspec().hysterPar().activeWag()) { + // if (numSatRegions != eclState.getWagHysteresis().size()) + // throw std::runtime_error("Inconsistent Wag-hysteresis data"); + // wagHystersisConfig_.resize(numSatRegions); + // for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { + // wagHystersisConfig_[satRegionIdx] = std::make_shared(eclState.getWagHysteresis()[satRegionIdx]); + // } + // } +} + +template +void EclMaterialLawManagerSimple:: +initParamsForElements(const EclipseState& eclState, size_t numCompressedElems, + const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + InitParams initParams {*this, eclState, numCompressedElems}; + initParams.run(fieldPropIntOnLeafAssigner, lookupIdxOnLevelZeroAssigner); +} + +// TODO: Better (proper?) handling of mixed wettability systems - see ecl kw OPTIONS switch 74 +// Note: Without OPTIONS[74] the negative part of the Pcow curve is not scaled +template +std::pair EclMaterialLawManagerSimple:: +applySwatinit(unsigned elemIdx, + Scalar pcow, + Scalar Sw) +{ + // Default is no SWATINIT scaling of the negative part of the Pcow curve, so look up saturation using the input Pcow curve + if (pcow <= 0.0) { + return {Sw, /*newSwatInit*/ true}; + } + + auto& elemScaledEpsInfo = oilWaterScaledEpsInfoDrainage_[elemIdx]; + 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.0e-6; //Pascal + + // avoid division by very small number and avoid negative PCW at connate Sw + // (look up saturation on input Pcow curve in this case) + if (pcowAtSw < pcowAtSwThreshold) { + return {Sw, /*newSwatInit*/ true}; + } + + // Sufficiently positive value, continue with max. capillary pressure (PCW) scaling to honor SWATINIT value + Scalar newMaxPcow = elemScaledEpsInfo.maxPcow * (pcow/pcowAtSw); + + // Limit max. capillary pressure with PPCWMAX + bool newSwatInit = false; + int satRegionIdx = satnumRegionIdx(elemIdx); + if (enablePpcwmax() && (newMaxPcow > maxAllowPc_[satRegionIdx])) { + // Two options in PPCWMAX to modify connate Sw or not. In both cases, init. Sw needs to be + // re-calculated (done in opm-simulators) + newSwatInit = true; + if (modifySwl_[satRegionIdx] == false) { + // Max. cap. pressure set to PCWO in PPCWMAX + elemScaledEpsInfo.maxPcow = maxAllowPc_[satRegionIdx]; + } + else { + // Max. cap. pressure remains unscaled and connate Sw is set to SWATINIT value + elemScaledEpsInfo.Swl = Sw; + } + } + // Max. cap. pressure adjusted from SWATINIT data + else + elemScaledEpsInfo.maxPcow = newMaxPcow; + + // auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx); + // elemEclEpsScalingPoints.init(elemScaledEpsInfo, + // *oilWaterEclEpsConfig_, + // EclTwoPhaseSystemType::OilWater); + + return {Sw, newSwatInit}; +} + +template +void +EclMaterialLawManagerSimple::applyRestartSwatInit(const unsigned elemIdx, + const Scalar maxPcow) +{ + // Maximum capillary pressure adjusted from SWATINIT data. + + auto& elemScaledEpsInfo = + this->oilWaterScaledEpsInfoDrainage_[elemIdx]; + + elemScaledEpsInfo.maxPcow = maxPcow; + + // this->oilWaterScaledEpsPointsDrainage(elemIdx) + // .init(elemScaledEpsInfo, + // *this->oilWaterEclEpsConfig_, + // EclTwoPhaseSystemType::OilWater); +} + +template +const typename EclMaterialLawManagerSimple::MaterialLawParams& +EclMaterialLawManagerSimple:: +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"); + + auto& realParams = mlp; + if (realParams.approach() == EclTwoPhaseApproach::GasOil) { + realParams.setGasOilParams(gasOilEffectiveParamVector_[satRegionIdx]); + } + else if (realParams.approach() == EclTwoPhaseApproach::GasWater) { + realParams.setGasWaterParams(gasWaterEffectiveParamVector_[satRegionIdx]); + } + else if (realParams.approach() == EclTwoPhaseApproach::OilWater) { + realParams.setOilWaterParams(oilWaterEffectiveParamVector_[satRegionIdx]); + } + return mlp; +} + +template +int EclMaterialLawManagerSimple:: +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 EclMaterialLawManagerSimple:: +oilWaterHysteresisParams(Scalar& soMax, + Scalar& swMax, + Scalar& swMin, + unsigned elemIdx) const +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled."); + + const auto& params = materialLawParams(elemIdx); + MaterialLaw::oilWaterHysteresisParams(soMax, swMax, swMin, params); +} + +template +void EclMaterialLawManagerSimple:: +setOilWaterHysteresisParams(const Scalar& soMax, + const Scalar& swMax, + const Scalar& swMin, + unsigned elemIdx) +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled."); + + auto& params = materialLawParams(elemIdx); + MaterialLaw::setOilWaterHysteresisParams(soMax, swMax, swMin, params); +} + +template +void EclMaterialLawManagerSimple:: +gasOilHysteresisParams(Scalar& sgmax, + Scalar& shmax, + Scalar& somin, + unsigned elemIdx) const +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled."); + + const auto& params = materialLawParams(elemIdx); + MaterialLaw::gasOilHysteresisParams(sgmax, shmax, somin, params); +} + +template +void EclMaterialLawManagerSimple:: +setGasOilHysteresisParams(const Scalar& sgmax, + const Scalar& shmax, + const Scalar& somin, + unsigned elemIdx) +{ + if (!enableHysteresis()) + throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled."); + + auto& params = materialLawParams(elemIdx); + MaterialLaw::setGasOilHysteresisParams(sgmax, shmax, somin, params); +} + +// template +// EclEpsScalingPoints& +// EclMaterialLawManagerSimple:: +// oilWaterScaledEpsPointsDrainage(unsigned elemIdx) +// { +// auto& materialParams = materialLawParams_[elemIdx]; +// auto& realParams = materialParams; +// return realParams.oilWaterParams().drainageParams().scaledPoints(); +// } + +template +const typename EclMaterialLawManagerSimple::MaterialLawParams& EclMaterialLawManagerSimple:: +materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const +{ + using Dir = FaceDir::DirEnum; + if (dirMaterialLawParams_) { + switch(facedir) { + case Dir::XMinus: + case Dir::XPlus: + return dirMaterialLawParams_->materialLawParamsX_[elemIdx]; + case Dir::YMinus: + case Dir::YPlus: + return dirMaterialLawParams_->materialLawParamsY_[elemIdx]; + case Dir::ZMinus: + case Dir::ZPlus: + return dirMaterialLawParams_->materialLawParamsZ_[elemIdx]; + default: + throw std::runtime_error("Unexpected face direction"); + } + } + else { + return materialLawParams_[elemIdx]; + } +} + +// template +// void EclMaterialLawManagerSimple:: +// readGlobalEpsOptions_(const EclipseState& eclState) +// { +// oilWaterEclEpsConfig_ = std::make_shared(); +// oilWaterEclEpsConfig_->initFromState(eclState, EclTwoPhaseSystemType::OilWater); + +// enableEndPointScaling_ = eclState.getTableManager().hasTables("ENKRVD"); +// } + +// template +// void EclMaterialLawManagerSimple:: +// readGlobalHysteresisOptions_(const EclipseState& state) +// { +// hysteresisConfig_ = std::make_shared(); +// hysteresisConfig_->initFromState(state.runspec()); +// } + +template +void EclMaterialLawManagerSimple:: +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 class EclMaterialLawManagerSimple>; +template class EclMaterialLawManagerSimple>; + +} // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp new file mode 100644 index 00000000000..a6ed36ddc4c --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp @@ -0,0 +1,503 @@ +// -*- 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. +*/ +/*! + * \file + * \copydoc Opm::EclMaterialLawManagerSimple + */ +#if ! HAVE_ECL_INPUT +#error "Eclipse input support in opm-common is required to use the ECL material manager!" +#endif + +#ifndef OPM_ECL_MATERIAL_LAW_MANAGER_SIMPLE_HPP +#define OPM_ECL_MATERIAL_LAW_MANAGER_SIMPLE_HPP + +#include +#include +#include "EclTwoPhaseMaterial.hpp" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm { + +class EclipseState; +class EclEpsConfig; +class EclEpsGridProperties; +template class EclEpsScalingPoints; +template struct EclEpsScalingPointsInfo; +class EclHysteresisConfig; +enum class EclTwoPhaseSystemType; +class FieldPropsManager; +class Runspec; +class SgfnTable; +class SgofTable; +class SlgofTable; +class TableColumn; + +/*! + * \ingroup fluidmatrixinteractions + * + * \brief Provides an simple way to create and manage the material law objects + * for a complete ECL deck. + */ +template +class EclMaterialLawManagerSimple +{ +private: + using Traits = TraitsT; + using Scalar = typename Traits::Scalar; + enum { waterPhaseIdx = Traits::wettingPhaseIdx }; + enum { oilPhaseIdx = Traits::nonWettingPhaseIdx }; + enum { gasPhaseIdx = Traits::gasPhaseIdx }; + enum { numPhases = Traits::numPhases }; + + using GasOilTraits = TwoPhaseMaterialTraits; + using OilWaterTraits = TwoPhaseMaterialTraits; + using GasWaterTraits = TwoPhaseMaterialTraits; + + // the two-phase material law which is defined on effective (unscaled) saturations + using GasOilEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial; + using OilWaterEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial; + using GasWaterEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial; + + using GasOilEffectiveTwoPhaseParams = typename GasOilEffectiveTwoPhaseLaw::Params; + using OilWaterEffectiveTwoPhaseParams = typename OilWaterEffectiveTwoPhaseLaw::Params; + using GasWaterEffectiveTwoPhaseParams = typename GasWaterEffectiveTwoPhaseLaw::Params; + + // the two-phase material law which is defined on absolute (scaled) saturations + // using GasOilEpsTwoPhaseLaw = EclEpsTwoPhaseLaw; + // using OilWaterEpsTwoPhaseLaw = EclEpsTwoPhaseLaw; + // using GasWaterEpsTwoPhaseLaw = EclEpsTwoPhaseLaw; + // using GasOilEpsTwoPhaseParams = typename GasOilEpsTwoPhaseLaw::Params; + // using OilWaterEpsTwoPhaseParams = typename OilWaterEpsTwoPhaseLaw::Params; + // using GasWaterEpsTwoPhaseParams = typename GasWaterEpsTwoPhaseLaw::Params; + + // the scaled two-phase material laws with hystersis + // using GasOilTwoPhaseLaw = EclHysteresisTwoPhaseLaw; + // using OilWaterTwoPhaseLaw = EclHysteresisTwoPhaseLaw; + // using GasWaterTwoPhaseLaw = EclHysteresisTwoPhaseLaw; + // using GasOilTwoPhaseHystParams = typename GasOilTwoPhaseLaw::Params; + // using OilWaterTwoPhaseHystParams = typename OilWaterTwoPhaseLaw::Params; + // using GasWaterTwoPhaseHystParams = typename GasWaterTwoPhaseLaw::Params; + + using GasOilTwoPhaseLaw = GasOilEffectiveTwoPhaseLaw; + using OilWaterTwoPhaseLaw = OilWaterEffectiveTwoPhaseLaw; + using GasWaterTwoPhaseLaw = GasWaterEffectiveTwoPhaseLaw; + +public: + // the three-phase material law used by the simulation + // using MaterialLaw = EclMultiplexerMaterial; + using MaterialLaw = EclTwoPhaseMaterial; + using MaterialLawParams = typename MaterialLaw::Params; + using DirectionalMaterialLawParamsPtr = std::unique_ptr>; + + EclMaterialLawManagerSimple(); + ~EclMaterialLawManagerSimple(); + +private: + // internal typedefs + using GasOilEffectiveParamVector = std::vector>; + using OilWaterEffectiveParamVector = std::vector>; + using GasWaterEffectiveParamVector = std::vector>; + + // using GasOilScalingPointsVector = std::vector>>; + // using OilWaterScalingPointsVector = std::vector>>; + // using GasWaterScalingPointsVector = std::vector>>; + using OilWaterScalingInfoVector = std::vector>; + // using GasOilParamVector = std::vector>; + // using OilWaterParamVector = std::vector>; + // using GasWaterParamVector = std::vector>; + using MaterialLawParamsVector = std::vector>; + + // helper classes + + // This class' implementation is defined in "EclMaterialLawManagerSimpleInitParams.cpp" + class InitParams { + public: + InitParams(EclMaterialLawManagerSimple& parent, const EclipseState& eclState, size_t numCompressedElems); + // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup + // field properties of cells on the leaf grid view for CpGrid with local grid refinement. + // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each + // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). + void run(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, + const std::function& lookupIdxOnLevelZeroAssigner); + private: + // class HystParams; + // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup + // field properties of cells on the leaf grid view for CpGrid with local grid refinement. + void copySatnumArrays_(const std::function(const FieldPropsManager&, const std::string&, bool)>& + fieldPropIntOnLeafAssigner); + // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup + // field properties of cells on the leaf grid view for CpGrid with local grid refinement. + void copyIntArray_(std::vector& dest, const std::string keyword, + const std::function(const FieldPropsManager&, const std::string&, bool)>& + fieldPropIntOnLeafAssigner); + unsigned imbRegion_(std::vector& array, unsigned elemIdx); + void initArrays_( + std::vector*>& satnumArray, + std::vector*>& imbnumArray, + std::vector*>& mlpArray); + void initMaterialLawParamVectors_(); + void initOilWaterScaledEpsInfo_(); + // \brief Function argument 'fieldProptOnLeadAssigner' needed to lookup + // field properties of cells on the leaf grid view for CpGrid with local grid refinement. + void initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& + fieldPropIntOnLeafAssigner); + void initThreePhaseParams_( + // HystParams &hystParams, + MaterialLawParams& materialParams, + unsigned satRegionIdx, + unsigned elemIdx); + void readEffectiveParameters_(); + // void readUnscaledEpsPointsVectors_(); + // template + // void readUnscaledEpsPoints_(Container& dest, std::shared_ptr config, EclTwoPhaseSystemType system_type); + unsigned satRegion_(std::vector& array, unsigned elemIdx); + unsigned satOrImbRegion_(std::vector& array, std::vector& default_vec, unsigned elemIdx); +/* + // This class' implementation is defined in "EclMaterialLawManagerSimpleHystParams.cpp" + class HystParams { + public: + HystParams(EclMaterialLawManagerSimple::InitParams& init_params); + void finalize(); + std::shared_ptr getGasOilParams(); + std::shared_ptr getOilWaterParams(); + std::shared_ptr getGasWaterParams(); + void setConfig(unsigned satRegionIdx); + // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each + // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). + void setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner); + void setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner); + void setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner); + void setImbibitionParamsOilGas(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner); + void setImbibitionParamsOilWater(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner); + void setImbibitionParamsGasWater(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner); + private: + bool hasGasWater_(); + bool hasGasOil_(); + bool hasOilWater_(); + + // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each + // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). + std::tuple, EclEpsScalingPoints> + readScaledEpsPoints_(const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, EclTwoPhaseSystemType type, + const std::function& lookupIdxOnLevelZeroAssigner); + std::tuple, EclEpsScalingPoints> + readScaledEpsPointsDrainage_(unsigned elemIdx, EclTwoPhaseSystemType type, + const std::function& lookupIdxOnLevelZeroAssigner); + std::tuple, EclEpsScalingPoints> + readScaledEpsPointsImbibition_(unsigned elemIdx, EclTwoPhaseSystemType type, + const std::function& lookupIdxOnLevelZeroAssigner); + + EclMaterialLawManagerSimple::InitParams& init_params_; + EclMaterialLawManagerSimple& parent_; + const EclipseState& eclState_; + std::shared_ptr gasOilParams_; + std::shared_ptr oilWaterParams_; + std::shared_ptr gasWaterParams_; + }; +*/ + + // This class' implementation is defined in "EclMaterialLawManagerSimpleReadEffectiveParams.cpp" + class ReadEffectiveParams { + public: + ReadEffectiveParams(EclMaterialLawManagerSimple::InitParams& init_params); + void read(); + private: + std::vector normalizeKrValues_(const double tolcrit, const TableColumn& krValues) const; + void readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx); + template + void readGasOilFamily2_( + GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const TableType& sofTable, + const SgfnTable& sgfnTable, + const std::string& columnName); + void readGasOilSgof_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const SgofTable& sgofTable); + + void readGasOilSlgof_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const SlgofTable& slgofTable); + void readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionIdx); + void readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionIdx); + + EclMaterialLawManagerSimple::InitParams& init_params_; + EclMaterialLawManagerSimple& parent_; + const EclipseState& eclState_; + }; // end of "class ReadEffectiveParams" + + EclMaterialLawManagerSimple& parent_; + const EclipseState& eclState_; + size_t numCompressedElems_; + + std::unique_ptr epsImbGridProperties_; //imbibition + std::unique_ptr epsGridProperties_; // drainage + + }; // end of "class InitParams" + +public: + void initFromState(const EclipseState& eclState); + + // \brief Function argument 'fieldPropIntOnLeadAssigner' needed to lookup + // field properties of cells on the leaf grid view for CpGrid with local grid refinement. + // Function argument 'lookupIdxOnLevelZeroAssigner' is added to lookup, for each + // leaf gridview cell with index 'elemIdx', its 'lookupIdx' (index of the parent/equivalent cell on level zero). + void initParamsForElements(const EclipseState& eclState, size_t numCompressedElems, + const std::function(const FieldPropsManager&, const std::string&, bool)>& + fieldPropIntOnLeafAssigner, + const std::function& lookupIdxOnLevelZeroAssigner); + + /*! + * \brief Modify the initial condition according to the SWATINIT keyword. + * + * The method returns the water saturation which yields a givenn capillary + * pressure. The reason this method is not folded directly into initFromState() is + * that the capillary pressure given depends on the particuars of how the simulator + * calculates its initial condition. + */ + std::pair + applySwatinit(unsigned elemIdx, + Scalar pcow, + Scalar Sw); + + /// Apply SWATINIT-like scaling of oil/water capillary pressure curve at + /// simulation restart. + /// + /// \param[in] elemIdx Active cell index + /// + /// \param[in] maxPcow Scaled maximum oil/water capillary pressure. + /// Typically the PPCW restart file array's entry for the + /// corresponding cell. + void applyRestartSwatInit(const unsigned elemIdx, const Scalar maxPcow); + + bool enableEndPointScaling() const + { return enableEndPointScaling_; } + + bool enablePpcwmax() const + { return enablePpcwmax_; } + + bool enableHysteresis() const + //{ return hysteresisConfig_->enableHysteresis(); } + { return false; } + + bool enablePCHysteresis() const + //{ return (enableHysteresis() && hysteresisConfig_->pcHysteresisModel() >= 0); } + { return false; } + + bool enableWettingHysteresis() const + //{ return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 4); } + { return false; } + + bool enableNonWettingHysteresis() const + //{ return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 0); } + { return false; } + + MaterialLawParams& materialLawParams(unsigned elemIdx) + { + assert(elemIdx < materialLawParams_.size()); + return materialLawParams_[elemIdx]; + } + + const MaterialLawParams& materialLawParams(unsigned elemIdx) const + { + assert(elemIdx < materialLawParams_.size()); + return materialLawParams_[elemIdx]; + } + + const MaterialLawParams& materialLawParams(unsigned elemIdx, FaceDir::DirEnum facedir) const + { + return materialLawParamsFunc_(elemIdx, facedir); + } + + MaterialLawParams& materialLawParams(unsigned elemIdx, FaceDir::DirEnum facedir) + { + return const_cast(materialLawParamsFunc_(elemIdx, facedir)); + } + + /*! + * \brief Returns a material parameter object for a given element and saturation region. + * + * This method changes the saturation table idx in the original material law parameter object. + * In the context of ECL reservoir simulators, this is required to properly handle + * 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; + + int satnumRegionIdx(unsigned elemIdx) const + { return satnumRegionArray_[elemIdx]; } + + int getKrnumSatIdx(unsigned elemIdx, FaceDir::DirEnum facedir) const; + + bool hasDirectionalRelperms() const + { + return !krnumXArray_.empty() || !krnumYArray_.empty() || !krnumZArray_.empty(); + } + + bool hasDirectionalImbnum() const { + if (imbnumXArray_.size() > 0 || imbnumYArray_.size() > 0 || imbnumZArray_.size() > 0) { + return true; + } + return false; + } + + int imbnumRegionIdx(unsigned elemIdx) const + { return imbnumRegionArray_[elemIdx]; } + + template + bool updateHysteresis(const FluidState& fluidState, unsigned elemIdx) + { + OPM_TIMEFUNCTION_LOCAL(); + if (!enableHysteresis()) + return false; + bool changed = MaterialLaw::updateHysteresis(materialLawParams(elemIdx), fluidState); + if (hasDirectionalRelperms() || hasDirectionalImbnum()) { + using Dir = FaceDir::DirEnum; + constexpr int ndim = 3; + Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus}; + for (int i = 0; i& oilWaterScaledEpsPointsDrainage(unsigned elemIdx); + + const EclEpsScalingPointsInfo& oilWaterScaledEpsInfoDrainage(size_t elemIdx) const + { return oilWaterScaledEpsInfoDrainage_[elemIdx]; } + + template + void serializeOp(Serializer& serializer) + { + // This is for restart serialization. + // Only dynamic state in the parameters need to be stored. + // For that reason we do not serialize the vector + // as that would recreate the objects inside. + for (auto& mat : materialLawParams_) { + serializer(mat); + } + } + +private: + const MaterialLawParams& materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const; + + // void readGlobalEpsOptions_(const EclipseState& eclState); + + // void readGlobalHysteresisOptions_(const EclipseState& state); + + void readGlobalThreePhaseOptions_(const Runspec& runspec); + + bool enableEndPointScaling_; + // std::shared_ptr hysteresisConfig_; + // std::vector> wagHystersisConfig_; + + + std::shared_ptr oilWaterEclEpsConfig_; + std::vector> unscaledEpsInfo_; + OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_; + + // std::shared_ptr gasWaterEclEpsConfig_; + + // GasOilScalingPointsVector gasOilUnscaledPointsVector_; + // OilWaterScalingPointsVector oilWaterUnscaledPointsVector_; + // GasWaterScalingPointsVector gasWaterUnscaledPointsVector_; + + GasOilEffectiveParamVector gasOilEffectiveParamVector_; + OilWaterEffectiveParamVector oilWaterEffectiveParamVector_; + GasWaterEffectiveParamVector gasWaterEffectiveParamVector_; + + EclMultiplexerApproach threePhaseApproach_ = EclMultiplexerApproach::Default; + // this attribute only makes sense for twophase simulations! + enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::GasWater; + + std::vector materialLawParams_; + DirectionalMaterialLawParamsPtr dirMaterialLawParams_; + + std::vector satnumRegionArray_; + std::vector krnumXArray_; + std::vector krnumYArray_; + std::vector krnumZArray_; + std::vector imbnumXArray_; + std::vector imbnumYArray_; + std::vector imbnumZArray_; + std::vector imbnumRegionArray_; + std::vector stoneEtas_; + + bool enablePpcwmax_; + std::vector maxAllowPc_; + std::vector modifySwl_; + + bool hasGas; + bool hasOil; + bool hasWater; + + std::shared_ptr gasOilConfig_; + std::shared_ptr oilWaterConfig_; + std::shared_ptr gasWaterConfig_; +}; +} // namespace Opm + +#endif diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp new file mode 100644 index 00000000000..b9df445069a --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp @@ -0,0 +1,314 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2022 Equinor ASA. + + 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 3 of the License, or + + 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 . +*/ + +#if 0 + +#include +#include +#include + +namespace Opm { + +/* constructors*/ +template +EclMaterialLawManagerSimple::InitParams::HystParams:: +HystParams(EclMaterialLawManagerSimple::InitParams& init_params) : + init_params_{init_params}, parent_{init_params_.parent_}, + eclState_{init_params_.eclState_} +{ + gasOilParams_ = std::make_shared(); + oilWaterParams_ = std::make_shared(); + gasWaterParams_ = std::make_shared(); +} + +/* public methods, alphabetically sorted */ + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +finalize() +{ + if (hasGasOil_()) + this->gasOilParams_->finalize(); + if (hasOilWater_()) + this->oilWaterParams_->finalize(); + if (hasGasWater_()) + this->gasWaterParams_->finalize(); +} + +template +std::shared_ptr::GasOilTwoPhaseHystParams> +EclMaterialLawManagerSimple::InitParams::HystParams:: +getGasOilParams() +{ + return gasOilParams_; +} + +template +std::shared_ptr::OilWaterTwoPhaseHystParams> +EclMaterialLawManagerSimple::InitParams::HystParams:: +getOilWaterParams() +{ + return oilWaterParams_; +} + +template +std::shared_ptr::GasWaterTwoPhaseHystParams> +EclMaterialLawManagerSimple::InitParams::HystParams:: +getGasWaterParams() +{ + return gasWaterParams_; +} + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setConfig(unsigned satRegionIdx) +{ + this->gasOilParams_->setConfig(this->parent_.hysteresisConfig_); + this->oilWaterParams_->setConfig(this->parent_.hysteresisConfig_); + this->gasWaterParams_->setConfig(this->parent_.hysteresisConfig_); + + if (this->parent_.hysteresisConfig_->enableWagHysteresis()) { + this->gasOilParams_->setWagConfig(this->parent_.wagHystersisConfig_[satRegionIdx]); + this->oilWaterParams_->setWagConfig(this->parent_.wagHystersisConfig_[satRegionIdx]); + this->gasWaterParams_->setWagConfig(this->parent_.wagHystersisConfig_[satRegionIdx]); + } + +} // namespace Opm + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + if (hasGasWater_()) { + auto [gasWaterScaledInfo, gasWaterScaledPoints] + = readScaledEpsPointsDrainage_(elemIdx, EclTwoPhaseSystemType::GasWater, lookupIdxOnLevelZeroAssigner); + GasWaterEpsTwoPhaseParams gasWaterDrainParams; + gasWaterDrainParams.setConfig(this->parent_.gasWaterConfig_); + gasWaterDrainParams.setUnscaledPoints(this->parent_.gasWaterUnscaledPointsVector_[satRegionIdx]); + gasWaterDrainParams.setScaledPoints(gasWaterScaledPoints); + gasWaterDrainParams.setEffectiveLawParams(this->parent_.gasWaterEffectiveParamVector_[satRegionIdx]); + gasWaterDrainParams.finalize(); + this->gasWaterParams_->setDrainageParams(gasWaterDrainParams, gasWaterScaledInfo, EclTwoPhaseSystemType::GasWater); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + if (hasGasOil_()) { + auto [gasOilScaledInfo, gasOilScaledPoints] + = readScaledEpsPointsDrainage_(elemIdx, EclTwoPhaseSystemType::GasOil, lookupIdxOnLevelZeroAssigner); + GasOilEpsTwoPhaseParams gasOilDrainParams; + gasOilDrainParams.setConfig(this->parent_.gasOilConfig_); + gasOilDrainParams.setUnscaledPoints(this->parent_.gasOilUnscaledPointsVector_[satRegionIdx]); + gasOilDrainParams.setScaledPoints(gasOilScaledPoints); + gasOilDrainParams.setEffectiveLawParams(this->parent_.gasOilEffectiveParamVector_[satRegionIdx]); + gasOilDrainParams.finalize(); + this->gasOilParams_->setDrainageParams(gasOilDrainParams, gasOilScaledInfo, EclTwoPhaseSystemType::GasOil); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + // We need to compute the oil-water scaled info even if we are running a two-phase case without + // water (e.g. gas-oil). The reason is that the oil-water scaled info is used when computing + // the initial condition see e.g. equilibrationhelpers.cc and initstateequil.cc + // Therefore, the below 7 lines should not be put inside the if(hasOilWater_){} below. + auto [oilWaterScaledInfo, oilWaterScaledPoints] + = readScaledEpsPointsDrainage_(elemIdx, EclTwoPhaseSystemType::OilWater, lookupIdxOnLevelZeroAssigner); + // TODO: This will reassign the same EclEpsScalingPointsInfo for each facedir + // since we currently does not support facedir for the scaling points info + // When such support is added, we need to extend the below vector which has info for each cell + // to include three more vectors, one with info for each facedir of a cell + this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx] = oilWaterScaledInfo; + if (hasOilWater_()) { + OilWaterEpsTwoPhaseParams oilWaterDrainParams; + oilWaterDrainParams.setConfig(this->parent_.oilWaterConfig_); + oilWaterDrainParams.setUnscaledPoints(this->parent_.oilWaterUnscaledPointsVector_[satRegionIdx]); + oilWaterDrainParams.setScaledPoints(oilWaterScaledPoints); + oilWaterDrainParams.setEffectiveLawParams(this->parent_.oilWaterEffectiveParamVector_[satRegionIdx]); + oilWaterDrainParams.finalize(); + oilWaterParams_->setDrainageParams(oilWaterDrainParams, oilWaterScaledInfo, EclTwoPhaseSystemType::OilWater); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setImbibitionParamsGasWater(unsigned elemIdx, unsigned imbRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + if (hasGasWater_()) { + auto [gasWaterScaledInfo, gasWaterScaledPoints] + = readScaledEpsPointsImbibition_(elemIdx, EclTwoPhaseSystemType::GasWater, lookupIdxOnLevelZeroAssigner); + GasWaterEpsTwoPhaseParams gasWaterImbParamsHyst; + gasWaterImbParamsHyst.setConfig(this->parent_.gasWaterConfig_); + gasWaterImbParamsHyst.setUnscaledPoints(this->parent_.gasWaterUnscaledPointsVector_[imbRegionIdx]); + gasWaterImbParamsHyst.setScaledPoints(gasWaterScaledPoints); + gasWaterImbParamsHyst.setEffectiveLawParams(this->parent_.gasWaterEffectiveParamVector_[imbRegionIdx]); + gasWaterImbParamsHyst.finalize(); + this->gasWaterParams_->setImbibitionParams(gasWaterImbParamsHyst, + gasWaterScaledInfo, + EclTwoPhaseSystemType::GasWater); + + } +} + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setImbibitionParamsOilGas(unsigned elemIdx, unsigned imbRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + if (hasGasOil_()) { + auto [gasOilScaledInfo, gasOilScaledPoints] + = readScaledEpsPointsImbibition_(elemIdx, EclTwoPhaseSystemType::GasOil, lookupIdxOnLevelZeroAssigner); + + GasOilEpsTwoPhaseParams gasOilImbParamsHyst; + gasOilImbParamsHyst.setConfig(this->parent_.gasOilConfig_); + gasOilImbParamsHyst.setUnscaledPoints(this->parent_.gasOilUnscaledPointsVector_[imbRegionIdx]); + gasOilImbParamsHyst.setScaledPoints(gasOilScaledPoints); + gasOilImbParamsHyst.setEffectiveLawParams(this->parent_.gasOilEffectiveParamVector_[imbRegionIdx]); + gasOilImbParamsHyst.finalize(); + this->gasOilParams_->setImbibitionParams(gasOilImbParamsHyst, + gasOilScaledInfo, + EclTwoPhaseSystemType::GasOil); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams::HystParams:: +setImbibitionParamsOilWater(unsigned elemIdx, unsigned imbRegionIdx, + const std::function& lookupIdxOnLevelZeroAssigner) +{ + if (hasOilWater_()) { + auto [oilWaterScaledInfo, oilWaterScaledPoints] + = readScaledEpsPointsImbibition_(elemIdx, EclTwoPhaseSystemType::OilWater, lookupIdxOnLevelZeroAssigner); + OilWaterEpsTwoPhaseParams oilWaterImbParamsHyst; + oilWaterImbParamsHyst.setConfig(this->parent_.oilWaterConfig_); + oilWaterImbParamsHyst.setUnscaledPoints(this->parent_.oilWaterUnscaledPointsVector_[imbRegionIdx]); + oilWaterImbParamsHyst.setScaledPoints(oilWaterScaledPoints); + oilWaterImbParamsHyst.setEffectiveLawParams(this->parent_.oilWaterEffectiveParamVector_[imbRegionIdx]); + oilWaterImbParamsHyst.finalize(); + this->oilWaterParams_->setImbibitionParams(oilWaterImbParamsHyst, + oilWaterScaledInfo, + EclTwoPhaseSystemType::OilWater); + + } +} + +/* private methods, alphabetically sorted */ + +template +bool +EclMaterialLawManagerSimple::InitParams::HystParams:: +hasGasOil_() +{ + return this->parent_.hasGas && this->parent_.hasOil; +} + +template +bool +EclMaterialLawManagerSimple::InitParams::HystParams:: +hasGasWater_() +{ + return this->parent_.hasGas && this->parent_.hasWater && !this->parent_.hasOil; +} + +template +bool +EclMaterialLawManagerSimple::InitParams::HystParams:: +hasOilWater_() +{ + return this->parent_.hasOil && this->parent_.hasWater; +} + +template +std::tuple< + EclEpsScalingPointsInfo::Scalar>, + EclEpsScalingPoints::Scalar> +> +EclMaterialLawManagerSimple::InitParams::HystParams:: +readScaledEpsPoints_(const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, EclTwoPhaseSystemType type, + const std::function& fieldPropIdxOnLevelZero) +{ + const EclEpsConfig& config = (type == EclTwoPhaseSystemType::OilWater)? *(this->parent_.oilWaterConfig_): *(this->parent_.gasOilConfig_); + // For CpGrids with LGRs, field prop is inherited from parent/equivalent cell from level 0. + // 'lookupIdx' is the index on level zero of the parent cell or the equivalent cell of the + // leaf grid view cell with index 'elemIdx'. + const auto lookupIdx = fieldPropIdxOnLevelZero(elemIdx); + unsigned satRegionIdx = epsGridProperties.satRegion( lookupIdx /* coincides with elemIdx when no LGRs */ ); + // Copy-construct a new instance of EclEpsScalingPointsInfo + EclEpsScalingPointsInfo destInfo(this->parent_.unscaledEpsInfo_[satRegionIdx]); + // TODO: currently epsGridProperties does not implement a face direction, e.g. SWLX, SWLY,... + // when these keywords get implemented, we need to use include facedir in the lookup + destInfo.extractScaled(this->eclState_, epsGridProperties, lookupIdx /* coincides with elemIdx when no LGRs */); + + EclEpsScalingPoints destPoint; + destPoint.init(destInfo, config, type); + return {destInfo, destPoint}; +} + +template +std::tuple< + EclEpsScalingPointsInfo::Scalar>, + EclEpsScalingPoints::Scalar> +> +EclMaterialLawManagerSimple::InitParams::HystParams:: +readScaledEpsPointsDrainage_(unsigned elemIdx, EclTwoPhaseSystemType type, + const std::function& fieldPropIdxOnLevelZero) +{ + const auto& epsGridProperties = this->init_params_.epsGridProperties_; + return readScaledEpsPoints_(*epsGridProperties, elemIdx, type, fieldPropIdxOnLevelZero); +} + +template +std::tuple< + EclEpsScalingPointsInfo::Scalar>, + EclEpsScalingPoints::Scalar> +> +EclMaterialLawManagerSimple::InitParams::HystParams:: +readScaledEpsPointsImbibition_(unsigned elemIdx, EclTwoPhaseSystemType type, + const std::function& fieldPropIdxOnLevelZero) +{ + const auto& epsGridProperties = this->init_params_.epsImbGridProperties_; + return readScaledEpsPoints_(*epsGridProperties, elemIdx, type, fieldPropIdxOnLevelZero); +} + +// Make some actual code, by realizing the previously defined templated class +template class EclMaterialLawManagerSimple>::InitParams::HystParams; +template class EclMaterialLawManagerSimple>::InitParams::HystParams; + + +} // namespace Opm + +#endif // 0 diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp new file mode 100644 index 00000000000..0ce62ed1230 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp @@ -0,0 +1,300 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2022 Equinor ASA. + + 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 3 of the License, or + + 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 . +*/ + +#include + +#include + +#include +#include + + +namespace Opm { + +/* constructors*/ + +template +EclMaterialLawManagerSimple::InitParams:: +InitParams(EclMaterialLawManagerSimple& parent, const EclipseState& eclState, size_t numCompressedElems) : + parent_{parent}, + eclState_{eclState}, + numCompressedElems_{numCompressedElems} +{ + // read end point scaling grid properties + // TODO: these objects might require some memory, can this be simplified? + if (this->parent_.enableHysteresis()) { + this->epsImbGridProperties_ + = std::make_unique(this->eclState_, /*useImbibition=*/true); + } + this->epsGridProperties_ + = std::make_unique(this->eclState_, /*useImbibition=*/false); +} + +/* public methods */ + +template +void +EclMaterialLawManagerSimple::InitParams:: +run(const std::function(const FieldPropsManager&, const std::string&, bool)>& + fieldPropIntOnLeafAssigner, + const std::function& lookupIdxOnLevelZeroAssigner) { + // readUnscaledEpsPointsVectors_(); + readEffectiveParameters_(); + initSatnumRegionArray_(fieldPropIntOnLeafAssigner); + copySatnumArrays_(fieldPropIntOnLeafAssigner); + initOilWaterScaledEpsInfo_(); + initMaterialLawParamVectors_(); + std::vector*> satnumArray; + std::vector*> imbnumArray; + std::vector*> mlpArray; + initArrays_(satnumArray, imbnumArray, mlpArray); + auto num_arrays = mlpArray.size(); + for (unsigned i=0; inumCompressedElems_; ++elemIdx) { + unsigned satRegionIdx = satRegion_(*satnumArray[i], elemIdx); + //unsigned satNumCell = this->parent_.satnumRegionArray_[elemIdx]; + // HystParams hystParams {*this}; + // hystParams.setConfig(satRegionIdx); + // hystParams.setDrainageParamsOilGas(elemIdx, satRegionIdx, lookupIdxOnLevelZeroAssigner); + // hystParams.setDrainageParamsOilWater(elemIdx, satRegionIdx, lookupIdxOnLevelZeroAssigner); + // hystParams.setDrainageParamsGasWater(elemIdx, satRegionIdx, lookupIdxOnLevelZeroAssigner); + // if (this->parent_.enableHysteresis()) { + // unsigned imbRegionIdx = imbRegion_(*imbnumArray[i], elemIdx); + // hystParams.setImbibitionParamsOilGas(elemIdx, imbRegionIdx, lookupIdxOnLevelZeroAssigner); + // hystParams.setImbibitionParamsOilWater(elemIdx, imbRegionIdx, lookupIdxOnLevelZeroAssigner); + // hystParams.setImbibitionParamsGasWater(elemIdx, imbRegionIdx, lookupIdxOnLevelZeroAssigner); + // } + // hystParams.finalize(); + initThreePhaseParams_(/*hystParams,*/ (*mlpArray[i])[elemIdx], satRegionIdx, elemIdx); + } + } +} + +/* private methods alphabetically sorted*/ + +template +void +EclMaterialLawManagerSimple::InitParams:: +copySatnumArrays_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) +{ + copyIntArray_(this->parent_.krnumXArray_, "KRNUMX", fieldPropIntOnLeafAssigner); + copyIntArray_(this->parent_.krnumYArray_, "KRNUMY", fieldPropIntOnLeafAssigner); + copyIntArray_(this->parent_.krnumZArray_, "KRNUMZ", fieldPropIntOnLeafAssigner); + copyIntArray_(this->parent_.imbnumXArray_, "IMBNUMX", fieldPropIntOnLeafAssigner); + copyIntArray_(this->parent_.imbnumYArray_, "IMBNUMY", fieldPropIntOnLeafAssigner); + copyIntArray_(this->parent_.imbnumZArray_, "IMBNUMZ", fieldPropIntOnLeafAssigner); + // create the information for the imbibition region (IMBNUM). By default this is + // the same as the saturation region (SATNUM) + this->parent_.imbnumRegionArray_ = this->parent_.satnumRegionArray_; + copyIntArray_(this->parent_.imbnumRegionArray_, "IMBNUM", fieldPropIntOnLeafAssigner); + assert(this->numCompressedElems_ == this->parent_.satnumRegionArray_.size()); + assert(!this->parent_.enableHysteresis() || this->numCompressedElems_ == this->parent_.imbnumRegionArray_.size()); +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +copyIntArray_(std::vector& dest, const std::string keyword, + const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) +{ + if (this->eclState_.fieldProps().has_int(keyword)) { + dest = fieldPropIntOnLeafAssigner(this->eclState_.fieldProps(), keyword, /*needsTranslation*/true); + } +} + +template +unsigned +EclMaterialLawManagerSimple::InitParams:: +imbRegion_(std::vector& array, unsigned elemIdx) +{ + std::vector& default_vec = this->parent_.imbnumRegionArray_; + return satOrImbRegion_(array, default_vec, elemIdx); +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +initArrays_( + std::vector*>& satnumArray, + std::vector*>& imbnumArray, + std::vector*>& mlpArray) +{ + satnumArray.push_back(&this->parent_.satnumRegionArray_); + imbnumArray.push_back(&this->parent_.imbnumRegionArray_); + mlpArray.push_back(&this->parent_.materialLawParams_); + if (this->parent_.dirMaterialLawParams_) { + if (this->parent_.hasDirectionalRelperms()) { + satnumArray.push_back(&this->parent_.krnumXArray_); + satnumArray.push_back(&this->parent_.krnumYArray_); + satnumArray.push_back(&this->parent_.krnumZArray_); + } + if (this->parent_.hasDirectionalImbnum()) { + imbnumArray.push_back(&this->parent_.imbnumXArray_); + imbnumArray.push_back(&this->parent_.imbnumYArray_); + imbnumArray.push_back(&this->parent_.imbnumZArray_); + } + mlpArray.push_back(&(this->parent_.dirMaterialLawParams_->materialLawParamsX_)); + mlpArray.push_back(&(this->parent_.dirMaterialLawParams_->materialLawParamsY_)); + mlpArray.push_back(&(this->parent_.dirMaterialLawParams_->materialLawParamsZ_)); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +initMaterialLawParamVectors_() +{ + this->parent_.materialLawParams_.resize(this->numCompressedElems_); + if (this->parent_.hasDirectionalImbnum() || this->parent_.hasDirectionalRelperms()) { + this->parent_.dirMaterialLawParams_ + = std::make_unique>(this->numCompressedElems_); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +initOilWaterScaledEpsInfo_() +{ + // This vector will be updated in the hystParams.setDrainageOilWater() in the run() method + this->parent_.oilWaterScaledEpsInfoDrainage_.resize(this->numCompressedElems_); +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner) +{ + // copy the SATNUM grid property. in some cases this is not necessary, but it + // should not require much memory anyway... + auto &satnumArray = this->parent_.satnumRegionArray_; + satnumArray.resize(this->numCompressedElems_); + if (this->eclState_.fieldProps().has_int("SATNUM")) { + satnumArray = fieldPropIntOnLeafAssigner(this->eclState_.fieldProps(), "SATNUM", /*needsTranslation*/true); + } + else { + std::fill(satnumArray.begin(), satnumArray.end(), 0); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +initThreePhaseParams_(//HystParams &hystParams, + MaterialLawParams& materialParams, + unsigned satRegionIdx, + [[maybe_unused]] unsigned elemIdx) +{ + auto& realParams = materialParams; + realParams.setApproach(this->parent_.twoPhaseApproach_); + + if (realParams.approach() == EclTwoPhaseApproach::GasOil) { + realParams.setGasOilParams(this->parent_.gasOilEffectiveParamVector_[satRegionIdx]); + } + else if (realParams.approach() == EclTwoPhaseApproach::GasWater) { + realParams.setGasWaterParams(this->parent_.gasWaterEffectiveParamVector_[satRegionIdx]); + } + else if (realParams.approach() == EclTwoPhaseApproach::OilWater) { + realParams.setOilWaterParams(this->parent_.oilWaterEffectiveParamVector_[satRegionIdx]); + } + realParams.finalize(); +} + +template +void +EclMaterialLawManagerSimple::InitParams:: +readEffectiveParameters_() +{ + ReadEffectiveParams effectiveReader {*this}; + // populates effective parameter vectors in the parent class (EclMaterialManager) + effectiveReader.read(); +} + +// template +// void +// EclMaterialLawManagerSimple::InitParams:: +// readUnscaledEpsPointsVectors_() +// { +// if (this->parent_.hasGas && this->parent_.hasOil) { +// readUnscaledEpsPoints_( +// this->parent_.gasOilUnscaledPointsVector_, +// this->parent_.gasOilConfig_, +// EclTwoPhaseSystemType::GasOil +// ); +// } +// if (this->parent_.hasOil && this->parent_.hasWater) { +// readUnscaledEpsPoints_( +// this->parent_.oilWaterUnscaledPointsVector_, +// this->parent_.oilWaterConfig_, +// EclTwoPhaseSystemType::OilWater +// ); +// } +// if (!this->parent_.hasOil) { +// readUnscaledEpsPoints_( +// this->parent_.gasWaterUnscaledPointsVector_, +// this->parent_.gasWaterConfig_, +// EclTwoPhaseSystemType::GasWater +// ); +// } +// } + +// template +// template +// void +// EclMaterialLawManagerSimple::InitParams:: +// readUnscaledEpsPoints_(Container& dest, std::shared_ptr config, EclTwoPhaseSystemType system_type) +// { +// const size_t numSatRegions = this->eclState_.runspec().tabdims().getNumSatTables(); +// dest.resize(numSatRegions); +// for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { +// dest[satRegionIdx] = std::make_shared >(); +// dest[satRegionIdx]->init(this->parent_.unscaledEpsInfo_[satRegionIdx], *config, system_type); +// } +// } + +template +unsigned +EclMaterialLawManagerSimple::InitParams:: +satRegion_(std::vector& array, unsigned elemIdx) +{ + std::vector& default_vec = this->parent_.satnumRegionArray_; + return satOrImbRegion_(array, default_vec, elemIdx); +} + +template +unsigned +EclMaterialLawManagerSimple::InitParams:: +satOrImbRegion_(std::vector& array, std::vector& default_vec, unsigned elemIdx) +{ + int value; + if (array.size() > 0) { + value = array[elemIdx]; + } + else { // use default value + value = default_vec[elemIdx]; + } + return static_cast(value); +} + +// Make some actual code, by realizing the previously defined templated class +template class EclMaterialLawManagerSimple>::InitParams; +template class EclMaterialLawManagerSimple>::InitParams; + +} // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp new file mode 100644 index 00000000000..417d845a39f --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp @@ -0,0 +1,380 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2022 Equinor ASA. + + 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 3 of the License, or + + 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 . +*/ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace Opm { + +/* constructors*/ +template +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +ReadEffectiveParams(EclMaterialLawManagerSimple::InitParams& init_params) : + init_params_{init_params}, parent_{init_params_.parent_}, + eclState_{init_params_.eclState_} +{ +} + +/* public methods */ +template +void +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +read() { + auto& gasOilVector = this->parent_.gasOilEffectiveParamVector_; + auto& oilWaterVector = this->parent_.oilWaterEffectiveParamVector_; + auto& gasWaterVector = this->parent_.gasWaterEffectiveParamVector_; + const size_t numSatRegions = this->eclState_.runspec().tabdims().getNumSatTables(); + gasOilVector.resize(numSatRegions); + oilWaterVector.resize(numSatRegions); + gasWaterVector.resize(numSatRegions); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { + readGasOilParameters_(gasOilVector, satRegionIdx); + readOilWaterParameters_(oilWaterVector, satRegionIdx); + readGasWaterParameters_(gasWaterVector, satRegionIdx); + } + +} + +/* private methods, alphabetically sorted*/ + +// Relative permeability values not strictly greater than 'tolcrit' treated as zero. +template +std::vector +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +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; +} + +template +void +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) +{ + if (!this->parent_.hasGas || !this->parent_.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 = this->parent_.unscaledEpsInfo_[satRegionIdx].Swl; + const auto tolcrit = this->eclState_.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto& tableManager = this->eclState_.getTableManager(); + + switch (this->eclState_.runspec().saturationFunctionControls().family()) { + case SatFuncControls::KeywordFamily::Family_I: + { + const TableContainer& sgofTables = tableManager.getSgofTables(); + const TableContainer& slgofTables = tableManager.getSlgofTables(); + if (!sgofTables.empty()) + readGasOilSgof_(effParams, Swco, tolcrit, sgofTables.template getTable(satRegionIdx)); + else if (!slgofTables.empty()) + readGasOilSlgof_(effParams, Swco, tolcrit, slgofTables.template getTable(satRegionIdx)); + else if ( !tableManager.getSgofletTable().empty() ) { + throw std::runtime_error("LET tables not supported!"); + } + break; + } + + case SatFuncControls::KeywordFamily::Family_II: + { + const SgfnTable& sgfnTable = tableManager.getSgfnTables().template getTable( satRegionIdx ); + if (!this->parent_.hasWater) { + // oil and gas case + const Sof2Table& sof2Table = tableManager.getSof2Tables().template getTable( satRegionIdx ); + readGasOilFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable, /*columnName=*/"KRO"); + } + else { + const Sof3Table& sof3Table = tableManager.getSof3Tables().template getTable( satRegionIdx ); + readGasOilFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable, /* columnName=*/"KROG"); + } + break; + } + + case SatFuncControls::KeywordFamily::Family_III: + { + throw std::domain_error("Saturation keyword family III is not applicable for a gas-oil system"); + } + + case SatFuncControls::KeywordFamily::Undefined: + throw std::domain_error("No valid saturation keyword family specified"); + } +} + +template +template +void +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +readGasOilFamily2_(GasOilEffectiveTwoPhaseParams& effParams, + const Scalar Swco, + const double tolcrit, + const TableType& sofTable, + const SgfnTable& sgfnTable, + const std::string& columnName) +{ + // convert the saturations of the SGFN keyword from gas to oil saturations + std::vector SoSamples(sgfnTable.numRows()); + std::vector SoColumn = sofTable.getColumn("SO").vectorCopy(); + for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) { + SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); + } + + auto& realParams = effParams; + + realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sofTable.getColumn(columnName))); + realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); + realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy()); + realParams.finalize(); +} + +template +void +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +readGasOilSgof_(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); + } + + auto& realParams = effParams; + + 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 +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +readGasOilSlgof_(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; + } + + auto& realParams = effParams; + + 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 +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionIdx) +{ + if (!this->parent_.hasGas || !this->parent_.hasWater || this->parent_.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 = this->eclState_.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto& tableManager = this->eclState_.getTableManager(); + + switch (this->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: + { + const TableContainer& sgwfnTables = tableManager.getSgwfnTables(); + auto& realParams = effParams; + if (!sgwfnTables.empty()){ + const SgwfnTable& sgwfnTable = tableManager.getSgwfnTables().template getTable( satRegionIdx ); + std::vector SwSamples(sgwfnTable.numRows()); + for (size_t sampleIdx = 0; sampleIdx < sgwfnTable.numRows(); ++ sampleIdx) + SwSamples[sampleIdx] = 1 - sgwfnTable.get("SG", sampleIdx); + realParams.setKrwSamples(SwSamples, normalizeKrValues_(tolcrit, sgwfnTable.getColumn("KRGW"))); + realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgwfnTable.getColumn("KRG"))); + realParams.setPcnwSamples(SwSamples, sgwfnTable.getColumn("PCGW").vectorCopy()); + } + else { + const SgfnTable& sgfnTable = tableManager.getSgfnTables().template getTable( satRegionIdx ); + const SwfnTable& swfnTable = tableManager.getSwfnTables().template getTable( satRegionIdx ); + + 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::Family_III: + { + const GsfTable& gsfTable = tableManager.getGsfTables().template getTable( satRegionIdx ); + const WsfTable& wsfTable = tableManager.getWsfTables().template getTable( satRegionIdx ); + + auto& realParams = effParams; + + std::vector SwColumn = wsfTable.getColumn("SW").vectorCopy(); + + realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, wsfTable.getColumn("KRW"))); + std::vector SwSamples(gsfTable.numRows()); + for (size_t sampleIdx = 0; sampleIdx < gsfTable.numRows(); ++ sampleIdx) + SwSamples[sampleIdx] = 1 - gsfTable.get("SG", sampleIdx); + realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, gsfTable.getColumn("KRG"))); + //Capillary pressure is read from GSF. + realParams.setPcnwSamples(SwSamples, gsfTable.getColumn("PCGW").vectorCopy()); + realParams.finalize(); + + break; + } + case SatFuncControls::KeywordFamily::Undefined: + throw std::domain_error("No valid saturation keyword family specified"); + } +} + +template +void +EclMaterialLawManagerSimple::InitParams::ReadEffectiveParams:: +readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionIdx) +{ + if (!this->parent_.hasOil || !this->parent_.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 = this->eclState_.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); + + const auto& tableManager = this->eclState_.getTableManager(); + auto& effParams = *dest[satRegionIdx]; + + switch (this->eclState_.runspec().saturationFunctionControls().family()) { + case SatFuncControls::KeywordFamily::Family_I: + { + if (tableManager.hasTables("SWOF")) { + const auto& swofTable = tableManager.getSwofTables().template getTable(satRegionIdx); + const std::vector SwColumn = swofTable.getColumn("SW").vectorCopy(); + + auto& realParams = effParams; + + 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() ) { + throw std::runtime_error("LET tables not supported!"); + } + break; + } + + case SatFuncControls::KeywordFamily::Family_II: + { + const auto& swfnTable = tableManager.getSwfnTables().template getTable(satRegionIdx); + const std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); + + auto& realParams = effParams; + + realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); + realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy()); + + if (!this->parent_.hasGas) { + const auto& sof2Table = tableManager.getSof2Tables().template 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().template 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::Family_III: + { + throw std::domain_error("Saturation keyword family III is not applicable for a oil-water system"); + } + + case SatFuncControls::KeywordFamily::Undefined: + throw std::domain_error("No valid saturation keyword family specified"); + } +} + +// Make some actual code, by realizing the previously defined templated class +template class EclMaterialLawManagerSimple>::InitParams::ReadEffectiveParams; +template class EclMaterialLawManagerSimple>::InitParams::ReadEffectiveParams; + + +} // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp index 1042fcf957e..29d771b9892 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp @@ -33,6 +33,7 @@ #include #include +#include namespace Opm { @@ -111,6 +112,11 @@ class EclTwoPhaseMaterial : public TraitsT //! are dependent on the phase composition static constexpr bool isCompositionDependent = false; + static constexpr bool isHysteresisDependent = + GasOilMaterialLaw::isHysteresisDependent + || OilWaterMaterialLaw::isHysteresisDependent + || GasWaterMaterialLaw::isHysteresisDependent; + template static Scalar relpermOilInOilGasSystem(const Params& /*params*/, const FluidState& /*fluidState*/) { @@ -192,12 +198,14 @@ class EclTwoPhaseMaterial : public TraitsT Scalar& swMin, const Params& params) { - soMax = 1.0 - params.oilWaterParams().krnSwMdc(); - swMax = params.oilWaterParams().krwSwMdc(); - swMin = params.oilWaterParams().pcSwMdc(); - Valgrind::CheckDefined(soMax); - Valgrind::CheckDefined(swMax); - Valgrind::CheckDefined(swMin); + if constexpr (isHysteresisDependent) { + soMax = 1.0 - params.oilWaterParams().krnSwMdc(); + swMax = params.oilWaterParams().krwSwMdc(); + swMin = params.oilWaterParams().pcSwMdc(); + Valgrind::CheckDefined(soMax); + Valgrind::CheckDefined(swMax); + Valgrind::CheckDefined(swMin); + } } /* @@ -212,7 +220,9 @@ class EclTwoPhaseMaterial : public TraitsT const Scalar& swMin, Params& params) { - params.oilWaterParams().update(swMin, swMax, 1.0 - soMax); + if constexpr (isHysteresisDependent) { + params.oilWaterParams().update(swMin, swMax, 1.0 - soMax); + } } @@ -228,12 +238,14 @@ class EclTwoPhaseMaterial : public TraitsT Scalar& somin, const Params& params) { - sgmax = 1.0 - params.gasOilParams().krnSwMdc(); - shmax = params.gasOilParams().krwSwMdc(); - somin = params.gasOilParams().pcSwMdc(); - Valgrind::CheckDefined(sgmax); - Valgrind::CheckDefined(shmax); - Valgrind::CheckDefined(somin); + if constexpr (isHysteresisDependent) { + sgmax = 1.0 - params.gasOilParams().krnSwMdc(); + shmax = params.gasOilParams().krwSwMdc(); + somin = params.gasOilParams().pcSwMdc(); + Valgrind::CheckDefined(sgmax); + Valgrind::CheckDefined(shmax); + Valgrind::CheckDefined(somin); + } } /* @@ -247,7 +259,9 @@ class EclTwoPhaseMaterial : public TraitsT const Scalar& somin, Params& params) { - params.gasOilParams().update(somin , shmax, 1.0 - sgmax); + if constexpr (isHysteresisDependent) { + params.gasOilParams().update(somin , shmax, 1.0 - sgmax); + } } static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping){ diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp index fec9c75ca0f..c72bfae315f 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp @@ -93,6 +93,8 @@ class PiecewiseLinearTwoPhaseMaterial : public TraitsT //! are dependent on the phase composition static constexpr bool isCompositionDependent = false; + static constexpr bool isHysteresisDependent = false; + /*! * \brief The capillary pressure-saturation curve. */ diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp index 9c26f01de26..7b00b09e7ef 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp @@ -214,6 +214,33 @@ class PiecewiseLinearTwoPhaseMaterialParams : public EnsureFinalized std::copy(values.begin(), values.end(), krnSamples_.begin()); } + template + void serializeOp([[maybe_unused]] Serializer& serializer) + { + // only serializes dynamic state - see update() and updateDynamic_() + // therefore: this is a no-op! + } + + Scalar SnTrapped([[maybe_unused]] bool maximumTrapping) const + { + return 0.0; + } + + Scalar SnStranded([[maybe_unused]] Scalar sg, [[maybe_unused]] Scalar krg) const + { + return 0.0; + } + + Scalar SwTrapped() const + { + return 0.0; + } + + bool update([[maybe_unused]] Scalar pcSw, [[maybe_unused]] Scalar krwSw, [[maybe_unused]] Scalar krnSw) + { + return false; + } + private: void swapOrderIfPossibleThrowOtherwise_(ValueVector& swValues, ValueVector& values) const {