From 254e1524982dda9aa131a7915c6cf003a27e8632 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Wed, 14 Aug 2024 11:41:47 +0200 Subject: [PATCH 1/7] compiling, weird numerical issue --- CMakeLists_files.cmake | 5 + .../EclMaterialLawManagerSimple.cpp | 525 ++++++++++++++++++ .../EclMaterialLawManagerSimple.hpp | 492 ++++++++++++++++ .../EclMaterialLawManagerSimpleHystParams.cpp | 310 +++++++++++ .../EclMaterialLawManagerSimpleInitParams.cpp | 342 ++++++++++++ ...ialLawManagerSimpleReadEffectiveParams.cpp | 459 +++++++++++++++ 6 files changed, 2133 insertions(+) create mode 100644 opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp create mode 100644 opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp create mode 100644 opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp create mode 100644 opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp create mode 100644 opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp 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/EclMaterialLawManagerSimple.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp new file mode 100644 index 00000000000..2bf8bb934f8 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp @@ -0,0 +1,525 @@ +// -*- 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"); + // 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(); + if (realParams.approach() == EclTwoPhaseApproach::GasOil) { + realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); + realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); + } + else if (realParams.approach() == EclTwoPhaseApproach::GasWater) { + realParams.gasWaterParams().drainageParams().setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]); + realParams.gasWaterParams().drainageParams().setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]); + } + else if (realParams.approach() == EclTwoPhaseApproach::OilWater) { + realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); + realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[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 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]; + 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 +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..c049d115330 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp @@ -0,0 +1,492 @@ +// -*- 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 +#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 = SatCurveMultiplexer; + using OilWaterEffectiveTwoPhaseLaw = SatCurveMultiplexer; + using GasWaterEffectiveTwoPhaseLaw = SatCurveMultiplexer; + + 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; + +public: + // the three-phase material law used by the simulation + using MaterialLaw = EclMultiplexerMaterial; + 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(); } + + bool enablePCHysteresis() const + { return (enableHysteresis() && hysteresisConfig_->pcHysteresisModel() >= 0); } + + bool enableWettingHysteresis() const + { return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 4); } + + bool enableNonWettingHysteresis() const + { return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 0); } + + 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::GasOil; + + 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..8bd225551fe --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp @@ -0,0 +1,310 @@ +// -*- 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 + +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 diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp new file mode 100644 index 00000000000..a77f374cb07 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp @@ -0,0 +1,342 @@ +// -*- 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, + unsigned elemIdx) +{ + const auto& epsInfo = this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx]; + + auto oilWaterParams = hystParams.getOilWaterParams(); + auto gasOilParams = hystParams.getGasOilParams(); + auto gasWaterParams = hystParams.getGasWaterParams(); + materialParams.setApproach(this->parent_.threePhaseApproach_); + switch (materialParams.approach()) { + case EclMultiplexerApproach::Stone1: { + auto& realParams = materialParams.template getRealParams(); + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setSwl(epsInfo.Swl); + + if (!this->parent_.stoneEtas_.empty()) { + realParams.setEta(this->parent_.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(this->parent_.twoPhaseApproach_); + realParams.finalize(); + break; + } + + case EclMultiplexerApproach::OnePhase: { + // Nothing to do, no parameters. + break; + } + } // end switch() +} + +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..ac09b0b2e58 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp @@ -0,0 +1,459 @@ +// -*- 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() ) { + 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().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); + } + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + 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); + } + + 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 +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; + } + + 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 +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(); + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + 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 ); + + effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); + auto& realParams = effParams.template getRealParams(); + + 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(); + + 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().template 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 (!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 From a5ea43844cb3256f8e586d346d5cd24614738d78 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Wed, 14 Aug 2024 15:27:24 +0200 Subject: [PATCH 2/7] remove 3phase multiplexer level --- .../EclMaterialLawManagerSimple.cpp | 101 ++---------------- .../EclMaterialLawManagerSimple.hpp | 4 +- .../EclMaterialLawManagerSimpleInitParams.cpp | 62 ++--------- .../EclTwoPhaseMaterial.hpp | 1 + 4 files changed, 21 insertions(+), 147 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp index 2bf8bb934f8..d354f0a4ab9 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp @@ -242,84 +242,20 @@ connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const 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]); + auto& realParams = mlp; + if (realParams.approach() == EclTwoPhaseApproach::GasOil) { 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]); -// } + else if (realParams.approach() == EclTwoPhaseApproach::GasWater) { + realParams.gasWaterParams().drainageParams().setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]); + realParams.gasWaterParams().drainageParams().setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]); } - break; - - case EclMultiplexerApproach::Default: { - auto& realParams = mlp.template getRealParams(); + else if (realParams.approach() == EclTwoPhaseApproach::OilWater) { 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(); - if (realParams.approach() == EclTwoPhaseApproach::GasOil) { - realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); - } - else if (realParams.approach() == EclTwoPhaseApproach::GasWater) { - realParams.gasWaterParams().drainageParams().setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]); - realParams.gasWaterParams().drainageParams().setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]); - } - else if (realParams.approach() == EclTwoPhaseApproach::OilWater) { - realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[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; } @@ -412,29 +348,8 @@ EclMaterialLawManagerSimple:: 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!"); - } + auto& realParams = materialParams; + return realParams.oilWaterParams().drainageParams().scaledPoints(); } template diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp index c049d115330..28d5df05960 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp @@ -33,6 +33,7 @@ #include #include +#include "EclTwoPhaseMaterial.hpp" #include #include @@ -111,7 +112,8 @@ class EclMaterialLawManagerSimple public: // the three-phase material law used by the simulation - using MaterialLaw = EclMultiplexerMaterial; + // using MaterialLaw = EclMultiplexerMaterial; + using MaterialLaw = EclTwoPhaseMaterial; using MaterialLawParams = typename MaterialLaw::Params; using DirectionalMaterialLawParamsPtr = std::unique_ptr>; diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp index a77f374cb07..08b4c76558c 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp @@ -199,64 +199,20 @@ void EclMaterialLawManagerSimple::InitParams:: initThreePhaseParams_(HystParams &hystParams, MaterialLawParams& materialParams, - unsigned satRegionIdx, - unsigned elemIdx) + [[maybe_unused]] unsigned satRegionIdx, + [[maybe_unused]] unsigned elemIdx) { - const auto& epsInfo = this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx]; + // const auto& epsInfo = this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx]; auto oilWaterParams = hystParams.getOilWaterParams(); auto gasOilParams = hystParams.getGasOilParams(); auto gasWaterParams = hystParams.getGasWaterParams(); - materialParams.setApproach(this->parent_.threePhaseApproach_); - switch (materialParams.approach()) { - case EclMultiplexerApproach::Stone1: { - auto& realParams = materialParams.template getRealParams(); - realParams.setGasOilParams(gasOilParams); - realParams.setOilWaterParams(oilWaterParams); - realParams.setSwl(epsInfo.Swl); - - if (!this->parent_.stoneEtas_.empty()) { - realParams.setEta(this->parent_.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(this->parent_.twoPhaseApproach_); - realParams.finalize(); - break; - } - - case EclMultiplexerApproach::OnePhase: { - // Nothing to do, no parameters. - break; - } - } // end switch() + auto& realParams = materialParams; + realParams.setGasOilParams(gasOilParams); + realParams.setOilWaterParams(oilWaterParams); + realParams.setGasWaterParams(gasWaterParams); + realParams.setApproach(this->parent_.twoPhaseApproach_); + realParams.finalize(); } template diff --git a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp index 1042fcf957e..d0fbfdcb71e 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp @@ -33,6 +33,7 @@ #include #include +#include namespace Opm { From 63f9b107c91db6e4b10e7452fdfa2ed117e7f757 Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Wed, 14 Aug 2024 15:35:18 +0200 Subject: [PATCH 3/7] remove LET/linearInterp multiplexing layer --- .../EclMaterialLawManagerSimple.hpp | 6 +- ...ialLawManagerSimpleReadEffectiveParams.cpp | 97 ++----------------- 2 files changed, 12 insertions(+), 91 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp index 28d5df05960..1054ffe8d55 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp @@ -86,9 +86,9 @@ class EclMaterialLawManagerSimple using GasWaterTraits = TwoPhaseMaterialTraits; // the two-phase material law which is defined on effective (unscaled) saturations - using GasOilEffectiveTwoPhaseLaw = SatCurveMultiplexer; - using OilWaterEffectiveTwoPhaseLaw = SatCurveMultiplexer; - using GasWaterEffectiveTwoPhaseLaw = SatCurveMultiplexer; + using GasOilEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial; + using OilWaterEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial; + using GasWaterEffectiveTwoPhaseLaw = PiecewiseLinearTwoPhaseMaterial; using GasOilEffectiveTwoPhaseParams = typename GasOilEffectiveTwoPhaseLaw::Params; using OilWaterEffectiveTwoPhaseParams = typename OilWaterEffectiveTwoPhaseLaw::Params; diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp index ac09b0b2e58..417d845a39f 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleReadEffectiveParams.cpp @@ -113,43 +113,7 @@ readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx) else if (!slgofTables.empty()) readGasOilSlgof_(effParams, Swco, tolcrit, slgofTables.template 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(); + throw std::runtime_error("LET tables not supported!"); } break; } @@ -197,8 +161,7 @@ readGasOilFamily2_(GasOilEffectiveTwoPhaseParams& effParams, SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx); } - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sofTable.getColumn(columnName))); realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG"))); @@ -220,8 +183,7 @@ readGasOilSgof_(GasOilEffectiveTwoPhaseParams& effParams, SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx); } - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; realParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KROG"))); realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KRG"))); @@ -243,8 +205,7 @@ readGasOilSlgof_(GasOilEffectiveTwoPhaseParams& effParams, SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco; } - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; realParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KROG"))); realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KRG"))); @@ -279,8 +240,7 @@ readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionId case SatFuncControls::KeywordFamily::Family_II: { const TableContainer& sgwfnTables = tableManager.getSgwfnTables(); - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; if (!sgwfnTables.empty()){ const SgwfnTable& sgwfnTable = tableManager.getSgwfnTables().template getTable( satRegionIdx ); std::vector SwSamples(sgwfnTable.numRows()); @@ -314,8 +274,7 @@ readGasWaterParameters_(GasWaterEffectiveParamVector& dest, unsigned satRegionId const GsfTable& gsfTable = tableManager.getGsfTables().template getTable( satRegionIdx ); const WsfTable& wsfTable = tableManager.getWsfTables().template getTable( satRegionIdx ); - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; std::vector SwColumn = wsfTable.getColumn("SW").vectorCopy(); @@ -359,8 +318,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId const auto& swofTable = tableManager.getSwofTables().template getTable(satRegionIdx); const std::vector SwColumn = swofTable.getColumn("SW").vectorCopy(); - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KRW"))); realParams.setKrnSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KROW"))); @@ -368,43 +326,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId 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(); + throw std::runtime_error("LET tables not supported!"); } break; } @@ -414,8 +336,7 @@ readOilWaterParameters_(OilWaterEffectiveParamVector& dest, unsigned satRegionId const auto& swfnTable = tableManager.getSwfnTables().template getTable(satRegionIdx); const std::vector SwColumn = swfnTable.getColumn("SW").vectorCopy(); - effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear); - auto& realParams = effParams.template getRealParams(); + auto& realParams = effParams; realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW"))); realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy()); From b97ec8bac7f9700a365cfe133f0c254543ca5f13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 28 Aug 2024 11:20:42 +0200 Subject: [PATCH 4/7] Add static isHysteresisDependent flag. --- .../EclHysteresisTwoPhaseLaw.hpp | 2 + .../EclTwoPhaseMaterial.hpp | 41 ++++++++++++------- .../PiecewiseLinearTwoPhaseMaterial.hpp | 2 + 3 files changed, 31 insertions(+), 14 deletions(-) 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/EclTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp index d0fbfdcb71e..29d771b9892 100644 --- a/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp @@ -112,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*/) { @@ -193,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); + } } /* @@ -213,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); + } } @@ -229,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); + } } /* @@ -248,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. */ From 78df2b27efef9f0c16c6adbbdc158c50a8b47cc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 28 Aug 2024 11:21:44 +0200 Subject: [PATCH 5/7] [HACK] add some dummy methods. --- .../PiecewiseLinearTwoPhaseMaterialParams.hpp | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) 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 { From 049d5d150d9c492d67f71623a5b58b59ab8f1d97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 28 Aug 2024 11:22:41 +0200 Subject: [PATCH 6/7] Remove hysteresis and scaling layers from simple law. --- .../EclMaterialLawManagerSimple.cpp | 111 ++++++++-------- .../EclMaterialLawManagerSimple.hpp | 85 ++++++------ .../EclMaterialLawManagerSimpleHystParams.cpp | 4 + .../EclMaterialLawManagerSimpleInitParams.cpp | 122 +++++++++--------- 4 files changed, 166 insertions(+), 156 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp index d354f0a4ab9..8b9a956db15 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.cpp @@ -54,17 +54,17 @@ initFromState(const EclipseState& eclState) this->hasOil = ph.active(Phase::OIL); this->hasWater = ph.active(Phase::WATER); - readGlobalEpsOptions_(eclState); - readGlobalHysteresisOptions_(eclState); + // 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); + // 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(); @@ -118,15 +118,15 @@ initFromState(const EclipseState& eclState) } // 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]); - } - } + // 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 @@ -207,10 +207,10 @@ applySwatinit(unsigned elemIdx, else elemScaledEpsInfo.maxPcow = newMaxPcow; - auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx); - elemEclEpsScalingPoints.init(elemScaledEpsInfo, - *oilWaterEclEpsConfig_, - EclTwoPhaseSystemType::OilWater); + // auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx); + // elemEclEpsScalingPoints.init(elemScaledEpsInfo, + // *oilWaterEclEpsConfig_, + // EclTwoPhaseSystemType::OilWater); return {Sw, newSwatInit}; } @@ -227,10 +227,10 @@ EclMaterialLawManagerSimple::applyRestartSwatInit(const unsigned elemId elemScaledEpsInfo.maxPcow = maxPcow; - this->oilWaterScaledEpsPointsDrainage(elemIdx) - .init(elemScaledEpsInfo, - *this->oilWaterEclEpsConfig_, - EclTwoPhaseSystemType::OilWater); + // this->oilWaterScaledEpsPointsDrainage(elemIdx) + // .init(elemScaledEpsInfo, + // *this->oilWaterEclEpsConfig_, + // EclTwoPhaseSystemType::OilWater); } template @@ -245,16 +245,13 @@ connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const auto& realParams = mlp; if (realParams.approach() == EclTwoPhaseApproach::GasOil) { - realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]); - realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]); + realParams.setGasOilParams(gasOilEffectiveParamVector_[satRegionIdx]); } else if (realParams.approach() == EclTwoPhaseApproach::GasWater) { - realParams.gasWaterParams().drainageParams().setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]); - realParams.gasWaterParams().drainageParams().setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]); + realParams.setGasWaterParams(gasWaterEffectiveParamVector_[satRegionIdx]); } else if (realParams.approach() == EclTwoPhaseApproach::OilWater) { - realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]); - realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]); + realParams.setOilWaterParams(oilWaterEffectiveParamVector_[satRegionIdx]); } return mlp; } @@ -342,15 +339,15 @@ setGasOilHysteresisParams(const Scalar& sgmax, 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 +// EclEpsScalingPoints& +// EclMaterialLawManagerSimple:: +// oilWaterScaledEpsPointsDrainage(unsigned elemIdx) +// { +// auto& materialParams = materialLawParams_[elemIdx]; +// auto& realParams = materialParams; +// return realParams.oilWaterParams().drainageParams().scaledPoints(); +// } template const typename EclMaterialLawManagerSimple::MaterialLawParams& EclMaterialLawManagerSimple:: @@ -377,23 +374,23 @@ materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const } } -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:: +// 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:: diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp index 1054ffe8d55..a6ed36ddc4c 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimple.hpp @@ -95,20 +95,24 @@ class EclMaterialLawManagerSimple 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; + // 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 = 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 @@ -126,13 +130,13 @@ class EclMaterialLawManagerSimple using OilWaterEffectiveParamVector = std::vector>; using GasWaterEffectiveParamVector = std::vector>; - using GasOilScalingPointsVector = std::vector>>; - using OilWaterScalingPointsVector = std::vector>>; - using GasWaterScalingPointsVector = 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 GasOilParamVector = std::vector>; + // using OilWaterParamVector = std::vector>; + // using GasWaterParamVector = std::vector>; using MaterialLawParamsVector = std::vector>; // helper classes @@ -148,7 +152,7 @@ class EclMaterialLawManagerSimple void run(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, const std::function& lookupIdxOnLevelZeroAssigner); private: - class HystParams; + // 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)>& @@ -170,17 +174,17 @@ class EclMaterialLawManagerSimple void initSatnumRegionArray_(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner); void initThreePhaseParams_( - HystParams &hystParams, + // HystParams &hystParams, MaterialLawParams& materialParams, unsigned satRegionIdx, unsigned elemIdx); void readEffectiveParameters_(); - void readUnscaledEpsPointsVectors_(); - template - void readUnscaledEpsPoints_(Container& dest, std::shared_ptr config, EclTwoPhaseSystemType system_type); + // 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: @@ -228,6 +232,7 @@ class EclMaterialLawManagerSimple std::shared_ptr oilWaterParams_; std::shared_ptr gasWaterParams_; }; +*/ // This class' implementation is defined in "EclMaterialLawManagerSimpleReadEffectiveParams.cpp" class ReadEffectiveParams { @@ -313,16 +318,20 @@ class EclMaterialLawManagerSimple { return enablePpcwmax_; } bool enableHysteresis() const - { return hysteresisConfig_->enableHysteresis(); } + //{ return hysteresisConfig_->enableHysteresis(); } + { return false; } bool enablePCHysteresis() const - { return (enableHysteresis() && hysteresisConfig_->pcHysteresisModel() >= 0); } + //{ return (enableHysteresis() && hysteresisConfig_->pcHysteresisModel() >= 0); } + { return false; } bool enableWettingHysteresis() const - { return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 4); } + //{ return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 4); } + { return false; } bool enableNonWettingHysteresis() const - { return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 0); } + //{ return (enableHysteresis() && hysteresisConfig_->krHysteresisModel() >= 0); } + { return false; } MaterialLawParams& materialLawParams(unsigned elemIdx) { @@ -415,7 +424,7 @@ class EclMaterialLawManagerSimple const Scalar& somin, unsigned elemIdx); - EclEpsScalingPoints& oilWaterScaledEpsPointsDrainage(unsigned elemIdx); + // EclEpsScalingPoints& oilWaterScaledEpsPointsDrainage(unsigned elemIdx); const EclEpsScalingPointsInfo& oilWaterScaledEpsInfoDrainage(size_t elemIdx) const { return oilWaterScaledEpsInfoDrainage_[elemIdx]; } @@ -435,26 +444,26 @@ class EclMaterialLawManagerSimple private: const MaterialLawParams& materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const; - void readGlobalEpsOptions_(const EclipseState& eclState); + // void readGlobalEpsOptions_(const EclipseState& eclState); - void readGlobalHysteresisOptions_(const EclipseState& state); + // void readGlobalHysteresisOptions_(const EclipseState& state); void readGlobalThreePhaseOptions_(const Runspec& runspec); bool enableEndPointScaling_; - std::shared_ptr hysteresisConfig_; - std::vector> wagHystersisConfig_; + // std::shared_ptr hysteresisConfig_; + // std::vector> wagHystersisConfig_; std::shared_ptr oilWaterEclEpsConfig_; std::vector> unscaledEpsInfo_; OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_; - std::shared_ptr gasWaterEclEpsConfig_; + // std::shared_ptr gasWaterEclEpsConfig_; - GasOilScalingPointsVector gasOilUnscaledPointsVector_; - OilWaterScalingPointsVector oilWaterUnscaledPointsVector_; - GasWaterScalingPointsVector gasWaterUnscaledPointsVector_; + // GasOilScalingPointsVector gasOilUnscaledPointsVector_; + // OilWaterScalingPointsVector oilWaterUnscaledPointsVector_; + // GasWaterScalingPointsVector gasWaterUnscaledPointsVector_; GasOilEffectiveParamVector gasOilEffectiveParamVector_; OilWaterEffectiveParamVector oilWaterEffectiveParamVector_; @@ -462,7 +471,7 @@ class EclMaterialLawManagerSimple EclMultiplexerApproach threePhaseApproach_ = EclMultiplexerApproach::Default; // this attribute only makes sense for twophase simulations! - enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::GasOil; + enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::GasWater; std::vector materialLawParams_; DirectionalMaterialLawParamsPtr dirMaterialLawParams_; diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp index 8bd225551fe..b9df445069a 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleHystParams.cpp @@ -18,6 +18,8 @@ along with OPM. If not, see . */ +#if 0 + #include #include #include @@ -308,3 +310,5 @@ template class EclMaterialLawManagerSimple } // namespace Opm + +#endif // 0 diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp index 08b4c76558c..1dd907a781b 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp @@ -55,7 +55,7 @@ EclMaterialLawManagerSimple::InitParams:: run(const std::function(const FieldPropsManager&, const std::string&, bool)>& fieldPropIntOnLeafAssigner, const std::function& lookupIdxOnLevelZeroAssigner) { - readUnscaledEpsPointsVectors_(); + // readUnscaledEpsPointsVectors_(); readEffectiveParameters_(); initSatnumRegionArray_(fieldPropIntOnLeafAssigner); copySatnumArrays_(fieldPropIntOnLeafAssigner); @@ -70,19 +70,19 @@ run(const std::function(const FieldPropsManager&, const std::st for (unsigned elemIdx = 0; elemIdx < this->numCompressedElems_; ++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); + // 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); } } } @@ -197,20 +197,20 @@ initSatnumRegionArray_(const std::function(const FieldPropsMana template void EclMaterialLawManagerSimple::InitParams:: -initThreePhaseParams_(HystParams &hystParams, +initThreePhaseParams_(//HystParams &hystParams, MaterialLawParams& materialParams, [[maybe_unused]] unsigned satRegionIdx, [[maybe_unused]] unsigned elemIdx) { // const auto& epsInfo = this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx]; - auto oilWaterParams = hystParams.getOilWaterParams(); - auto gasOilParams = hystParams.getGasOilParams(); - auto gasWaterParams = hystParams.getGasWaterParams(); + // auto oilWaterParams = hystParams.getOilWaterParams(); + // auto gasOilParams = hystParams.getGasOilParams(); + // auto gasWaterParams = hystParams.getGasWaterParams(); auto& realParams = materialParams; - realParams.setGasOilParams(gasOilParams); - realParams.setOilWaterParams(oilWaterParams); - realParams.setGasWaterParams(gasWaterParams); + // realParams.setGasOilParams(gasOilParams); + // realParams.setOilWaterParams(oilWaterParams); + // realParams.setGasWaterParams(gasWaterParams); realParams.setApproach(this->parent_.twoPhaseApproach_); realParams.finalize(); } @@ -225,47 +225,47 @@ readEffectiveParameters_() 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 +// 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 +// 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 From bb2cb135841b78356fcf5841c28ae77153be14c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 28 Aug 2024 11:26:19 +0200 Subject: [PATCH 7/7] Set the actual parameter pointers. --- .../EclMaterialLawManagerSimpleInitParams.cpp | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp index 1dd907a781b..0ce62ed1230 100644 --- a/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp +++ b/opm/material/fluidmatrixinteractions/EclMaterialLawManagerSimpleInitParams.cpp @@ -199,19 +199,21 @@ void EclMaterialLawManagerSimple::InitParams:: initThreePhaseParams_(//HystParams &hystParams, MaterialLawParams& materialParams, - [[maybe_unused]] unsigned satRegionIdx, + unsigned satRegionIdx, [[maybe_unused]] unsigned elemIdx) { - // const auto& epsInfo = this->parent_.oilWaterScaledEpsInfoDrainage_[elemIdx]; - - // auto oilWaterParams = hystParams.getOilWaterParams(); - // auto gasOilParams = hystParams.getGasOilParams(); - // auto gasWaterParams = hystParams.getGasWaterParams(); auto& realParams = materialParams; - // realParams.setGasOilParams(gasOilParams); - // realParams.setOilWaterParams(oilWaterParams); - // realParams.setGasWaterParams(gasWaterParams); 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(); }