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

Merge pull request #3314 from akva2/more_cpp_pvt

Browse Source 
Move more pvt code to compile units
This commit is contained in:
Arne Morten Kvarving 2023-01-05 10:40:43 +01:00 committed by GitHub
commit e440b6b008
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
15 changed files with 830 additions and 680 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
CMakeLists_files.cmake
View File

@ -52,6 +52,13 @@ list (APPEND MAIN_SOURCE_FILES
src/opm/material/components/CO2.cpp src/opm/material/components/CO2.cpp
src/opm/material/densead/Evaluation.cpp src/opm/material/densead/Evaluation.cpp
src/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.cpp src/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.cpp
src/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/DryGasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/DryHumidGasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/SolventPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/WetGasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/WetHumidGasPvt.cpp
) )
if(ENABLE_ECL_INPUT) if(ENABLE_ECL_INPUT)
list(APPEND MAIN_SOURCE_FILES list(APPEND MAIN_SOURCE_FILES
@ -253,19 +260,12 @@ if(ENABLE_ECL_INPUT)
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.cpp src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.cpp
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.cpp src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.cpp src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/DryGasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/DryHumidGasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.cpp src/opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.cpp
src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp
src/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.cpp src/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.cpp
src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp
src/opm/material/fluidsystems/blackoilpvt/SolventPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.cpp src/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.cpp
src/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.cpp src/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.cpp
src/opm/material/fluidsystems/blackoilpvt/WetGasPvt.cpp
src/opm/material/fluidsystems/blackoilpvt/WetHumidGasPvt.cpp
) )

35
opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp
View File

@ -53,13 +53,7 @@ public:
void initFromState(const EclipseState& eclState, const Schedule&); void initFromState(const EclipseState& eclState, const Schedule&);
#endif // HAVE_ECL_INPUT #endif // HAVE_ECL_INPUT
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
oilReferenceDensity_.resize(numRegions);
inverseOilB_.resize(numRegions);
inverseOilBMu_.resize(numRegions);
oilMu_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference densities of all fluids for a given PVT region * \brief Initialize the reference densities of all fluids for a given PVT region
@ -96,32 +90,7 @@ public:
/*! /*!
* \brief Finish initializing the oil phase PVT properties. * \brief Finish initializing the oil phase PVT properties.
*/ */
void initEnd() void initEnd();
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = oilMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the oil
// formation volume factor and the oil viscosity
const auto& oilMu = oilMu_[regionIdx];
const auto& invOilB = inverseOilB_[regionIdx];
assert(oilMu.numSamples() == invOilB.numSamples());
std::vector<Scalar> invBMuColumn;
std::vector<Scalar> pressureColumn;
invBMuColumn.resize(oilMu.numSamples());
pressureColumn.resize(oilMu.numSamples());
for (unsigned pIdx = 0; pIdx < oilMu.numSamples(); ++pIdx) {
pressureColumn[pIdx] = invOilB.xAt(pIdx);
invBMuColumn[pIdx] = invOilB.valueAt(pIdx)*1/oilMu.valueAt(pIdx);
}
inverseOilBMu_[regionIdx].setXYArrays(pressureColumn.size(),
pressureColumn,
invBMuColumn);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.

44
opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp
View File

@ -61,14 +61,7 @@ public:
void initFromState(const EclipseState& eclState, const Schedule&); void initFromState(const EclipseState& eclState, const Schedule&);
#endif #endif
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions);
inverseGasBMu_.resize(numRegions);
gasMu_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference densities of all fluids for a given PVT region * \brief Initialize the reference densities of all fluids for a given PVT region
@ -103,42 +96,13 @@ public:
* *
* \param samplePoints A container of \f$(p_g, B_g)\f$ values * \param samplePoints A container of \f$(p_g, B_g)\f$ values
*/ */
void setGasFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints) void setGasFormationVolumeFactor(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
SamplingPoints tmp(samplePoints);
auto it = tmp.begin();
const auto& endIt = tmp.end();
for (; it != endIt; ++ it)
std::get<1>(*it) = 1.0/std::get<1>(*it);
inverseGasB_[regionIdx].setContainerOfTuples(tmp);
assert(inverseGasB_[regionIdx].monotonic());
}
/*! /*!
* \brief Finish initializing the oil phase PVT properties. * \brief Finish initializing the oil phase PVT properties.
*/ */
void initEnd() void initEnd();
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMu = gasMu_[regionIdx];
const auto& invGasB = inverseGasB_[regionIdx];
assert(gasMu.numSamples() == invGasB.numSamples());
std::vector<Scalar> pressureValues(gasMu.numSamples());
std::vector<Scalar> invGasBMuValues(gasMu.numSamples());
for (unsigned pIdx = 0; pIdx < gasMu.numSamples(); ++pIdx) {
pressureValues[pIdx] = invGasB.xAt(pIdx);
invGasBMuValues[pIdx] = invGasB.valueAt(pIdx) * (1.0/gasMu.valueAt(pIdx));
}
inverseGasBMu_[regionIdx].setXYContainers(pressureValues, invGasBMuValues);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.

126
opm/material/fluidsystems/blackoilpvt/DryHumidGasPvt.hpp
View File

@ -72,18 +72,7 @@ private:
public: public:
#endif // HAVE_ECL_INPUT #endif // HAVE_ECL_INPUT
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
waterReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedWaterVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference densities of all fluids for a given PVT region * \brief Initialize the reference densities of all fluids for a given PVT region
@ -91,11 +80,7 @@ public:
void setReferenceDensities(unsigned regionIdx, void setReferenceDensities(unsigned regionIdx,
Scalar /*rhoRefOil*/, Scalar /*rhoRefOil*/,
Scalar rhoRefGas, Scalar rhoRefGas,
Scalar rhoRefWater) Scalar rhoRefWater);
{
waterReferenceDensity_[regionIdx] = rhoRefWater;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
/*! /*!
* \brief Initialize the function for the oil vaporization factor \f$R_v\f$ * \brief Initialize the function for the oil vaporization factor \f$R_v\f$
@ -135,85 +120,13 @@ public:
* requires the viscosity of oil-saturated gas (which only depends on pressure) while * requires the viscosity of oil-saturated gas (which only depends on pressure) while
* there is assumed to be no dependence on the gas mass fraction... * there is assumed to be no dependence on the gas mass fraction...
*/ */
void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints& samplePoints ) void setSaturatedGasViscosity(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
constexpr const Scalar RwMin = 0.0;
Scalar RwMax = waterVaporizationFac.eval(saturatedWaterVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
Scalar poMin = samplePoints.front().first;
Scalar poMax = samplePoints.back().first;
constexpr const size_t nRw = 20;
size_t nP = samplePoints.size()*2;
TabulatedOneDFunction mugTable;
mugTable.setContainerOfTuples(samplePoints);
// calculate a table of estimated densities depending on pressure and gas mass
// fraction
for (size_t RwIdx = 0; RwIdx < nRw; ++RwIdx) {
Scalar Rw = RwMin + (RwMax - RwMin)*RwIdx/nRw;
gasMu_[regionIdx].appendXPos(Rw);
for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
Scalar mug = mugTable.eval(pg, /*extrapolate=*/true);
gasMu_[regionIdx].appendSamplePoint(RwIdx, pg, mug);
}
}
}
/*! /*!
* \brief Finish initializing the gas phase PVT properties. * \brief Finish initializing the gas phase PVT properties.
*/ */
void initEnd() void initEnd();
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMu = gasMu_[regionIdx];
const auto& invGasB = inverseGasB_[regionIdx];
assert(gasMu.numX() == invGasB.numX());
auto& invGasBMu = inverseGasBMu_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
invGasBMu.appendXPos(gasMu.xAt(pIdx));
assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
size_t numRw = gasMu.numY(pIdx);
for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
invGasBMu.appendSamplePoint(pIdx,
gasMu.yAt(pIdx, RwIdx),
invGasB.valueAt(pIdx, RwIdx)
/ gasMu.valueAt(pIdx, RwIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rw value)
satPressuresArray.push_back(gasMu.xAt(pIdx));
invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRw - 1));
invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRw - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
updateSaturationPressure_(regionIdx);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.
@ -441,34 +354,7 @@ public:
} }
private: private:
void updateSaturationPressure_(unsigned regionIdx) void updateSaturationPressure_(unsigned regionIdx);
{
typedef std::pair<Scalar, Scalar> Pair;
const auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rw
size_t n = waterVaporizationFac.numSamples();
Scalar delta = (waterVaporizationFac.xMax() - waterVaporizationFac.xMin())/Scalar(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rw = 0;
for (size_t i = 0; i <= n; ++ i) {
Scalar pSat = waterVaporizationFac.xMin() + Scalar(i)*delta;
Rw = saturatedWaterVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
Pair val(Rw, pSat);
pSatSamplePoints.push_back(val);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
std::vector<Scalar> gasReferenceDensity_; std::vector<Scalar> gasReferenceDensity_;
std::vector<Scalar> waterReferenceDensity_; std::vector<Scalar> waterReferenceDensity_;

148
opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp
View File

@ -70,19 +70,7 @@ private:
public: public:
#endif // HAVE_ECL_INPUT #endif // HAVE_ECL_INPUT
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseOilBTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
inverseOilBMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
inverseSaturatedOilBTable_.resize(numRegions);
inverseSaturatedOilBMuTable_.resize(numRegions);
oilMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
saturatedOilMuTable_.resize(numRegions);
saturatedGasDissolutionFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference densities of all fluids for a given PVT region * \brief Initialize the reference densities of all fluids for a given PVT region
@ -90,11 +78,7 @@ public:
void setReferenceDensities(unsigned regionIdx, void setReferenceDensities(unsigned regionIdx,
Scalar rhoRefOil, Scalar rhoRefOil,
Scalar rhoRefGas, Scalar rhoRefGas,
Scalar /*rhoRefWater*/) Scalar /*rhoRefWater*/);
{
oilReferenceDensity_[regionIdx] = rhoRefOil;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
/*! /*!
* \brief Initialize the function for the gas dissolution factor \f$R_s\f$ * \brief Initialize the function for the gas dissolution factor \f$R_s\f$
@ -113,29 +97,8 @@ public:
* only depends on pressure) while the dependence on the gas mass fraction is * only depends on pressure) while the dependence on the gas mass fraction is
* guesstimated. * guesstimated.
*/ */
void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints) void setSaturatedOilFormationVolumeFactor(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
constexpr const Scalar T = 273.15 + 15.56; // [K]
auto& invOilB = inverseOilBTable_[regionIdx];
updateSaturationPressure_(regionIdx);
// calculate a table of estimated densities of undersatured gas
for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
Scalar p1 = std::get<0>(samplePoints[pIdx]);
Scalar p2 = p1 * 2.0;
Scalar Bo1 = std::get<1>(samplePoints[pIdx]);
Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
Scalar Bo2 = Bo1/(1.0 + (p2 - p1)*drhoo_dp);
Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
invOilB.appendXPos(Rs);
invOilB.appendSamplePoint(pIdx, p1, 1.0/Bo1);
invOilB.appendSamplePoint(pIdx, p2, 1.0/Bo2);
}
}
/*! /*!
* \brief Initialize the function for the oil formation volume factor * \brief Initialize the function for the oil formation volume factor
@ -167,79 +130,13 @@ public:
* requires the viscosity of gas-saturated oil (which only depends on pressure) while * requires the viscosity of gas-saturated oil (which only depends on pressure) while
* there is assumed to be no dependence on the gas mass fraction... * there is assumed to be no dependence on the gas mass fraction...
*/ */
void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints& samplePoints) void setSaturatedOilViscosity(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
constexpr const Scalar T = 273.15 + 15.56; // [K]
// update the table for the saturated oil
saturatedOilMuTable_[regionIdx].setContainerOfTuples(samplePoints);
// calculate a table of estimated viscosities depending on pressure and gas mass
// fraction for untersaturated oil to make the other code happy
for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
Scalar p1 = std::get<0>(samplePoints[pIdx]);
Scalar p2 = p1 * 2.0;
// no pressure dependence of the viscosity
Scalar mu1 = std::get<1>(samplePoints[pIdx]);
Scalar mu2 = mu1;
Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
oilMuTable_[regionIdx].appendXPos(Rs);
oilMuTable_[regionIdx].appendSamplePoint(pIdx, p1, mu1);
oilMuTable_[regionIdx].appendSamplePoint(pIdx, p2, mu2);
}
}
/*! /*!
* \brief Finish initializing the oil phase PVT properties. * \brief Finish initializing the oil phase PVT properties.
*/ */
void initEnd() void initEnd();
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = oilMuTable_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the oil
// formation volume factor and the oil viscosity
const auto& oilMu = oilMuTable_[regionIdx];
const auto& satOilMu = saturatedOilMuTable_[regionIdx];
const auto& invOilB = inverseOilBTable_[regionIdx];
assert(oilMu.numX() == invOilB.numX());
auto& invOilBMu = inverseOilBMuTable_[regionIdx];
auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
auto& invSatOilBMu = inverseSaturatedOilBMuTable_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatOilBArray;
std::vector<Scalar> invSatOilBMuArray;
for (unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
invOilBMu.appendXPos(oilMu.xAt(rsIdx));
assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
size_t numPressures = oilMu.numY(rsIdx);
for (unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
invOilBMu.appendSamplePoint(rsIdx,
oilMu.yAt(rsIdx, pIdx),
invOilB.valueAt(rsIdx, pIdx)
/ oilMu.valueAt(rsIdx, pIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated oil is the first one
// (i.e., the one for the lowest pressure value)
satPressuresArray.push_back(oilMu.yAt(rsIdx, 0));
invSatOilBArray.push_back(invOilB.valueAt(rsIdx, 0));
invSatOilBMuArray.push_back(invSatOilBArray.back()/satOilMu.valueAt(rsIdx));
}
invSatOilB.setXYContainers(satPressuresArray, invSatOilBArray);
invSatOilBMu.setXYContainers(satPressuresArray, invSatOilBMuArray);
updateSaturationPressure_(regionIdx);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.
@ -454,36 +351,7 @@ public:
{ return vapPar2_; } { return vapPar2_; }
private: private:
void updateSaturationPressure_(unsigned regionIdx) void updateSaturationPressure_(unsigned regionIdx);
{
typedef std::pair<Scalar, Scalar> Pair;
const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
// create the function representing saturation pressure depending of the mass
// fraction in gas
size_t n = gasDissolutionFac.numSamples();
Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/Scalar(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rs = 0;
for (size_t i=0; i <= n; ++ i) {
Scalar pSat = gasDissolutionFac.xMin() + Scalar(i)*delta;
Rs = saturatedGasDissolutionFactor(regionIdx,
/*temperature=*/Scalar(1e30),
pSat);
Pair val(Rs, pSat);
pSatSamplePoints.push_back(val);
}
//Prune duplicate Rs values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
std::vector<Scalar> gasReferenceDensity_; std::vector<Scalar> gasReferenceDensity_;
std::vector<Scalar> oilReferenceDensity_; std::vector<Scalar> oilReferenceDensity_;

43
opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp
View File

@ -59,13 +59,7 @@ public:
void initFromState(const EclipseState& eclState, const Schedule&); void initFromState(const EclipseState& eclState, const Schedule&);
#endif #endif
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
solventReferenceDensity_.resize(numRegions);
inverseSolventB_.resize(numRegions);
inverseSolventBMu_.resize(numRegions);
solventMu_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference density of the solvent gas for a given PVT region * \brief Initialize the reference density of the solvent gas for a given PVT region
@ -86,42 +80,13 @@ public:
* *
* \param samplePoints A container of \f$(p_g, B_s)\f$ values * \param samplePoints A container of \f$(p_g, B_s)\f$ values
*/ */
void setSolventFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints) void setSolventFormationVolumeFactor(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
SamplingPoints tmp(samplePoints);
auto it = tmp.begin();
const auto& endIt = tmp.end();
for (; it != endIt; ++ it)
std::get<1>(*it) = 1.0/std::get<1>(*it);
inverseSolventB_[regionIdx].setContainerOfTuples(tmp);
assert(inverseSolventB_[regionIdx].monotonic());
}
/*! /*!
* \brief Finish initializing the oil phase PVT properties. * \brief Finish initializing the oil phase PVT properties.
*/ */
void initEnd() void initEnd();
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = solventMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the solvent
// formation volume factor and its viscosity
const auto& solventMu = solventMu_[regionIdx];
const auto& invSolventB = inverseSolventB_[regionIdx];
assert(solventMu.numSamples() == invSolventB.numSamples());
std::vector<Scalar> pressureValues(solventMu.numSamples());
std::vector<Scalar> invSolventBMuValues(solventMu.numSamples());
for (unsigned pIdx = 0; pIdx < solventMu.numSamples(); ++pIdx) {
pressureValues[pIdx] = invSolventB.xAt(pIdx);
invSolventBMuValues[pIdx] = invSolventB.valueAt(pIdx) * (1.0/solventMu.valueAt(pIdx));
}
inverseSolventBMu_[regionIdx].setXYContainers(pressureValues, invSolventBMuValues);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.

173
opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp
View File

@ -72,18 +72,7 @@ private:
public: public:
#endif // HAVE_ECL_INPUT #endif // HAVE_ECL_INPUT
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedOilVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference densities of all fluids for a given PVT region * \brief Initialize the reference densities of all fluids for a given PVT region
@ -91,11 +80,7 @@ public:
void setReferenceDensities(unsigned regionIdx, void setReferenceDensities(unsigned regionIdx,
Scalar rhoRefOil, Scalar rhoRefOil,
Scalar rhoRefGas, Scalar rhoRefGas,
Scalar /*rhoRefWater*/) Scalar /*rhoRefWater*/);
{
oilReferenceDensity_[regionIdx] = rhoRefOil;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
/*! /*!
* \brief Initialize the function for the oil vaporization factor \f$R_v\f$ * \brief Initialize the function for the oil vaporization factor \f$R_v\f$
@ -114,53 +99,8 @@ public:
* only depends on pressure) while the dependence on the oil mass fraction is * only depends on pressure) while the dependence on the oil mass fraction is
* guesstimated... * guesstimated...
*/ */
void setSaturatedGasFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints) void setSaturatedGasFormationVolumeFactor(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
auto& invGasB = inverseGasB_[regionIdx];
const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
constexpr const Scalar T = 273.15 + 15.56; // [K]
constexpr const Scalar RvMin = 0.0;
Scalar RvMax = RvTable.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
Scalar poMin = samplePoints.front().first;
Scalar poMax = samplePoints.back().first;
constexpr const size_t nRv = 20;
size_t nP = samplePoints.size()*2;
Scalar rhooRef = oilReferenceDensity_[regionIdx];
TabulatedOneDFunction gasFormationVolumeFactor;
gasFormationVolumeFactor.setContainerOfTuples(samplePoints);
updateSaturationPressure_(regionIdx);
// calculate a table of estimated densities depending on pressure and gas mass
// fraction. note that this assumes oil of constant compressibility. (having said
// that, if only the saturated gas densities are available, there's not much
// choice.)
for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
invGasB.appendXPos(Rv);
for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
Scalar poSat = saturationPressure(regionIdx, T, Rv);
Scalar BgSat = gasFormationVolumeFactor.eval(poSat, /*extrapolate=*/true);
Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
Scalar rhoo = rhooRef/BgSat*(1 + drhoo_dp*(pg - poSat));
Scalar Bg = rhooRef/rhoo;
invGasB.appendSamplePoint(RvIdx, pg, 1.0/Bg);
}
}
}
/*! /*!
* \brief Initialize the function for the gas formation volume factor * \brief Initialize the function for the gas formation volume factor
@ -192,85 +132,13 @@ public:
* requires the viscosity of oil-saturated gas (which only depends on pressure) while * requires the viscosity of oil-saturated gas (which only depends on pressure) while
* there is assumed to be no dependence on the gas mass fraction... * there is assumed to be no dependence on the gas mass fraction...
*/ */
void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints& samplePoints ) void setSaturatedGasViscosity(unsigned regionIdx,
{ const SamplingPoints& samplePoints);
auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
constexpr const Scalar RvMin = 0.0;
Scalar RvMax = oilVaporizationFac.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
Scalar poMin = samplePoints.front().first;
Scalar poMax = samplePoints.back().first;
constexpr const size_t nRv = 20;
size_t nP = samplePoints.size()*2;
TabulatedOneDFunction mugTable;
mugTable.setContainerOfTuples(samplePoints);
// calculate a table of estimated densities depending on pressure and gas mass
// fraction
for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
gasMu_[regionIdx].appendXPos(Rv);
for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
Scalar mug = mugTable.eval(pg, /*extrapolate=*/true);
gasMu_[regionIdx].appendSamplePoint(RvIdx, pg, mug);
}
}
}
/*! /*!
* \brief Finish initializing the gas phase PVT properties. * \brief Finish initializing the gas phase PVT properties.
*/ */
void initEnd() void initEnd();
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMu = gasMu_[regionIdx];
const auto& invGasB = inverseGasB_[regionIdx];
assert(gasMu.numX() == invGasB.numX());
auto& invGasBMu = inverseGasBMu_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
invGasBMu.appendXPos(gasMu.xAt(pIdx));
assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
size_t numRv = gasMu.numY(pIdx);
for (size_t rvIdx = 0; rvIdx < numRv; ++rvIdx)
invGasBMu.appendSamplePoint(pIdx,
gasMu.yAt(pIdx, rvIdx),
invGasB.valueAt(pIdx, rvIdx)
/ gasMu.valueAt(pIdx, rvIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rv value)
satPressuresArray.push_back(gasMu.xAt(pIdx));
invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRv - 1));
invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRv - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
updateSaturationPressure_(regionIdx);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.
@ -509,32 +377,7 @@ public:
} }
private: private:
void updateSaturationPressure_(unsigned regionIdx) void updateSaturationPressure_(unsigned regionIdx);
{
const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rv
size_t n = oilVaporizationFac.numSamples();
Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rv = 0;
for (size_t i = 0; i <= n; ++ i) {
Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
pSatSamplePoints.emplace_back(Rv, pSat);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
std::vector<Scalar> gasReferenceDensity_; std::vector<Scalar> gasReferenceDensity_;
std::vector<Scalar> oilReferenceDensity_; std::vector<Scalar> oilReferenceDensity_;

140
opm/material/fluidsystems/blackoilpvt/WetHumidGasPvt.hpp
View File

@ -77,24 +77,7 @@ private:
public: public:
#endif // HAVE_ECL_INPUT #endif // HAVE_ECL_INPUT
void setNumRegions(size_t numRegions) void setNumRegions(size_t numRegions);
{
waterReferenceDensity_.resize(numRegions);
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasBRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMuRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMuRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMuRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
gasMuRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedWaterVaporizationFactorTable_.resize(numRegions);
saturatedWaterVaporizationSaltFactorTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedOilVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
/*! /*!
* \brief Initialize the reference densities of all fluids for a given PVT region * \brief Initialize the reference densities of all fluids for a given PVT region
@ -102,12 +85,7 @@ public:
void setReferenceDensities(unsigned regionIdx, void setReferenceDensities(unsigned regionIdx,
Scalar rhoRefOil, Scalar rhoRefOil,
Scalar rhoRefGas, Scalar rhoRefGas,
Scalar rhoRefWater) Scalar rhoRefWater);
{
waterReferenceDensity_[regionIdx] = rhoRefWater;
oilReferenceDensity_[regionIdx] = rhoRefOil;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
/*! /*!
* \brief Initialize the function for the water vaporization factor \f$R_v\f$ * \brief Initialize the function for the water vaporization factor \f$R_v\f$
@ -129,93 +107,7 @@ public:
/*! /*!
* \brief Finish initializing the gas phase PVT properties. * \brief Finish initializing the gas phase PVT properties.
*/ */
void initEnd() void initEnd();
{
//PVTGW
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMuRvSat_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMuRvSat = gasMuRvSat_[regionIdx];
const auto& invGasBRvSat = inverseGasBRvSat_[regionIdx];
assert(gasMuRvSat.numX() == invGasBRvSat.numX());
auto& invGasBMuRvSat = inverseGasBMuRvSat_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMuRvSat.numX(); ++pIdx) {
invGasBMuRvSat.appendXPos(gasMuRvSat.xAt(pIdx));
assert(gasMuRvSat.numY(pIdx) == invGasBRvSat.numY(pIdx));
size_t numRw = gasMuRvSat.numY(pIdx);
for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
invGasBMuRvSat.appendSamplePoint(pIdx,
gasMuRvSat.yAt(pIdx, RwIdx),
invGasBRvSat.valueAt(pIdx, RwIdx)
/ gasMuRvSat.valueAt(pIdx, RwIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rw value)
satPressuresArray.push_back(gasMuRvSat.xAt(pIdx));
invSatGasBArray.push_back(invGasBRvSat.valueAt(pIdx, numRw - 1));
invSatGasBMuArray.push_back(invGasBMuRvSat.valueAt(pIdx, numRw - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
}
//PVTG
// calculate the final 2D functions which are used for interpolation.
//size_t numRegions = gasMuRvwSat_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMuRvwSat = gasMuRvwSat_[regionIdx];
const auto& invGasBRvwSat = inverseGasBRvwSat_[regionIdx];
assert(gasMuRvwSat.numX() == invGasBRvwSat.numX());
auto& invGasBMuRvwSat = inverseGasBMuRvwSat_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMuRvwSat.numX(); ++pIdx) {
invGasBMuRvwSat.appendXPos(gasMuRvwSat.xAt(pIdx));
assert(gasMuRvwSat.numY(pIdx) == invGasBRvwSat.numY(pIdx));
size_t numRw = gasMuRvwSat.numY(pIdx);
for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
invGasBMuRvwSat.appendSamplePoint(pIdx,
gasMuRvwSat.yAt(pIdx, RwIdx),
invGasBRvwSat.valueAt(pIdx, RwIdx)
/ gasMuRvwSat.valueAt(pIdx, RwIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rw value)
satPressuresArray.push_back(gasMuRvwSat.xAt(pIdx));
invSatGasBArray.push_back(invGasBRvwSat.valueAt(pIdx, numRw - 1));
invSatGasBMuArray.push_back(invGasBMuRvwSat.valueAt(pIdx, numRw - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
updateSaturationPressure_(regionIdx);
}
}
/*! /*!
* \brief Return the number of PVT regions which are considered by this PVT-object. * \brief Return the number of PVT regions which are considered by this PVT-object.
@ -497,31 +389,7 @@ public:
} }
private: private:
void updateSaturationPressure_(unsigned regionIdx) void updateSaturationPressure_(unsigned regionIdx);
{
const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rv
size_t n = oilVaporizationFac.numSamples();
Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rv = 0;
for (size_t i = 0; i <= n; ++ i) {
Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
pSatSamplePoints.emplace_back(Rv, pSat);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
std::vector<Scalar> gasReferenceDensity_; std::vector<Scalar> gasReferenceDensity_;
std::vector<Scalar> oilReferenceDensity_; std::vector<Scalar> oilReferenceDensity_;

41
src/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.cpp
View File

@ -26,14 +26,17 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PvdoTable.hpp> #include <opm/input/eclipse/EclipseState/Tables/PvdoTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void DeadOilPvt<Scalar>:: void DeadOilPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule&) initFromState(const EclipseState& eclState, const Schedule&)
@ -74,6 +77,44 @@ initFromState(const EclipseState& eclState, const Schedule&)
initEnd(); initEnd();
} }
#endif
template<class Scalar>
void DeadOilPvt<Scalar>::setNumRegions(size_t numRegions)
{
oilReferenceDensity_.resize(numRegions);
inverseOilB_.resize(numRegions);
inverseOilBMu_.resize(numRegions);
oilMu_.resize(numRegions);
}
template<class Scalar>
void DeadOilPvt<Scalar>::initEnd()
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = oilMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the oil
// formation volume factor and the oil viscosity
const auto& oilMu = oilMu_[regionIdx];
const auto& invOilB = inverseOilB_[regionIdx];
assert(oilMu.numSamples() == invOilB.numSamples());
std::vector<Scalar> invBMuColumn;
std::vector<Scalar> pressureColumn;
invBMuColumn.resize(oilMu.numSamples());
pressureColumn.resize(oilMu.numSamples());
for (unsigned pIdx = 0; pIdx < oilMu.numSamples(); ++pIdx) {
pressureColumn[pIdx] = invOilB.xAt(pIdx);
invBMuColumn[pIdx] = invOilB.valueAt(pIdx) / oilMu.valueAt(pIdx);
}
inverseOilBMu_[regionIdx].setXYArrays(pressureColumn.size(),
pressureColumn,
invBMuColumn);
}
}
template class DeadOilPvt<double>; template class DeadOilPvt<double>;
template class DeadOilPvt<float>; template class DeadOilPvt<float>;

51
src/opm/material/fluidsystems/blackoilpvt/DryGasPvt.cpp
View File

@ -26,14 +26,17 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PvdgTable.hpp> #include <opm/input/eclipse/EclipseState/Tables/PvdgTable.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void DryGasPvt<Scalar>:: void DryGasPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule&) initFromState(const EclipseState& eclState, const Schedule&)
@ -84,6 +87,54 @@ initFromState(const EclipseState& eclState, const Schedule&)
initEnd(); initEnd();
} }
#endif
template<class Scalar>
void DryGasPvt<Scalar>::setNumRegions(size_t numRegions)
{
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions);
inverseGasBMu_.resize(numRegions);
gasMu_.resize(numRegions);
}
template<class Scalar>
void DryGasPvt<Scalar>::
setGasFormationVolumeFactor(unsigned regionIdx,
const SamplingPoints& samplePoints)
{
SamplingPoints tmp(samplePoints);
auto it = tmp.begin();
const auto& endIt = tmp.end();
for (; it != endIt; ++ it)
std::get<1>(*it) = 1.0/std::get<1>(*it);
inverseGasB_[regionIdx].setContainerOfTuples(tmp);
assert(inverseGasB_[regionIdx].monotonic());
}
template<class Scalar>
void DryGasPvt<Scalar>::initEnd()
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMu = gasMu_[regionIdx];
const auto& invGasB = inverseGasB_[regionIdx];
assert(gasMu.numSamples() == invGasB.numSamples());
std::vector<Scalar> pressureValues(gasMu.numSamples());
std::vector<Scalar> invGasBMuValues(gasMu.numSamples());
for (unsigned pIdx = 0; pIdx < gasMu.numSamples(); ++pIdx) {
pressureValues[pIdx] = invGasB.xAt(pIdx);
invGasBMuValues[pIdx] = invGasB.valueAt(pIdx) * (1.0/gasMu.valueAt(pIdx));
}
inverseGasBMu_[regionIdx].setXYContainers(pressureValues, invGasBMuValues);
}
}
template class DryGasPvt<double>; template class DryGasPvt<double>;
template class DryGasPvt<float>; template class DryGasPvt<float>;

140
src/opm/material/fluidsystems/blackoilpvt/DryHumidGasPvt.cpp
View File

@ -26,13 +26,16 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void DryHumidGasPvt<Scalar>:: void DryHumidGasPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule&) initFromState(const EclipseState& eclState, const Schedule&)
@ -233,6 +236,143 @@ extendPvtgwTable_(unsigned regionIdx,
gasMu.appendSamplePoint(xIdx, newRw, newMug); gasMu.appendSamplePoint(xIdx, newRw, newMug);
} }
} }
#endif
template<class Scalar>
void DryHumidGasPvt<Scalar>::setNumRegions(size_t numRegions)
{
waterReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedWaterVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
template<class Scalar>
void DryHumidGasPvt<Scalar>::
setReferenceDensities(unsigned regionIdx,
Scalar /*rhoRefOil*/,
Scalar rhoRefGas,
Scalar rhoRefWater)
{
waterReferenceDensity_[regionIdx] = rhoRefWater;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
template<class Scalar>
void DryHumidGasPvt<Scalar>::
setSaturatedGasViscosity(unsigned regionIdx,
const SamplingPoints& samplePoints )
{
auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
constexpr const Scalar RwMin = 0.0;
Scalar RwMax = waterVaporizationFac.eval(saturatedWaterVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
Scalar poMin = samplePoints.front().first;
Scalar poMax = samplePoints.back().first;
constexpr const size_t nRw = 20;
size_t nP = samplePoints.size()*2;
TabulatedOneDFunction mugTable;
mugTable.setContainerOfTuples(samplePoints);
// calculate a table of estimated densities depending on pressure and gas mass
// fraction
for (size_t RwIdx = 0; RwIdx < nRw; ++RwIdx) {
Scalar Rw = RwMin + (RwMax - RwMin)*RwIdx/nRw;
gasMu_[regionIdx].appendXPos(Rw);
for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
Scalar mug = mugTable.eval(pg, /*extrapolate=*/true);
gasMu_[regionIdx].appendSamplePoint(RwIdx, pg, mug);
}
}
}
template<class Scalar>
void DryHumidGasPvt<Scalar>::initEnd()
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMu = gasMu_[regionIdx];
const auto& invGasB = inverseGasB_[regionIdx];
assert(gasMu.numX() == invGasB.numX());
auto& invGasBMu = inverseGasBMu_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
invGasBMu.appendXPos(gasMu.xAt(pIdx));
assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
size_t numRw = gasMu.numY(pIdx);
for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
invGasBMu.appendSamplePoint(pIdx,
gasMu.yAt(pIdx, RwIdx),
invGasB.valueAt(pIdx, RwIdx)
/ gasMu.valueAt(pIdx, RwIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rw value)
satPressuresArray.push_back(gasMu.xAt(pIdx));
invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRw - 1));
invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRw - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
updateSaturationPressure_(regionIdx);
}
}
template<class Scalar>
void DryHumidGasPvt<Scalar>::
updateSaturationPressure_(unsigned regionIdx)
{
const auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rw
size_t n = waterVaporizationFac.numSamples();
const Scalar delta = (waterVaporizationFac.xMax() -
waterVaporizationFac.xMin()) / Scalar(n + 1);
SamplingPoints pSatSamplePoints;
for (size_t i = 0; i <= n; ++ i) {
const Scalar pSat = waterVaporizationFac.xMin() + i*delta;
const Scalar Rw = saturatedWaterVaporizationFactor(regionIdx,
Scalar(1e30),
pSat);
pSatSamplePoints.emplace_back(Rw, pSat);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
template class DryHumidGasPvt<double>; template class DryHumidGasPvt<double>;
template class DryHumidGasPvt<float>; template class DryHumidGasPvt<float>;

161
src/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.cpp
View File

@ -26,15 +26,18 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp> #include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp> #include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void LiveOilPvt<Scalar>:: void LiveOilPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule& schedule) initFromState(const EclipseState& eclState, const Schedule& schedule)
@ -210,6 +213,164 @@ extendPvtoTable_(unsigned regionIdx,
oilMu.appendSamplePoint(xIdx, newPo, newMuo); oilMu.appendSamplePoint(xIdx, newPo, newMuo);
} }
} }
#endif
template<class Scalar>
void LiveOilPvt<Scalar>::setNumRegions(size_t numRegions)
{
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseOilBTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
inverseOilBMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
inverseSaturatedOilBTable_.resize(numRegions);
inverseSaturatedOilBMuTable_.resize(numRegions);
oilMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
saturatedOilMuTable_.resize(numRegions);
saturatedGasDissolutionFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
template<class Scalar>
void LiveOilPvt<Scalar>::
setReferenceDensities(unsigned regionIdx,
Scalar rhoRefOil,
Scalar rhoRefGas,
Scalar)
{
oilReferenceDensity_[regionIdx] = rhoRefOil;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
template<class Scalar>
void LiveOilPvt<Scalar>::
setSaturatedOilFormationVolumeFactor(unsigned regionIdx,
const SamplingPoints& samplePoints)
{
constexpr const Scalar T = 273.15 + 15.56; // [K]
auto& invOilB = inverseOilBTable_[regionIdx];
updateSaturationPressure_(regionIdx);
// calculate a table of estimated densities of undersatured gas
for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
Scalar p1 = std::get<0>(samplePoints[pIdx]);
Scalar p2 = p1 * 2.0;
Scalar Bo1 = std::get<1>(samplePoints[pIdx]);
Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
Scalar Bo2 = Bo1/(1.0 + (p2 - p1)*drhoo_dp);
Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
invOilB.appendXPos(Rs);
invOilB.appendSamplePoint(pIdx, p1, 1.0/Bo1);
invOilB.appendSamplePoint(pIdx, p2, 1.0/Bo2);
}
}
template<class Scalar>
void LiveOilPvt<Scalar>::
setSaturatedOilViscosity(unsigned regionIdx,
const SamplingPoints& samplePoints)
{
constexpr const Scalar T = 273.15 + 15.56; // [K]
// update the table for the saturated oil
saturatedOilMuTable_[regionIdx].setContainerOfTuples(samplePoints);
// calculate a table of estimated viscosities depending on pressure and gas mass
// fraction for untersaturated oil to make the other code happy
for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
Scalar p1 = std::get<0>(samplePoints[pIdx]);
Scalar p2 = p1 * 2.0;
// no pressure dependence of the viscosity
Scalar mu1 = std::get<1>(samplePoints[pIdx]);
Scalar mu2 = mu1;
Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
oilMuTable_[regionIdx].appendXPos(Rs);
oilMuTable_[regionIdx].appendSamplePoint(pIdx, p1, mu1);
oilMuTable_[regionIdx].appendSamplePoint(pIdx, p2, mu2);
}
}
template<class Scalar>
void LiveOilPvt<Scalar>::initEnd()
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = oilMuTable_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the oil
// formation volume factor and the oil viscosity
const auto& oilMu = oilMuTable_[regionIdx];
const auto& satOilMu = saturatedOilMuTable_[regionIdx];
const auto& invOilB = inverseOilBTable_[regionIdx];
assert(oilMu.numX() == invOilB.numX());
auto& invOilBMu = inverseOilBMuTable_[regionIdx];
auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
auto& invSatOilBMu = inverseSaturatedOilBMuTable_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatOilBArray;
std::vector<Scalar> invSatOilBMuArray;
for (unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
invOilBMu.appendXPos(oilMu.xAt(rsIdx));
assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
size_t numPressures = oilMu.numY(rsIdx);
for (unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
invOilBMu.appendSamplePoint(rsIdx,
oilMu.yAt(rsIdx, pIdx),
invOilB.valueAt(rsIdx, pIdx)
/ oilMu.valueAt(rsIdx, pIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated oil is the first one
// (i.e., the one for the lowest pressure value)
satPressuresArray.push_back(oilMu.yAt(rsIdx, 0));
invSatOilBArray.push_back(invOilB.valueAt(rsIdx, 0));
invSatOilBMuArray.push_back(invSatOilBArray.back()/satOilMu.valueAt(rsIdx));
}
invSatOilB.setXYContainers(satPressuresArray, invSatOilBArray);
invSatOilBMu.setXYContainers(satPressuresArray, invSatOilBMuArray);
updateSaturationPressure_(regionIdx);
}
}
template<class Scalar>
void LiveOilPvt<Scalar>::updateSaturationPressure_(unsigned regionIdx)
{
const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
// create the function representing saturation pressure depending of the mass
// fraction in gas
size_t n = gasDissolutionFac.numSamples();
const Scalar delta = (gasDissolutionFac.xMax() -
gasDissolutionFac.xMin()) / Scalar(n + 1);
SamplingPoints pSatSamplePoints;
for (size_t i = 0; i <= n; ++ i) {
const Scalar pSat = gasDissolutionFac.xMin() + i*delta;
const Scalar Rs = saturatedGasDissolutionFactor(regionIdx,
Scalar(1e30),
pSat);
pSatSamplePoints.emplace_back(Rs, pSat);
}
//Prune duplicate Rs values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
template class LiveOilPvt<double>; template class LiveOilPvt<double>;
template class LiveOilPvt<float>; template class LiveOilPvt<float>;

51
src/opm/material/fluidsystems/blackoilpvt/SolventPvt.cpp
View File

@ -26,14 +26,17 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PvdsTable.hpp> #include <opm/input/eclipse/EclipseState/Tables/PvdsTable.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void SolventPvt<Scalar>:: void SolventPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule&) initFromState(const EclipseState& eclState, const Schedule&)
@ -72,6 +75,54 @@ initFromState(const EclipseState& eclState, const Schedule&)
initEnd(); initEnd();
} }
#endif
template<class Scalar>
void SolventPvt<Scalar>::setNumRegions(size_t numRegions)
{
solventReferenceDensity_.resize(numRegions);
inverseSolventB_.resize(numRegions);
inverseSolventBMu_.resize(numRegions);
solventMu_.resize(numRegions);
}
template<class Scalar>
void SolventPvt<Scalar>::
setSolventFormationVolumeFactor(unsigned regionIdx,
const SamplingPoints& samplePoints)
{
SamplingPoints tmp(samplePoints);
auto it = tmp.begin();
const auto& endIt = tmp.end();
for (; it != endIt; ++ it)
std::get<1>(*it) = 1.0/std::get<1>(*it);
inverseSolventB_[regionIdx].setContainerOfTuples(tmp);
assert(inverseSolventB_[regionIdx].monotonic());
}
template<class Scalar>
void SolventPvt<Scalar>::initEnd()
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = solventMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the solvent
// formation volume factor and its viscosity
const auto& solventMu = solventMu_[regionIdx];
const auto& invSolventB = inverseSolventB_[regionIdx];
assert(solventMu.numSamples() == invSolventB.numSamples());
std::vector<Scalar> pressureValues(solventMu.numSamples());
std::vector<Scalar> invSolventBMuValues(solventMu.numSamples());
for (unsigned pIdx = 0; pIdx < solventMu.numSamples(); ++pIdx) {
pressureValues[pIdx] = invSolventB.xAt(pIdx);
invSolventBMuValues[pIdx] = invSolventB.valueAt(pIdx) * (1.0/solventMu.valueAt(pIdx));
}
inverseSolventBMu_[regionIdx].setXYContainers(pressureValues, invSolventBMuValues);
}
}
template class SolventPvt<double>; template class SolventPvt<double>;
template class SolventPvt<float>; template class SolventPvt<float>;

190
src/opm/material/fluidsystems/blackoilpvt/WetGasPvt.cpp
View File

@ -26,15 +26,18 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp> #include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp> #include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void WetGasPvt<Scalar>:: void WetGasPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule& schedule) initFromState(const EclipseState& eclState, const Schedule& schedule)
@ -211,6 +214,193 @@ extendPvtgTable_(unsigned regionIdx,
gasMu.appendSamplePoint(xIdx, newRv, newMug); gasMu.appendSamplePoint(xIdx, newRv, newMug);
} }
} }
#endif
template<class Scalar>
void WetGasPvt<Scalar>::setNumRegions(size_t numRegions)
{
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedOilVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
template<class Scalar>
void WetGasPvt<Scalar>::
setReferenceDensities(unsigned regionIdx,
Scalar rhoRefOil,
Scalar rhoRefGas,
Scalar)
{
oilReferenceDensity_[regionIdx] = rhoRefOil;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
template<class Scalar>
void WetGasPvt<Scalar>::
setSaturatedGasFormationVolumeFactor(unsigned regionIdx,
const SamplingPoints& samplePoints)
{
auto& invGasB = inverseGasB_[regionIdx];
const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
constexpr const Scalar T = 273.15 + 15.56; // [K]
constexpr const Scalar RvMin = 0.0;
Scalar RvMax = RvTable.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
Scalar poMin = samplePoints.front().first;
Scalar poMax = samplePoints.back().first;
constexpr const size_t nRv = 20;
size_t nP = samplePoints.size()*2;
Scalar rhooRef = oilReferenceDensity_[regionIdx];
TabulatedOneDFunction gasFormationVolumeFactor;
gasFormationVolumeFactor.setContainerOfTuples(samplePoints);
updateSaturationPressure_(regionIdx);
// calculate a table of estimated densities depending on pressure and gas mass
// fraction. note that this assumes oil of constant compressibility. (having said
// that, if only the saturated gas densities are available, there's not much
// choice.)
for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
invGasB.appendXPos(Rv);
for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
Scalar poSat = saturationPressure(regionIdx, T, Rv);
Scalar BgSat = gasFormationVolumeFactor.eval(poSat, /*extrapolate=*/true);
Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
Scalar rhoo = rhooRef/BgSat*(1 + drhoo_dp*(pg - poSat));
Scalar Bg = rhooRef/rhoo;
invGasB.appendSamplePoint(RvIdx, pg, 1.0/Bg);
}
}
}
template<class Scalar>
void WetGasPvt<Scalar>::
setSaturatedGasViscosity(unsigned regionIdx,
const SamplingPoints& samplePoints)
{
auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
constexpr const Scalar RvMin = 0.0;
Scalar RvMax = oilVaporizationFac.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
Scalar poMin = samplePoints.front().first;
Scalar poMax = samplePoints.back().first;
constexpr const size_t nRv = 20;
size_t nP = samplePoints.size()*2;
TabulatedOneDFunction mugTable;
mugTable.setContainerOfTuples(samplePoints);
// calculate a table of estimated densities depending on pressure and gas mass
// fraction
for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
gasMu_[regionIdx].appendXPos(Rv);
for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
Scalar mug = mugTable.eval(pg, /*extrapolate=*/true);
gasMu_[regionIdx].appendSamplePoint(RvIdx, pg, mug);
}
}
}
template<class Scalar>
void WetGasPvt<Scalar>::initEnd()
{
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMu_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMu = gasMu_[regionIdx];
const auto& invGasB = inverseGasB_[regionIdx];
assert(gasMu.numX() == invGasB.numX());
auto& invGasBMu = inverseGasBMu_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
invGasBMu.appendXPos(gasMu.xAt(pIdx));
assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
size_t numRv = gasMu.numY(pIdx);
for (size_t rvIdx = 0; rvIdx < numRv; ++rvIdx)
invGasBMu.appendSamplePoint(pIdx,
gasMu.yAt(pIdx, rvIdx),
invGasB.valueAt(pIdx, rvIdx)
/ gasMu.valueAt(pIdx, rvIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rv value)
satPressuresArray.push_back(gasMu.xAt(pIdx));
invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRv - 1));
invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRv - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
updateSaturationPressure_(regionIdx);
}
}
template<class Scalar>
void WetGasPvt<Scalar>::
updateSaturationPressure_(unsigned regionIdx)
{
const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rv
size_t n = oilVaporizationFac.numSamples();
Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rv = 0;
for (size_t i = 0; i <= n; ++ i) {
Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
pSatSamplePoints.emplace_back(Rv, pSat);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
template class WetGasPvt<double>; template class WetGasPvt<double>;
template class WetGasPvt<float>; template class WetGasPvt<float>;

153
src/opm/material/fluidsystems/blackoilpvt/WetHumidGasPvt.cpp
View File

@ -26,14 +26,17 @@
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp> #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp> #include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp> #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar> template<class Scalar>
void WetHumidGasPvt<Scalar>:: void WetHumidGasPvt<Scalar>::
initFromState(const EclipseState& eclState, const Schedule& schedule) initFromState(const EclipseState& eclState, const Schedule& schedule)
@ -392,6 +395,156 @@ extendPvtgTable_(unsigned regionIdx,
gasMuRvwSat.appendSamplePoint(xIdx, newRv, newMug); gasMuRvwSat.appendSamplePoint(xIdx, newRv, newMug);
} }
} }
#endif
template<class Scalar>
void WetHumidGasPvt<Scalar>::setNumRegions(size_t numRegions)
{
waterReferenceDensity_.resize(numRegions);
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasBRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMuRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMuRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMuRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
gasMuRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedWaterVaporizationFactorTable_.resize(numRegions);
saturatedWaterVaporizationSaltFactorTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedOilVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}
template<class Scalar>
void WetHumidGasPvt<Scalar>::
setReferenceDensities(unsigned regionIdx,
Scalar rhoRefOil,
Scalar rhoRefGas,
Scalar rhoRefWater)
{
waterReferenceDensity_[regionIdx] = rhoRefWater;
oilReferenceDensity_[regionIdx] = rhoRefOil;
gasReferenceDensity_[regionIdx] = rhoRefGas;
}
template<class Scalar>
void WetHumidGasPvt<Scalar>::initEnd()
{
//PVTGW
// calculate the final 2D functions which are used for interpolation.
size_t numRegions = gasMuRvSat_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMuRvSat = gasMuRvSat_[regionIdx];
const auto& invGasBRvSat = inverseGasBRvSat_[regionIdx];
assert(gasMuRvSat.numX() == invGasBRvSat.numX());
auto& invGasBMuRvSat = inverseGasBMuRvSat_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMuRvSat.numX(); ++pIdx) {
invGasBMuRvSat.appendXPos(gasMuRvSat.xAt(pIdx));
assert(gasMuRvSat.numY(pIdx) == invGasBRvSat.numY(pIdx));
size_t numRw = gasMuRvSat.numY(pIdx);
for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
invGasBMuRvSat.appendSamplePoint(pIdx,
gasMuRvSat.yAt(pIdx, RwIdx),
invGasBRvSat.valueAt(pIdx, RwIdx)
/ gasMuRvSat.valueAt(pIdx, RwIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rw value)
satPressuresArray.push_back(gasMuRvSat.xAt(pIdx));
invSatGasBArray.push_back(invGasBRvSat.valueAt(pIdx, numRw - 1));
invSatGasBMuArray.push_back(invGasBMuRvSat.valueAt(pIdx, numRw - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
}
//PVTG
// calculate the final 2D functions which are used for interpolation.
//size_t numRegions = gasMuRvwSat_.size();
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
// calculate the table which stores the inverse of the product of the gas
// formation volume factor and the gas viscosity
const auto& gasMuRvwSat = gasMuRvwSat_[regionIdx];
const auto& invGasBRvwSat = inverseGasBRvwSat_[regionIdx];
assert(gasMuRvwSat.numX() == invGasBRvwSat.numX());
auto& invGasBMuRvwSat = inverseGasBMuRvwSat_[regionIdx];
auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
std::vector<Scalar> satPressuresArray;
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
for (size_t pIdx = 0; pIdx < gasMuRvwSat.numX(); ++pIdx) {
invGasBMuRvwSat.appendXPos(gasMuRvwSat.xAt(pIdx));
assert(gasMuRvwSat.numY(pIdx) == invGasBRvwSat.numY(pIdx));
size_t numRw = gasMuRvwSat.numY(pIdx);
for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
invGasBMuRvwSat.appendSamplePoint(pIdx,
gasMuRvwSat.yAt(pIdx, RwIdx),
invGasBRvwSat.valueAt(pIdx, RwIdx)
/ gasMuRvwSat.valueAt(pIdx, RwIdx));
// the sampling points in UniformXTabulated2DFunction are always sorted
// in ascending order. Thus, the value for saturated gas is the last one
// (i.e., the one with the largest Rw value)
satPressuresArray.push_back(gasMuRvwSat.xAt(pIdx));
invSatGasBArray.push_back(invGasBRvwSat.valueAt(pIdx, numRw - 1));
invSatGasBMuArray.push_back(invGasBMuRvwSat.valueAt(pIdx, numRw - 1));
}
invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
updateSaturationPressure_(regionIdx);
}
}
template<class Scalar>
void WetHumidGasPvt<Scalar>::
updateSaturationPressure_(unsigned regionIdx)
{
const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rv
size_t n = oilVaporizationFac.numSamples();
Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rv = 0;
for (size_t i = 0; i <= n; ++ i) {
Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
pSatSamplePoints.emplace_back(Rv, pSat);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}
template class WetHumidGasPvt<double>; template class WetHumidGasPvt<double>;
template class WetHumidGasPvt<float>; template class WetHumidGasPvt<float>;