OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #3364 from hakonhagland/imbnum3

Browse Source 
Adds support for directional relative permeabilities with hysteresis
This commit is contained in:
Tor Harald Sandve
2023-01-30 12:33:36 +01:00
committed by GitHub
parent b83f10d2d7 5430dc46a4
commit 9b3e1f92bf
9 changed files with 1295 additions and 834 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
CMakeLists_files.cmake
View File

@@ -259,6 +259,9 @@ if(ENABLE_ECL_INPUT)
src/opm/material/fluidmatrixinteractions/EclEpsGridProperties.cpp
src/opm/material/fluidmatrixinteractions/EclHysteresisConfig.cpp
src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp
src/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.cpp
src/opm/material/fluidsystems/blackoilpvt/Co2GasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.cpp
@@ -859,6 +862,8 @@ list( APPEND PUBLIC_HEADER_FILES
opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp
opm/material/fluidmatrixinteractions/TwoPhaseLETCurves.hpp
opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp
opm/material/fluidmatrixinteractions/DirectionalMaterialLawParams.hpp
opm/material/fluidmatrixinteractions/DirectionalMaterialLawParams.hpp
opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp
opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp
opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp

3
opm/input/eclipse/EclipseState/Grid/FieldProps.hpp
View File

@@ -247,6 +247,9 @@ static const std::unordered_map<std::string, keyword_info<int>> int_keywords = {
{"KRNUMX", keyword_info<int>{}},
{"KRNUMY", keyword_info<int>{}},
{"KRNUMZ", keyword_info<int>{}},
{"IMBNUMX", keyword_info<int>{}},
{"IMBNUMY", keyword_info<int>{}},
{"IMBNUMZ", keyword_info<int>{}},
};
}

54
opm/material/fluidmatrixinteractions/DirectionalMaterialLawParams.hpp Normal file
View File

@@ -0,0 +1,54 @@
/*
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
(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 <http://www.gnu.org/licenses/>.
*/
/*!
* \file
*
* \brief This file contains definitions related to directional material law parameters
*/
#ifndef OPM_DIRECTIONAL_MATERIAL_LAW_PARAMS_HH
#define OPM_DIRECTIONAL_MATERIAL_LAW_PARAMS_HH
#include <cstddef>
#include <vector>
namespace Opm {
template <class MaterialLawParams>
struct DirectionalMaterialLawParams {
using vector_type = std::vector<MaterialLawParams>;
DirectionalMaterialLawParams() : materialLawParamsX_{}, materialLawParamsY_{}, materialLawParamsZ_{} {}
DirectionalMaterialLawParams(std::size_t size) :
materialLawParamsX_(size), materialLawParamsY_(size), materialLawParamsZ_(size) {}
vector_type& getArray(int index) {
switch(index) {
case 0:
return materialLawParamsX_;
case 1:
return materialLawParamsY_;
case 2:
return materialLawParamsZ_;
default:
throw std::runtime_error("Unexpected mobility array index");
}
}
vector_type materialLawParamsX_;
vector_type materialLawParamsY_;
vector_type materialLawParamsZ_;
};
} // namespace Opm
#endif

225
opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp
View File

@@ -38,6 +38,9 @@
#include <opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp>
#include <opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
#include <opm/material/fluidmatrixinteractions/DirectionalMaterialLawParams.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableColumn.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <cassert>
#include <memory>
@@ -56,8 +59,6 @@ class Runspec;
class SgfnTable;
class SgofTable;
class SlgofTable;
class Sof2Table;
class Sof3Table;
/*!
* \ingroup fluidmatrixinteractions
@@ -109,6 +110,7 @@ public:
// the three-phase material law used by the simulation
using MaterialLaw = EclMultiplexerMaterial<Traits, GasOilTwoPhaseLaw, OilWaterTwoPhaseLaw, GasWaterTwoPhaseLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
using DirectionalMaterialLawParamsPtr = std::unique_ptr<DirectionalMaterialLawParams<MaterialLawParams>>;
EclMaterialLawManager();
~EclMaterialLawManager();
@@ -128,6 +130,113 @@ private:
using GasWaterParamVector = std::vector<std::shared_ptr<GasWaterTwoPhaseHystParams>>;
using MaterialLawParamsVector = std::vector<std::shared_ptr<MaterialLawParams>>;
// helper classes
// This class' implementation is defined in "EclMaterialLawManagerInitParams.cpp"
class InitParams {
public:
InitParams(EclMaterialLawManager<TraitsT>& parent, const EclipseState& eclState, size_t numCompressedElems);
void run();
private:
class HystParams;
void copySatnumArrays_();
void copyIntArray_(std::vector<int>& dest, const std::string keyword);
unsigned imbRegion_(std::vector<int>& array, unsigned elemIdx);
void initArrays_(
std::vector<std::vector<int>*>& satnumArray,
std::vector<std::vector<int>*>& imbnumArray,
std::vector<std::vector<MaterialLawParams>*>& mlpArray);
void initMaterialLawParamVectors_();
void initOilWaterScaledEpsInfo_();
void initSatnumRegionArray_();
void initThreePhaseParams_(
HystParams &hystParams,
MaterialLawParams& materialParams,
unsigned satRegionIdx,
unsigned elemIdx);
void readEffectiveParameters_();
void readUnscaledEpsPointsVectors_();
template <class Container>
void readUnscaledEpsPoints_(Container& dest, std::shared_ptr<EclEpsConfig> config, EclTwoPhaseSystemType system_type);
unsigned satRegion_(std::vector<int>& array, unsigned elemIdx);
unsigned satOrImbRegion_(std::vector<int>& array, std::vector<int>& default_vec, unsigned elemIdx);
// This class' implementation is defined in "EclMaterialLawManagerHystParams.cpp"
class HystParams {
public:
HystParams(EclMaterialLawManager<TraitsT>::InitParams& init_params);
void finalize();
std::shared_ptr<GasOilTwoPhaseHystParams> getGasOilParams();
std::shared_ptr<OilWaterTwoPhaseHystParams> getOilWaterParams();
std::shared_ptr<GasWaterTwoPhaseHystParams> getGasWaterParams();
void setConfig();
void setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx);
void setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx);
void setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx);
void setImbibitionParamsOilGas(unsigned elemIdx, unsigned satRegionIdx);
void setImbibitionParamsOilWater(unsigned elemIdx, unsigned satRegionIdx);
void setImbibitionParamsGasWater(unsigned elemIdx, unsigned satRegionIdx);
private:
bool hasGasWater_();
bool hasGasOil_();
bool hasOilWater_();
std::tuple<EclEpsScalingPointsInfo<Scalar>, EclEpsScalingPoints<Scalar>> readScaledEpsPoints_(
const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, EclTwoPhaseSystemType type);
std::tuple<EclEpsScalingPointsInfo<Scalar>, EclEpsScalingPoints<Scalar>>
readScaledEpsPointsDrainage_(unsigned elemIdx, EclTwoPhaseSystemType type);
std::tuple<EclEpsScalingPointsInfo<Scalar>, EclEpsScalingPoints<Scalar>>
readScaledEpsPointsImbibition_(unsigned elemIdx, EclTwoPhaseSystemType type);
EclMaterialLawManager<TraitsT>::InitParams& init_params_;
EclMaterialLawManager<TraitsT>& parent_;
const EclipseState& eclState_;
std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams_;
std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams_;
std::shared_ptr<GasWaterTwoPhaseHystParams> gasWaterParams_;
};
// This class' implementation is defined in "EclMaterialLawManagerReadEffectiveParams.cpp"
class ReadEffectiveParams {
public:
ReadEffectiveParams(EclMaterialLawManager<TraitsT>::InitParams& init_params);
void read();
private:
std::vector<double> normalizeKrValues_(const double tolcrit, const TableColumn& krValues) const;
void readGasOilParameters_(GasOilEffectiveParamVector& dest, unsigned satRegionIdx);
template <class TableType>
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);
EclMaterialLawManager<TraitsT>::InitParams& init_params_;
EclMaterialLawManager<TraitsT>& parent_;
const EclipseState& eclState_;
}; // end of "class ReadEffectiveParams"
EclMaterialLawManager<TraitsT>& parent_;
const EclipseState& eclState_;
size_t numCompressedElems_;
std::unique_ptr<EclEpsGridProperties> epsImbGridProperties_; //imbibition
std::unique_ptr<EclEpsGridProperties> epsGridProperties_; // drainage
}; // end of "class InitParams"
public:
void initFromState(const EclipseState& eclState);
@@ -163,6 +272,16 @@ public:
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<MaterialLawParams&>(materialLawParamsFunc_(elemIdx, facedir));
}
/*!
* \brief Returns a material parameter object for a given element and saturation region.
*
@@ -183,6 +302,13 @@ public:
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]; }
@@ -192,7 +318,15 @@ public:
if (!enableHysteresis())
return;
MaterialLaw::updateHysteresis(materialLawParams_[elemIdx], fluidState);
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<ndim; i++) {
MaterialLaw::updateHysteresis(materialLawParams(elemIdx, facedirs[i]), fluidState);
}
}
}
void oilWaterHysteresisParams(Scalar& pcSwMdc,
@@ -217,83 +351,14 @@ public:
{ return oilWaterScaledEpsInfoDrainage_[elemIdx]; }
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);
template <class Container>
void readGasOilEffectiveParameters_(Container& dest,
const EclipseState& eclState,
unsigned satRegionIdx);
void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const SgofTable& sgofTable);
void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const SlgofTable& slgofTable);
void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const Sof3Table& sof3Table,
const SgfnTable& sgfnTable);
void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const Sof2Table& sof2Table,
const SgfnTable& sgfnTable);
template <class Container>
void readOilWaterEffectiveParameters_(Container& dest,
const EclipseState& eclState,
unsigned satRegionIdx);
template <class Container>
void readGasWaterEffectiveParameters_(Container& dest,
const EclipseState& eclState,
unsigned satRegionIdx);
template <class Container>
void readGasOilUnscaledPoints_(Container& dest,
std::shared_ptr<EclEpsConfig> config,
const EclipseState& /* eclState */,
unsigned satRegionIdx);
template <class Container>
void readOilWaterUnscaledPoints_(Container& dest,
std::shared_ptr<EclEpsConfig> config,
const EclipseState& /* eclState */,
unsigned satRegionIdx);
template <class Container>
void readGasWaterUnscaledPoints_(Container& dest,
std::shared_ptr<EclEpsConfig> config,
const EclipseState& /* eclState */,
unsigned satRegionIdx);
std::tuple<EclEpsScalingPointsInfo<Scalar>,
EclEpsScalingPoints<Scalar>>
readScaledPoints_(const EclEpsConfig& config,
const EclipseState& eclState,
const EclEpsGridProperties& epsGridProperties,
unsigned elemIdx,
EclTwoPhaseSystemType type);
void initThreePhaseParams_(const EclipseState& /* eclState */,
MaterialLawParams& materialParams,
unsigned satRegionIdx,
const EclEpsScalingPointsInfo<Scalar>& epsInfo,
std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams,
std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams,
std::shared_ptr<GasWaterTwoPhaseHystParams> gasWaterParams);
bool enableEndPointScaling_;
std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;
@@ -316,21 +381,25 @@ private:
enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::GasOil;
std::vector<MaterialLawParams> materialLawParams_;
DirectionalMaterialLawParamsPtr dirMaterialLawParams_;
std::vector<int> satnumRegionArray_;
std::vector<int> krnumXArray_;
std::vector<int> krnumYArray_;
std::vector<int> krnumZArray_;
std::vector<int> imbnumXArray_;
std::vector<int> imbnumYArray_;
std::vector<int> imbnumZArray_;
std::vector<int> imbnumRegionArray_;
std::vector<Scalar> stoneEtas;
std::vector<Scalar> stoneEtas_;
bool hasGas;
bool hasOil;
bool hasWater;
std::shared_ptr<EclEpsConfig> gasOilConfig;
std::shared_ptr<EclEpsConfig> oilWaterConfig;
std::shared_ptr<EclEpsConfig> gasWaterConfig;
std::shared_ptr<EclEpsConfig> gasOilConfig_;
std::shared_ptr<EclEpsConfig> oilWaterConfig_;
std::shared_ptr<EclEpsConfig> gasWaterConfig_;
};
} // namespace Opm

6
src/opm/input/eclipse/share/keywords/000_Eclipse100/I/IMBNUM
View File

@@ -3,6 +3,12 @@
"sections": [
"REGIONS"
],
"deck_names": [
"IMBNUM",
"IMBNUMX",
"IMBNUMY",
"IMBNUMZ"
],
"data": {
"value_type": "INT"
}

789
src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp
View File

@@ -26,41 +26,13 @@
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SgfnTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SgofTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SlgofTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Sof2Table.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Sof3Table.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SwfnTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SwofTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableColumn.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
#include <algorithm>
namespace {
// Relative permeability values not strictly greater than 'tolcrit' treated as zero.
std::vector<double> normalizeKrValues_(const double tolcrit,
const Opm::TableColumn& krValues)
{
auto kr = krValues.vectorCopy();
std::transform(kr.begin(), kr.end(), kr.begin(),
[tolcrit](const double kri)
{
return (kri > tolcrit) ? kri : 0.0;
});
return kr;
}
}
namespace Opm {
template<class TraitsT>
@@ -87,12 +59,12 @@ initFromState(const EclipseState& eclState)
readGlobalThreePhaseOptions_(runspec);
// Read the end point scaling configuration (once per run).
gasOilConfig = std::make_shared<EclEpsConfig>();
oilWaterConfig = std::make_shared<EclEpsConfig>();
gasWaterConfig = std::make_shared<EclEpsConfig>();
gasOilConfig->initFromState(eclState, EclTwoPhaseSystemType::GasOil);
oilWaterConfig->initFromState(eclState, EclTwoPhaseSystemType::OilWater);
gasWaterConfig->initFromState(eclState, EclTwoPhaseSystemType::GasWater);
gasOilConfig_ = std::make_shared<EclEpsConfig>();
oilWaterConfig_ = std::make_shared<EclEpsConfig>();
gasWaterConfig_ = std::make_shared<EclEpsConfig>();
gasOilConfig_->initFromState(eclState, EclTwoPhaseSystemType::GasOil);
oilWaterConfig_->initFromState(eclState, EclTwoPhaseSystemType::OilWater);
gasWaterConfig_->initFromState(eclState, EclTwoPhaseSystemType::GasWater);
const auto& tables = eclState.getTableManager();
@@ -101,11 +73,11 @@ initFromState(const EclipseState& eclState)
const auto& stone1exTables = tables.getStone1exTable();
if (! stone1exTables.empty()) {
stoneEtas.clear();
stoneEtas.reserve(numSatRegions);
stoneEtas_.clear();
stoneEtas_.reserve(numSatRegions);
for (const auto& table : stone1exTables) {
stoneEtas.push_back(table.eta);
stoneEtas_.push_back(table.eta);
}
}
}
@@ -134,232 +106,8 @@ template<class TraitsT>
void EclMaterialLawManager<TraitsT>::
initParamsForElements(const EclipseState& eclState, size_t numCompressedElems)
{
// get the number of saturation regions
const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
// setup the saturation region specific parameters
gasOilUnscaledPointsVector_.resize(numSatRegions);
oilWaterUnscaledPointsVector_.resize(numSatRegions);
gasWaterUnscaledPointsVector_.resize(numSatRegions);
gasOilEffectiveParamVector_.resize(numSatRegions);
oilWaterEffectiveParamVector_.resize(numSatRegions);
gasWaterEffectiveParamVector_.resize(numSatRegions);
for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
// unscaled points for end-point scaling
readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, eclState, satRegionIdx);
readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, eclState, satRegionIdx);
readGasWaterUnscaledPoints_(gasWaterUnscaledPointsVector_, gasWaterConfig, eclState, satRegionIdx);
// the parameters for the effective two-phase matererial laws
readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, eclState, satRegionIdx);
readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, eclState, satRegionIdx);
readGasWaterEffectiveParameters_(gasWaterEffectiveParamVector_, eclState, satRegionIdx);
}
// copy the SATNUM grid property. in some cases this is not necessary, but it
// should not require much memory anyway...
satnumRegionArray_.resize(numCompressedElems);
if (eclState.fieldProps().has_int("SATNUM")) {
const auto& satnumRawData = eclState.fieldProps().get_int("SATNUM");
for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
satnumRegionArray_[elemIdx] = satnumRawData[elemIdx] - 1;
}
}
else {
std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0);
}
auto copy_krnum = [&eclState, numCompressedElems](std::vector<int>& dest, const std::string keyword) {
if (eclState.fieldProps().has_int(keyword)) {
dest.resize(numCompressedElems);
const auto& satnumRawData = eclState.fieldProps().get_int(keyword);
for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
dest[elemIdx] = satnumRawData[elemIdx] - 1;
}
}
};
copy_krnum(krnumXArray_, "KRNUMX");
copy_krnum(krnumYArray_, "KRNUMY");
copy_krnum(krnumZArray_, "KRNUMZ");
// create the information for the imbibition region (IMBNUM). By default this is
// the same as the saturation region (SATNUM)
imbnumRegionArray_ = satnumRegionArray_;
if (eclState.fieldProps().has_int("IMBNUM")) {
const auto& imbnumRawData = eclState.fieldProps().get_int("IMBNUM");
for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
imbnumRegionArray_[elemIdx] = imbnumRawData[elemIdx] - 1;
}
}
assert(numCompressedElems == satnumRegionArray_.size());
assert(!enableHysteresis() || numCompressedElems == imbnumRegionArray_.size());
// read the scaled end point scaling parameters which are specific for each
// element
oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems);
std::unique_ptr<EclEpsGridProperties> epsImbGridProperties;
if (enableHysteresis()) {
epsImbGridProperties = std::make_unique<EclEpsGridProperties>(eclState, true);
}
EclEpsGridProperties epsGridProperties(eclState, false);
materialLawParams_.resize(numCompressedElems);
for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
unsigned satRegionIdx = static_cast<unsigned>(satnumRegionArray_[elemIdx]);
auto gasOilParams = std::make_shared<GasOilTwoPhaseHystParams>();
auto oilWaterParams = std::make_shared<OilWaterTwoPhaseHystParams>();
auto gasWaterParams = std::make_shared<GasWaterTwoPhaseHystParams>();
gasOilParams->setConfig(hysteresisConfig_);
oilWaterParams->setConfig(hysteresisConfig_);
gasWaterParams->setConfig(hysteresisConfig_);
auto [gasOilScaledInfo, gasOilScaledPoint] =
readScaledPoints_(*gasOilConfig,
eclState,
epsGridProperties,
elemIdx,
EclTwoPhaseSystemType::GasOil);
auto [owinfo, oilWaterScaledEpsPointDrainage] =
readScaledPoints_(*oilWaterConfig,
eclState,
epsGridProperties,
elemIdx,
EclTwoPhaseSystemType::OilWater);
oilWaterScaledEpsInfoDrainage_[elemIdx] = owinfo;
auto [gasWaterScaledInfo, gasWaterScaledPoint] =
readScaledPoints_(*gasWaterConfig,
eclState,
epsGridProperties,
elemIdx,
EclTwoPhaseSystemType::GasWater);
if (hasGas && hasOil) {
GasOilEpsTwoPhaseParams gasOilDrainParams;
gasOilDrainParams.setConfig(gasOilConfig);
gasOilDrainParams.setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
gasOilDrainParams.setScaledPoints(gasOilScaledPoint);
gasOilDrainParams.setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
gasOilDrainParams.finalize();
gasOilParams->setDrainageParams(gasOilDrainParams,
gasOilScaledInfo,
EclTwoPhaseSystemType::GasOil);
}
if (hasOil && hasWater) {
OilWaterEpsTwoPhaseParams oilWaterDrainParams;
oilWaterDrainParams.setConfig(oilWaterConfig);
oilWaterDrainParams.setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
oilWaterDrainParams.setScaledPoints(oilWaterScaledEpsPointDrainage);
oilWaterDrainParams.setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
oilWaterDrainParams.finalize();
oilWaterParams->setDrainageParams(oilWaterDrainParams,
owinfo,
EclTwoPhaseSystemType::OilWater);
}
if (hasGas && hasWater && !hasOil) {
GasWaterEpsTwoPhaseParams gasWaterDrainParams;
gasWaterDrainParams.setConfig(gasWaterConfig);
gasWaterDrainParams.setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]);
gasWaterDrainParams.setScaledPoints(gasWaterScaledPoint);
gasWaterDrainParams.setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]);
gasWaterDrainParams.finalize();
gasWaterParams->setDrainageParams(gasWaterDrainParams,
gasWaterScaledInfo,
EclTwoPhaseSystemType::GasWater);
}
if (enableHysteresis()) {
auto [gasOilScaledImbInfo, gasOilScaledImbPoint] =
readScaledPoints_(*gasOilConfig,
eclState,
*epsImbGridProperties,
elemIdx,
EclTwoPhaseSystemType::GasOil);
auto [oilWaterScaledImbInfo, oilWaterScaledImbPoint] =
readScaledPoints_(*oilWaterConfig,
eclState,
*epsImbGridProperties,
elemIdx,
EclTwoPhaseSystemType::OilWater);
auto [gasWaterScaledImbInfo, gasWaterScaledImbPoint] =
readScaledPoints_(*gasWaterConfig,
eclState,
*epsImbGridProperties,
elemIdx,
EclTwoPhaseSystemType::GasWater);
unsigned imbRegionIdx = imbnumRegionArray_[elemIdx];
if (hasGas && hasOil) {
GasOilEpsTwoPhaseParams gasOilImbParamsHyst;
gasOilImbParamsHyst.setConfig(gasOilConfig);
gasOilImbParamsHyst.setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]);
gasOilImbParamsHyst.setScaledPoints(gasOilScaledImbPoint);
gasOilImbParamsHyst.setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]);
gasOilImbParamsHyst.finalize();
gasOilParams->setImbibitionParams(gasOilImbParamsHyst,
gasOilScaledImbInfo,
EclTwoPhaseSystemType::GasOil);
}
if (hasOil && hasWater) {
OilWaterEpsTwoPhaseParams oilWaterImbParamsHyst;
oilWaterImbParamsHyst.setConfig(oilWaterConfig);
oilWaterImbParamsHyst.setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]);
oilWaterImbParamsHyst.setScaledPoints(oilWaterScaledImbPoint);
oilWaterImbParamsHyst.setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]);
oilWaterImbParamsHyst.finalize();
oilWaterParams->setImbibitionParams(oilWaterImbParamsHyst,
oilWaterScaledImbInfo,
EclTwoPhaseSystemType::OilWater);
}
if (hasGas && hasWater && !hasOil) {
GasWaterEpsTwoPhaseParams gasWaterImbParamsHyst;
gasWaterImbParamsHyst.setConfig(gasWaterConfig);
gasWaterImbParamsHyst.setUnscaledPoints(gasWaterUnscaledPointsVector_[imbRegionIdx]);
gasWaterImbParamsHyst.setScaledPoints(gasWaterScaledImbPoint);
gasWaterImbParamsHyst.setEffectiveLawParams(gasWaterEffectiveParamVector_[imbRegionIdx]);
gasWaterImbParamsHyst.finalize();
gasWaterParams->setImbibitionParams(gasWaterImbParamsHyst,
gasWaterScaledImbInfo,
EclTwoPhaseSystemType::GasWater);
}
}
if (hasGas && hasOil)
gasOilParams->finalize();
if (hasOil && hasWater)
oilWaterParams->finalize();
if (hasGas && hasWater && !hasOil)
gasWaterParams->finalize();
initThreePhaseParams_(eclState,
materialLawParams_[elemIdx],
satRegionIdx,
oilWaterScaledEpsInfoDrainage_[elemIdx],
oilWaterParams,
gasOilParams,
gasWaterParams);
materialLawParams_[elemIdx].finalize();
}
InitParams initParams {*this, eclState, numCompressedElems};
initParams.run();
}
template<class TraitsT>
@@ -606,6 +354,28 @@ oilWaterScaledEpsPointsDrainage(unsigned elemIdx)
}
}
template<class TraitsT>
const typename EclMaterialLawManager<TraitsT>::MaterialLawParams& EclMaterialLawManager<TraitsT>::
materialLawParamsFunc_(unsigned elemIdx, FaceDir::DirEnum facedir) const
{
using Dir = FaceDir::DirEnum;
if (dirMaterialLawParams_) {
switch(facedir) {
case Dir::XPlus:
return dirMaterialLawParams_->materialLawParamsX_[elemIdx];
case Dir::YPlus:
return dirMaterialLawParams_->materialLawParamsY_[elemIdx];
case Dir::ZPlus:
return dirMaterialLawParams_->materialLawParamsZ_[elemIdx];
default:
throw std::runtime_error("Unexpected face direction");
}
}
else {
return materialLawParams_[elemIdx];
}
}
template<class TraitsT>
void EclMaterialLawManager<TraitsT>::
readGlobalEpsOptions_(const EclipseState& eclState)
@@ -662,499 +432,6 @@ readGlobalThreePhaseOptions_(const Runspec& runspec)
}
}
template<class TraitsT>
template <class Container>
void EclMaterialLawManager<TraitsT>::
readGasOilEffectiveParameters_(Container& dest,
const EclipseState& eclState,
unsigned satRegionIdx)
{
if (!hasGas || !hasOil)
// we don't read anything if either the gas or the oil phase is not active
return;
dest[satRegionIdx] = std::make_shared<GasOilEffectiveTwoPhaseParams>();
auto& effParams = *dest[satRegionIdx];
// the situation for the gas phase is complicated that all saturations are
// shifted by the connate water saturation.
const Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl;
const auto tolcrit = eclState.runspec().saturationFunctionControls()
.minimumRelpermMobilityThreshold();
const auto& tableManager = eclState.getTableManager();
switch (eclState.runspec().saturationFunctionControls().family()) {
case SatFuncControls::KeywordFamily::Family_I:
{
const TableContainer& sgofTables = tableManager.getSgofTables();
const TableContainer& slgofTables = tableManager.getSlgofTables();
if (!sgofTables.empty())
readGasOilEffectiveParametersSgof_(effParams, Swco, tolcrit,
sgofTables.getTable<SgofTable>(satRegionIdx));
else if (!slgofTables.empty())
readGasOilEffectiveParametersSlgof_(effParams, Swco, tolcrit,
slgofTables.getTable<SlgofTable>(satRegionIdx));
else if ( !tableManager.getSgofletTable().empty() ) {
const auto& letSgofTab = tableManager.getSgofletTable()[satRegionIdx];
const std::vector<Scalar> dum; // dummy arg to comform with existing interface
effParams.setApproach(SatCurveMultiplexerApproach::LET);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::LET>();
// 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<Scalar>& letCoeffsOil = {s_min_w, s_max_w,
static_cast<Scalar>(letSgofTab.l2_relperm),
static_cast<Scalar>(letSgofTab.e2_relperm),
static_cast<Scalar>(letSgofTab.t2_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsGas = {s_min_nw, s_max_nw,
static_cast<Scalar>(letSgofTab.l1_relperm),
static_cast<Scalar>(letSgofTab.e1_relperm),
static_cast<Scalar>(letSgofTab.t1_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsPc = {static_cast<Scalar>(letSgofTab.s2_residual),
static_cast<Scalar>(letSgofTab.s1_residual+Swco),
static_cast<Scalar>(letSgofTab.l_pc),
static_cast<Scalar>(letSgofTab.e_pc),
static_cast<Scalar>(letSgofTab.t_pc),
static_cast<Scalar>(letSgofTab.pcir_pc),
static_cast<Scalar>(letSgofTab.pct_pc)};
realParams.setPcnwSamples(letCoeffsPc, dum);
realParams.finalize();
}
break;
}
case SatFuncControls::KeywordFamily::Family_II:
{
const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
if (!hasWater) {
// oil and gas case
const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>( satRegionIdx );
readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable);
}
else {
const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable);
}
break;
}
case SatFuncControls::KeywordFamily::Undefined:
throw std::domain_error("No valid saturation keyword family specified");
}
}
template<class TraitsT>
void EclMaterialLawManager<TraitsT>::
readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const SgofTable& sgofTable)
{
// convert the saturations of the SGOF keyword from gas to oil saturations
std::vector<double> 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<SatCurveMultiplexerApproach::PiecewiseLinear>();
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<class TraitsT>
void EclMaterialLawManager<TraitsT>::
readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const SlgofTable& slgofTable)
{
// convert the saturations of the SLGOF keyword from "liquid" to oil saturations
std::vector<double> 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<SatCurveMultiplexerApproach::PiecewiseLinear>();
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<class TraitsT>
void EclMaterialLawManager<TraitsT>::
readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const Sof3Table& sof3Table,
const SgfnTable& sgfnTable)
{
// convert the saturations of the SGFN keyword from gas to oil saturations
std::vector<double> SoSamples(sgfnTable.numRows());
std::vector<double> SoColumn = sof3Table.getColumn("SO").vectorCopy();
for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
}
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROG")));
realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
realParams.finalize();
}
template<class TraitsT>
void EclMaterialLawManager<TraitsT>::
readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
const Scalar Swco,
const double tolcrit,
const Sof2Table& sof2Table,
const SgfnTable& sgfnTable)
{
// convert the saturations of the SGFN keyword from gas to oil saturations
std::vector<double> SoSamples(sgfnTable.numRows());
std::vector<double> SoColumn = sof2Table.getColumn("SO").vectorCopy();
for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
}
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
realParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO")));
realParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
realParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
realParams.finalize();
}
template<class TraitsT>
template <class Container>
void EclMaterialLawManager<TraitsT>::
readOilWaterEffectiveParameters_(Container& dest,
const EclipseState& eclState,
unsigned satRegionIdx)
{
if (!hasOil || !hasWater)
// we don't read anything if either the water or the oil phase is not active
return;
dest[satRegionIdx] = std::make_shared<OilWaterEffectiveTwoPhaseParams>();
const auto tolcrit = eclState.runspec().saturationFunctionControls()
.minimumRelpermMobilityThreshold();
const auto& tableManager = eclState.getTableManager();
auto& effParams = *dest[satRegionIdx];
switch (eclState.runspec().saturationFunctionControls().family()) {
case SatFuncControls::KeywordFamily::Family_I:
{
if (tableManager.hasTables("SWOF")) {
const auto& swofTable = tableManager.getSwofTables().getTable<SwofTable>(satRegionIdx);
const std::vector<double> SwColumn = swofTable.getColumn("SW").vectorCopy();
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
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<Scalar> dum; // dummy arg to conform with existing interface
effParams.setApproach(SatCurveMultiplexerApproach::LET);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::LET>();
// 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<Scalar>& letCoeffsWat = {s_min_w, s_max_w,
static_cast<Scalar>(letTab.l1_relperm),
static_cast<Scalar>(letTab.e1_relperm),
static_cast<Scalar>(letTab.t1_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsOil = {s_min_nw, s_max_nw,
static_cast<Scalar>(letTab.l2_relperm),
static_cast<Scalar>(letTab.e2_relperm),
static_cast<Scalar>(letTab.t2_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsPc = {static_cast<Scalar>(letTab.s1_residual),
static_cast<Scalar>(letTab.s2_residual),
static_cast<Scalar>(letTab.l_pc),
static_cast<Scalar>(letTab.e_pc),
static_cast<Scalar>(letTab.t_pc),
static_cast<Scalar>(letTab.pcir_pc),
static_cast<Scalar>(letTab.pct_pc)};
realParams.setPcnwSamples(letCoeffsPc, dum);
realParams.finalize();
}
break;
}
case SatFuncControls::KeywordFamily::Family_II:
{
const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
const std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
if (!hasGas) {
const auto& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>(satRegionIdx);
// convert the saturations of the SOF2 keyword from oil to water saturations
std::vector<double> SwSamples(sof2Table.numRows());
for (size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sof2Table.get("SO", sampleIdx);
realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO")));
} else {
const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
// convert the saturations of the SOF3 keyword from oil to water saturations
std::vector<double> SwSamples(sof3Table.numRows());
for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx);
realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROW")));
}
realParams.finalize();
break;
}
case SatFuncControls::KeywordFamily::Undefined:
throw std::domain_error("No valid saturation keyword family specified");
}
}
template<class TraitsT>
template <class Container>
void EclMaterialLawManager<TraitsT>::
readGasWaterEffectiveParameters_(Container& dest,
const EclipseState& eclState,
unsigned satRegionIdx)
{
if (!hasGas || !hasWater || hasOil)
// we don't read anything if either the gas or the water phase is not active or if oil is present
return;
dest[satRegionIdx] = std::make_shared<GasWaterEffectiveTwoPhaseParams>();
auto& effParams = *dest[satRegionIdx];
const auto tolcrit = eclState.runspec().saturationFunctionControls()
.minimumRelpermMobilityThreshold();
const auto& tableManager = eclState.getTableManager();
switch (eclState.runspec().saturationFunctionControls().family()) {
case SatFuncControls::KeywordFamily::Family_I:
{
throw std::domain_error("Saturation keyword family I is not applicable for a gas-water system");
}
case SatFuncControls::KeywordFamily::Family_II:
{
//Todo: allow also for Sgwfn table input as alternative to Sgfn and Swfn table input
const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
const SwfnTable& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>( satRegionIdx );
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
std::vector<double> SwSamples(sgfnTable.numRows());
for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx);
realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
//Capillary pressure is read from SWFN.
//For gas-water system the capillary pressure column values are set to 0 in SGFN
realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
realParams.finalize();
break;
}
case SatFuncControls::KeywordFamily::Undefined:
throw std::domain_error("No valid saturation keyword family specified");
}
}
template<class TraitsT>
template <class Container>
void EclMaterialLawManager<TraitsT>::
readGasOilUnscaledPoints_(Container& dest,
std::shared_ptr<EclEpsConfig> config,
const EclipseState& /* eclState */,
unsigned satRegionIdx)
{
if (!hasGas || !hasOil)
// we don't read anything if either the gas or the oil phase is not active
return;
dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx],
*config,
EclTwoPhaseSystemType::GasOil);
}
template<class TraitsT>
template <class Container>
void EclMaterialLawManager<TraitsT>::
readOilWaterUnscaledPoints_(Container& dest,
std::shared_ptr<EclEpsConfig> config,
const EclipseState& /* eclState */,
unsigned satRegionIdx)
{
if (!hasOil || !hasWater)
// we don't read anything if either the water or the oil phase is not active
return;
dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx],
*config,
EclTwoPhaseSystemType::OilWater);
}
template<class TraitsT>
template <class Container>
void EclMaterialLawManager<TraitsT>::
readGasWaterUnscaledPoints_(Container& dest,
std::shared_ptr<EclEpsConfig> config,
const EclipseState& /* eclState */,
unsigned satRegionIdx)
{
if (hasOil)
// we don't read anything if oil phase is active
return;
dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx],
*config,
EclTwoPhaseSystemType::GasWater);
}
template<class TraitsT>
std::tuple<EclEpsScalingPointsInfo<typename TraitsT::Scalar>,
EclEpsScalingPoints<typename TraitsT::Scalar>>
EclMaterialLawManager<TraitsT>::
readScaledPoints_(const EclEpsConfig& config,
const EclipseState& eclState,
const EclEpsGridProperties& epsGridProperties,
unsigned elemIdx,
EclTwoPhaseSystemType type)
{
unsigned satRegionIdx = epsGridProperties.satRegion( elemIdx );
EclEpsScalingPointsInfo<Scalar> destInfo(unscaledEpsInfo_[satRegionIdx]);
destInfo.extractScaled(eclState, epsGridProperties, elemIdx);
EclEpsScalingPoints<Scalar> destPoint;
destPoint.init(destInfo, config, type);
return {destInfo, destPoint};
}
template<class TraitsT>
void EclMaterialLawManager<TraitsT>::
initThreePhaseParams_(const EclipseState& /* eclState */,
MaterialLawParams& materialParams,
unsigned satRegionIdx,
const EclEpsScalingPointsInfo<Scalar>& epsInfo,
std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams,
std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams,
std::shared_ptr<GasWaterTwoPhaseHystParams> gasWaterParams)
{
materialParams.setApproach(threePhaseApproach_);
switch (materialParams.approach()) {
case EclMultiplexerApproach::Stone1: {
auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::Stone1>();
realParams.setGasOilParams(gasOilParams);
realParams.setOilWaterParams(oilWaterParams);
realParams.setSwl(epsInfo.Swl);
if (!stoneEtas.empty()) {
realParams.setEta(stoneEtas[satRegionIdx]);
}
else
realParams.setEta(1.0);
realParams.finalize();
break;
}
case EclMultiplexerApproach::Stone2: {
auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::Stone2>();
realParams.setGasOilParams(gasOilParams);
realParams.setOilWaterParams(oilWaterParams);
realParams.setSwl(epsInfo.Swl);
realParams.finalize();
break;
}
case EclMultiplexerApproach::Default: {
auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::Default>();
realParams.setGasOilParams(gasOilParams);
realParams.setOilWaterParams(oilWaterParams);
realParams.setSwl(epsInfo.Swl);
realParams.finalize();
break;
}
case EclMultiplexerApproach::TwoPhase: {
auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::TwoPhase>();
realParams.setGasOilParams(gasOilParams);
realParams.setOilWaterParams(oilWaterParams);
realParams.setGasWaterParams(gasWaterParams);
realParams.setApproach(twoPhaseApproach_);
realParams.finalize();
break;
}
case EclMultiplexerApproach::OnePhase: {
// Nothing to do, no parameters.
break;
}
}
}
template class EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>;
template class EclMaterialLawManager<ThreePhaseMaterialTraits<float,0,1,2>>;

291
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp Normal file
View File

@@ -0,0 +1,291 @@
// -*- 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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp>
namespace Opm {
/* constructors*/
template <class Traits>
EclMaterialLawManager<Traits>::InitParams::HystParams::
HystParams(EclMaterialLawManager<Traits>::InitParams& init_params) :
init_params_{init_params}, parent_{init_params_.parent_},
eclState_{init_params_.eclState_}
{
gasOilParams_ = std::make_shared<GasOilTwoPhaseHystParams>();
oilWaterParams_ = std::make_shared<OilWaterTwoPhaseHystParams>();
gasWaterParams_ = std::make_shared<GasWaterTwoPhaseHystParams>();
}
/* public methods, alphabetically sorted */
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
finalize()
{
if (hasGasOil_())
this->gasOilParams_->finalize();
if (hasOilWater_())
this->oilWaterParams_->finalize();
if (hasGasWater_())
this->gasWaterParams_->finalize();
}
template <class Traits>
std::shared_ptr<typename EclMaterialLawManager<Traits>::GasOilTwoPhaseHystParams>
EclMaterialLawManager<Traits>::InitParams::HystParams::
getGasOilParams()
{
return gasOilParams_;
}
template <class Traits>
std::shared_ptr<typename EclMaterialLawManager<Traits>::OilWaterTwoPhaseHystParams>
EclMaterialLawManager<Traits>::InitParams::HystParams::
getOilWaterParams()
{
return oilWaterParams_;
}
template <class Traits>
std::shared_ptr<typename EclMaterialLawManager<Traits>::GasWaterTwoPhaseHystParams>
EclMaterialLawManager<Traits>::InitParams::HystParams::
getGasWaterParams()
{
return gasWaterParams_;
}
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setConfig()
{
this->gasOilParams_->setConfig(this->parent_.hysteresisConfig_);
this->oilWaterParams_->setConfig(this->parent_.hysteresisConfig_);
this->gasWaterParams_->setConfig(this->parent_.hysteresisConfig_);
} // namespace Opm
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setDrainageParamsGasWater(unsigned elemIdx, unsigned satRegionIdx)
{
if (hasGasWater_()) {
auto [gasWaterScaledInfo, gasWaterScaledPoints]
= readScaledEpsPointsDrainage_(elemIdx, EclTwoPhaseSystemType::GasWater);
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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setDrainageParamsOilGas(unsigned elemIdx, unsigned satRegionIdx)
{
if (hasGasOil_()) {
auto [gasOilScaledInfo, gasOilScaledPoints]
= readScaledEpsPointsDrainage_(elemIdx, EclTwoPhaseSystemType::GasOil);
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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setDrainageParamsOilWater(unsigned elemIdx, unsigned satRegionIdx)
{
// 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);
// 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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setImbibitionParamsGasWater(unsigned elemIdx, unsigned imbRegionIdx)
{
if (hasGasWater_()) {
auto [gasWaterScaledInfo, gasWaterScaledPoints]
= readScaledEpsPointsImbibition_(elemIdx, EclTwoPhaseSystemType::GasWater);
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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setImbibitionParamsOilGas(unsigned elemIdx, unsigned imbRegionIdx)
{
if (hasGasOil_()) {
auto [gasOilScaledInfo, gasOilScaledPoints]
= readScaledEpsPointsImbibition_(elemIdx, EclTwoPhaseSystemType::GasOil);
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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::HystParams::
setImbibitionParamsOilWater(unsigned elemIdx, unsigned imbRegionIdx)
{
if (hasOilWater_()) {
auto [oilWaterScaledInfo, oilWaterScaledPoints]
= readScaledEpsPointsImbibition_(elemIdx, EclTwoPhaseSystemType::OilWater);
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 <class Traits>
bool
EclMaterialLawManager<Traits>::InitParams::HystParams::
hasGasOil_()
{
return this->parent_.hasGas && this->parent_.hasOil;
}
template <class Traits>
bool
EclMaterialLawManager<Traits>::InitParams::HystParams::
hasGasWater_()
{
return this->parent_.hasGas && this->parent_.hasWater && !this->parent_.hasOil;
}
template <class Traits>
bool
EclMaterialLawManager<Traits>::InitParams::HystParams::
hasOilWater_()
{
return this->parent_.hasOil && this->parent_.hasWater;
}
template <class Traits>
std::tuple<
EclEpsScalingPointsInfo<typename EclMaterialLawManager<Traits>::Scalar>,
EclEpsScalingPoints<typename EclMaterialLawManager<Traits>::Scalar>
>
EclMaterialLawManager<Traits>::InitParams::HystParams::
readScaledEpsPoints_(const EclEpsGridProperties& epsGridProperties, unsigned elemIdx, EclTwoPhaseSystemType type)
{
const EclEpsConfig& config = *(this->parent_.gasOilConfig_);
unsigned satRegionIdx = epsGridProperties.satRegion( elemIdx );
// Copy-construct a new instance of EclEpsScalingPointsInfo
EclEpsScalingPointsInfo<Scalar> 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, elemIdx);
EclEpsScalingPoints<Scalar> destPoint;
destPoint.init(destInfo, config, type);
return {destInfo, destPoint};
}
template <class Traits>
std::tuple<
EclEpsScalingPointsInfo<typename EclMaterialLawManager<Traits>::Scalar>,
EclEpsScalingPoints<typename EclMaterialLawManager<Traits>::Scalar>
>
EclMaterialLawManager<Traits>::InitParams::HystParams::
readScaledEpsPointsDrainage_(unsigned elemIdx, EclTwoPhaseSystemType type)
{
const auto& epsGridProperties = this->init_params_.epsGridProperties_;
return readScaledEpsPoints_(*epsGridProperties, elemIdx, type);
}
template <class Traits>
std::tuple<
EclEpsScalingPointsInfo<typename EclMaterialLawManager<Traits>::Scalar>,
EclEpsScalingPoints<typename EclMaterialLawManager<Traits>::Scalar>
>
EclMaterialLawManager<Traits>::InitParams::HystParams::
readScaledEpsPointsImbibition_(unsigned elemIdx, EclTwoPhaseSystemType type)
{
const auto& epsGridProperties = this->init_params_.epsImbGridProperties_;
return readScaledEpsPoints_(*epsGridProperties, elemIdx, type);
}
// Make some actual code, by realizing the previously defined templated class
template class EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>::InitParams::HystParams;
template class EclMaterialLawManager<ThreePhaseMaterialTraits<float,0,1,2>>::InitParams::HystParams;
} // namespace Opm

342
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerInitParams.cpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/fluidmatrixinteractions/EclEpsGridProperties.hpp>
namespace Opm {
/* constructors*/
template <class Traits>
EclMaterialLawManager<Traits>::InitParams::
InitParams(EclMaterialLawManager<Traits>& 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<EclEpsGridProperties>(this->eclState_, /*useImbibition=*/true);
}
this->epsGridProperties_
= std::make_unique<EclEpsGridProperties>(this->eclState_, /*useImbibition=*/false);
}
/* public methods */
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
run() {
readUnscaledEpsPointsVectors_();
readEffectiveParameters_();
initSatnumRegionArray_();
copySatnumArrays_();
initOilWaterScaledEpsInfo_();
initMaterialLawParamVectors_();
std::vector<std::vector<int>*> satnumArray;
std::vector<std::vector<int>*> imbnumArray;
std::vector<std::vector<MaterialLawParams>*> mlpArray;
initArrays_(satnumArray, imbnumArray, mlpArray);
auto num_arrays = mlpArray.size();
for (unsigned i=0; i<num_arrays; i++) {
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();
hystParams.setDrainageParamsOilGas(elemIdx, satRegionIdx);
hystParams.setDrainageParamsOilWater(elemIdx, satRegionIdx);
hystParams.setDrainageParamsGasWater(elemIdx, satRegionIdx);
if (this->parent_.enableHysteresis()) {
unsigned imbRegionIdx = imbRegion_(*imbnumArray[i], elemIdx);
hystParams.setImbibitionParamsOilGas(elemIdx, imbRegionIdx);
hystParams.setImbibitionParamsOilWater(elemIdx, imbRegionIdx);
hystParams.setImbibitionParamsGasWater(elemIdx, imbRegionIdx);
}
hystParams.finalize();
initThreePhaseParams_(hystParams, (*mlpArray[i])[elemIdx], satRegionIdx, elemIdx);
}
}
}
/* private methods alphabetically sorted*/
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
copySatnumArrays_()
{
copyIntArray_(this->parent_.krnumXArray_, "KRNUMX");
copyIntArray_(this->parent_.krnumYArray_, "KRNUMY");
copyIntArray_(this->parent_.krnumZArray_, "KRNUMZ");
copyIntArray_(this->parent_.imbnumXArray_, "IMBNUMX");
copyIntArray_(this->parent_.imbnumYArray_, "IMBNUMY");
copyIntArray_(this->parent_.imbnumZArray_, "IMBNUMZ");
// 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");
assert(this->numCompressedElems_ == this->parent_.satnumRegionArray_.size());
assert(!this->parent_.enableHysteresis() || this->numCompressedElems_ == this->parent_.imbnumRegionArray_.size());
}
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
copyIntArray_(std::vector<int>& dest, const std::string keyword)
{
if (this->eclState_.fieldProps().has_int(keyword)) {
dest.resize(this->numCompressedElems_);
const auto& satnumRawData = this->eclState_.fieldProps().get_int(keyword);
for (unsigned elemIdx = 0; elemIdx < this->numCompressedElems_; ++elemIdx) {
dest[elemIdx] = satnumRawData[elemIdx] - 1;
}
}
}
template <class Traits>
unsigned
EclMaterialLawManager<Traits>::InitParams::
imbRegion_(std::vector<int>& array, unsigned elemIdx)
{
std::vector<int>& default_vec = this->parent_.imbnumRegionArray_;
return satOrImbRegion_(array, default_vec, elemIdx);
}
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
initArrays_(
std::vector<std::vector<int>*>& satnumArray,
std::vector<std::vector<int>*>& imbnumArray,
std::vector<std::vector<MaterialLawParams>*>& 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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
initMaterialLawParamVectors_()
{
this->parent_.materialLawParams_.resize(this->numCompressedElems_);
if (this->parent_.hasDirectionalImbnum() || this->parent_.hasDirectionalRelperms()) {
this->parent_.dirMaterialLawParams_
= std::make_unique<DirectionalMaterialLawParams<MaterialLawParams>>(this->numCompressedElems_);
}
}
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
initOilWaterScaledEpsInfo_()
{
// This vector will be updated in the hystParams.setDrainageOilWater() in the run() method
this->parent_.oilWaterScaledEpsInfoDrainage_.resize(this->numCompressedElems_);
}
template <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
initSatnumRegionArray_()
{
// 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")) {
const auto& satnumRawData = this->eclState_.fieldProps().get_int("SATNUM");
for (unsigned elemIdx = 0; elemIdx < this->numCompressedElems_; ++elemIdx) {
satnumArray[elemIdx] = satnumRawData[elemIdx] - 1;
}
}
else {
std::fill(satnumArray.begin(), satnumArray.end(), 0);
}
}
template <class Traits>
void
EclMaterialLawManager<Traits>::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<EclMultiplexerApproach::Stone1>();
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<EclMultiplexerApproach::Stone2>();
realParams.setGasOilParams(gasOilParams);
realParams.setOilWaterParams(oilWaterParams);
realParams.setSwl(epsInfo.Swl);
realParams.finalize();
break;
}
case EclMultiplexerApproach::Default: {
auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::Default>();
realParams.setGasOilParams(gasOilParams);
realParams.setOilWaterParams(oilWaterParams);
realParams.setSwl(epsInfo.Swl);
realParams.finalize();
break;
}
case EclMultiplexerApproach::TwoPhase: {
auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::TwoPhase>();
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 <class Traits>
void
EclMaterialLawManager<Traits>::InitParams::
readEffectiveParameters_()
{
ReadEffectiveParams effectiveReader {*this};
// populates effective parameter vectors in the parent class (EclMaterialManager)
effectiveReader.read();
}
template <class Traits>
void
EclMaterialLawManager<Traits>::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 <class Traits>
template <class Container>
void
EclMaterialLawManager<Traits>::InitParams::
readUnscaledEpsPoints_(Container& dest, std::shared_ptr<EclEpsConfig> 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<EclEpsScalingPoints<Scalar> >();
dest[satRegionIdx]->init(this->parent_.unscaledEpsInfo_[satRegionIdx], *config, system_type);
}
}
template <class Traits>
unsigned
EclMaterialLawManager<Traits>::InitParams::
satRegion_(std::vector<int>& array, unsigned elemIdx)
{
std::vector<int>& default_vec = this->parent_.satnumRegionArray_;
return satOrImbRegion_(array, default_vec, elemIdx);
}
template <class Traits>
unsigned
EclMaterialLawManager<Traits>::InitParams::
satOrImbRegion_(std::vector<int>& array, std::vector<int>& default_vec, unsigned elemIdx)
{
int value;
if (array.size() > 0) {
value = array[elemIdx];
}
else { // use default value
value = default_vec[elemIdx];
}
return static_cast<unsigned>(value);
}
// Make some actual code, by realizing the previously defined templated class
template class EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>::InitParams;
template class EclMaterialLawManager<ThreePhaseMaterialTraits<float,0,1,2>>::InitParams;
} // namespace Opm

414
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerReadEffectiveParams.cpp Normal file
View File

@@ -0,0 +1,414 @@
// -*- 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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SgfnTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SgofTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SlgofTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Sof2Table.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Sof3Table.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SwfnTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SwofTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
namespace Opm {
/* constructors*/
template <class Traits>
EclMaterialLawManager<Traits>::InitParams::ReadEffectiveParams::
ReadEffectiveParams(EclMaterialLawManager<Traits>::InitParams& init_params) :
init_params_{init_params}, parent_{init_params_.parent_},
eclState_{init_params_.eclState_}
{
}
/* public methods */
template <class Traits>
void
EclMaterialLawManager<Traits>::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 <class Traits>
std::vector<double>
EclMaterialLawManager<Traits>::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 <class Traits>
void
EclMaterialLawManager<Traits>::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<GasOilEffectiveTwoPhaseParams>();
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.getTable<SgofTable>(satRegionIdx));
else if (!slgofTables.empty())
readGasOilSlgof_(effParams, Swco, tolcrit, slgofTables.getTable<SlgofTable>(satRegionIdx));
else if ( !tableManager.getSgofletTable().empty() ) {
const auto& letSgofTab = tableManager.getSgofletTable()[satRegionIdx];
const std::vector<Scalar> dum; // dummy arg to comform with existing interface
effParams.setApproach(SatCurveMultiplexerApproach::LET);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::LET>();
// 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<Scalar>& letCoeffsOil = {s_min_w, s_max_w,
static_cast<Scalar>(letSgofTab.l2_relperm),
static_cast<Scalar>(letSgofTab.e2_relperm),
static_cast<Scalar>(letSgofTab.t2_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsGas = {s_min_nw, s_max_nw,
static_cast<Scalar>(letSgofTab.l1_relperm),
static_cast<Scalar>(letSgofTab.e1_relperm),
static_cast<Scalar>(letSgofTab.t1_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsPc = {static_cast<Scalar>(letSgofTab.s2_residual),
static_cast<Scalar>(letSgofTab.s1_residual+Swco),
static_cast<Scalar>(letSgofTab.l_pc),
static_cast<Scalar>(letSgofTab.e_pc),
static_cast<Scalar>(letSgofTab.t_pc),
static_cast<Scalar>(letSgofTab.pcir_pc),
static_cast<Scalar>(letSgofTab.pct_pc)};
realParams.setPcnwSamples(letCoeffsPc, dum);
realParams.finalize();
}
break;
}
case SatFuncControls::KeywordFamily::Family_II:
{
const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
if (!this->parent_.hasWater) {
// oil and gas case
const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>( satRegionIdx );
readGasOilFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable, /*columnName=*/"KRO");
}
else {
const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
readGasOilFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable, /* columnName=*/"KROG");
}
break;
}
case SatFuncControls::KeywordFamily::Undefined:
throw std::domain_error("No valid saturation keyword family specified");
}
}
template <class Traits>
template <class TableType>
void
EclMaterialLawManager<Traits>::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<double> SoSamples(sgfnTable.numRows());
std::vector<double> 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<SatCurveMultiplexerApproach::PiecewiseLinear>();
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 <class Traits>
void
EclMaterialLawManager<Traits>::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<double> 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<SatCurveMultiplexerApproach::PiecewiseLinear>();
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 <class Traits>
void
EclMaterialLawManager<Traits>::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<double> 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<SatCurveMultiplexerApproach::PiecewiseLinear>();
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 <class Traits>
void
EclMaterialLawManager<Traits>::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<GasWaterEffectiveTwoPhaseParams>();
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:
{
//Todo: allow also for Sgwfn table input as alternative to Sgfn and Swfn table input
const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
const SwfnTable& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>( satRegionIdx );
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
std::vector<double> SwSamples(sgfnTable.numRows());
for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx);
realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
//Capillary pressure is read from SWFN.
//For gas-water system the capillary pressure column values are set to 0 in SGFN
realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
realParams.finalize();
break;
}
case SatFuncControls::KeywordFamily::Undefined:
throw std::domain_error("No valid saturation keyword family specified");
}
}
template <class Traits>
void
EclMaterialLawManager<Traits>::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<OilWaterEffectiveTwoPhaseParams>();
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().getTable<SwofTable>(satRegionIdx);
const std::vector<double> SwColumn = swofTable.getColumn("SW").vectorCopy();
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
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<Scalar> dum; // dummy arg to conform with existing interface
effParams.setApproach(SatCurveMultiplexerApproach::LET);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::LET>();
// 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<Scalar>& letCoeffsWat = {s_min_w, s_max_w,
static_cast<Scalar>(letTab.l1_relperm),
static_cast<Scalar>(letTab.e1_relperm),
static_cast<Scalar>(letTab.t1_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsOil = {s_min_nw, s_max_nw,
static_cast<Scalar>(letTab.l2_relperm),
static_cast<Scalar>(letTab.e2_relperm),
static_cast<Scalar>(letTab.t2_relperm),
static_cast<Scalar>(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<Scalar>& letCoeffsPc = {static_cast<Scalar>(letTab.s1_residual),
static_cast<Scalar>(letTab.s2_residual),
static_cast<Scalar>(letTab.l_pc),
static_cast<Scalar>(letTab.e_pc),
static_cast<Scalar>(letTab.t_pc),
static_cast<Scalar>(letTab.pcir_pc),
static_cast<Scalar>(letTab.pct_pc)};
realParams.setPcnwSamples(letCoeffsPc, dum);
realParams.finalize();
}
break;
}
case SatFuncControls::KeywordFamily::Family_II:
{
const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
const std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
effParams.setApproach(SatCurveMultiplexerApproach::PiecewiseLinear);
auto& realParams = effParams.template getRealParams<SatCurveMultiplexerApproach::PiecewiseLinear>();
realParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
realParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
if (!this->parent_.hasGas) {
const auto& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>(satRegionIdx);
// convert the saturations of the SOF2 keyword from oil to water saturations
std::vector<double> SwSamples(sof2Table.numRows());
for (size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sof2Table.get("SO", sampleIdx);
realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO")));
} else {
const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
// convert the saturations of the SOF3 keyword from oil to water saturations
std::vector<double> SwSamples(sof3Table.numRows());
for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx);
realParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROW")));
}
realParams.finalize();
break;
}
case SatFuncControls::KeywordFamily::Undefined:
throw std::domain_error("No valid saturation keyword family specified");
}
}
// Make some actual code, by realizing the previously defined templated class
template class EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>::InitParams::ReadEffectiveParams;
template class EclMaterialLawManager<ThreePhaseMaterialTraits<float,0,1,2>>::InitParams::ReadEffectiveParams;
} // namespace Opm