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

Always Store Three Saturation End-Points

Browse Source 
The two-point scalings will use [0] and [2] instead of [0] and [1].
This is in preparation of introducing 'KRWR' scaling which does not
necessarily require the alternative scaling method ('SCALECRS') for
the horizontal direction.
This commit is contained in:
Bård Skaflestad 2020-12-20 11:51:14 +01:00
parent bc70036de5
commit bb191a5381
2 changed files with 78 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp
View File

@ -198,7 +198,6 @@ struct EclEpsScalingPointsInfo
update(maxPcgo, epsProperties.pcg(activeIndex));
update(maxKrw, epsProperties.krw(activeIndex));
update(maxKrg, epsProperties.krg(activeIndex));
// quite likely that's wrong!
update(maxKrow, epsProperties.kro(activeIndex));
update(maxKrog, epsProperties.kro(activeIndex));
@ -298,32 +297,20 @@ public:
if (epsSystemType == EclOilWaterSystem) {
// saturation scaling for capillary pressure
saturationPcPoints_[0] = epsInfo.Swl;
saturationPcPoints_[1] = epsInfo.Swu;
saturationPcPoints_[2] = saturationPcPoints_[1] = epsInfo.Swu;
// krw saturation scaling endpoints
if (config.enableThreePointKrSatScaling()) {
saturationKrwPoints_[0] = epsInfo.Swcr;
saturationKrwPoints_[1] = 1.0 - epsInfo.Sowcr - epsInfo.Sgl;
saturationKrwPoints_[2] = epsInfo.Swu;
}
else {
saturationKrwPoints_[0] = epsInfo.Swcr;
saturationKrwPoints_[1] = epsInfo.Swu;
}
saturationKrwPoints_[0] = epsInfo.Swcr;
saturationKrwPoints_[1] = 1.0 - epsInfo.Sowcr - epsInfo.Sgl;
saturationKrwPoints_[2] = epsInfo.Swu;
// krn saturation scaling endpoints (with the non-wetting phase being oil).
// because opm-material specifies non-wetting phase relperms in terms of the
// wetting phase saturations, the code here uses 1 minus the values specified
// by the Eclipse TD and the order of the scaling points is reversed
if (config.enableThreePointKrSatScaling()) {
saturationKrnPoints_[2] = 1.0 - epsInfo.Sowcr;
saturationKrnPoints_[1] = epsInfo.Swcr + epsInfo.Sgl;
saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
}
else {
saturationKrnPoints_[1] = 1 - epsInfo.Sowcr;
saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
}
saturationKrnPoints_[2] = 1.0 - epsInfo.Sowcr;
saturationKrnPoints_[1] = epsInfo.Swcr + epsInfo.Sgl;
saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
if (config.enableLeverettScaling())
maxPcnwOrLeverettFactor_ = epsInfo.pcowLeverettFactor;
@ -337,32 +324,20 @@ public:
// saturation scaling for capillary pressure
saturationPcPoints_[0] = 1.0 - epsInfo.Sgu;
saturationPcPoints_[1] = 1.0 - epsInfo.Sgl;
saturationPcPoints_[2] = saturationPcPoints_[1] = 1.0 - epsInfo.Sgl;
// krw saturation scaling endpoints
if (config.enableThreePointKrSatScaling()) {
saturationKrwPoints_[0] = epsInfo.Sogcr;
saturationKrwPoints_[1] = 1 - epsInfo.Sgcr - epsInfo.Swl;
saturationKrwPoints_[2] = 1 - epsInfo.Swl - epsInfo.Sgl;
}
else {
saturationKrwPoints_[0] = epsInfo.Sogcr;
saturationKrwPoints_[1] = 1 - epsInfo.Swl - epsInfo.Sgl;
}
saturationKrwPoints_[0] = epsInfo.Sogcr;
saturationKrwPoints_[1] = 1 - epsInfo.Sgcr - epsInfo.Swl;
saturationKrwPoints_[2] = 1 - epsInfo.Swl - epsInfo.Sgl;
// krn saturation scaling endpoints (with the non-wetting phase being gas).
// because opm-material specifies non-wetting phase relperms in terms of the
// wetting phase saturations, the code here uses 1 minus the values specified
// by the Eclipse TD and the order of the scaling points is reversed
if (config.enableThreePointKrSatScaling()) {
saturationKrnPoints_[2] = 1.0 - epsInfo.Sgcr;
saturationKrnPoints_[1] = epsInfo.Sogcr + epsInfo.Swl;
saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
}
else {
saturationKrnPoints_[1] = 1.0 - epsInfo.Sgcr;
saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
}
saturationKrnPoints_[2] = 1.0 - epsInfo.Sgcr;
saturationKrnPoints_[1] = epsInfo.Sogcr + epsInfo.Swl;
saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
if (config.enableLeverettScaling())
maxPcnwOrLeverettFactor_ = epsInfo.pcgoLeverettFactor;
@ -383,7 +358,7 @@ public:
/*!
* \brief Returns the points used for capillary pressure saturation scaling
*/
const std::array<Scalar, 2>& saturationPcPoints() const
const std::array<Scalar, 3>& saturationPcPoints() const
{ return saturationPcPoints_; }
/*!
@ -466,17 +441,17 @@ public:
}
private:
// The the points used for the "y-axis" scaling of capillary pressure
// Points used for vertical scaling of capillary pressure
Scalar maxPcnwOrLeverettFactor_;
// The the points used for the "y-axis" scaling of wetting phase relative permability
// Maximum wetting phase relative permability value.
Scalar maxKrw_;
// The the points used for the "y-axis" scaling of non-wetting phase relative permability
// Maximum non-wetting phase relative permability value
Scalar maxKrn_;
// The the points used for saturation ("x-axis") scaling of capillary pressure
std::array<Scalar, 2> saturationPcPoints_;
std::array<Scalar, 3> saturationPcPoints_;
// The the points used for saturation ("x-axis") scaling of wetting phase relative permeability
std::array<Scalar, 3> saturationKrwPoints_;

97
opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp
View File

@ -31,6 +31,10 @@
#include <opm/material/fluidstates/SaturationOverlayFluidState.hpp>
#include <algorithm>
#include <cstddef>
#include <type_traits>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
@ -376,7 +380,7 @@ private:
return
unscaledSats[0]
+
(scaledSat - scaledSats[0])*((unscaledSats[1] - unscaledSats[0])/(scaledSats[1] - scaledSats[0]));
(scaledSat - scaledSats[0])*((unscaledSats[2] - unscaledSats[0])/(scaledSats[2] - scaledSats[0]));
}
template <class Evaluation, class PointsContainer>
@ -387,7 +391,7 @@ private:
return
scaledSats[0]
+
(unscaledSat - unscaledSats[0])*((scaledSats[1] - scaledSats[0])/(unscaledSats[1] - unscaledSats[0]));
(unscaledSat - unscaledSats[0])*((scaledSats[2] - scaledSats[0])/(unscaledSats[2] - unscaledSats[0]));
}
template <class Evaluation, class PointsContainer>
@ -395,28 +399,38 @@ private:
const PointsContainer& unscaledSats,
const PointsContainer& scaledSats)
{
if (unscaledSats[1] >= unscaledSats[2])
return scaledToUnscaledSatTwoPoint_(scaledSat, unscaledSats, scaledSats);
using UnscaledSat = std::remove_cv_t<std::remove_reference_t<decltype(unscaledSats[0])>>;
if (scaledSat < scaledSats[1]) {
Scalar delta = scaledSats[1] - scaledSats[0];
if (delta <= 1e-20)
delta = 1.0; // prevent division by zero for (possibly) incorrect input data
auto map = [&scaledSat, &unscaledSats, &scaledSats](const std::size_t i)
{
const auto distance = (scaledSat - scaledSats[i])
/ (scaledSats[i + 1] - scaledSats[i]);
return
unscaledSats[0]
+
(scaledSat - scaledSats[0])*((unscaledSats[1] - unscaledSats[0])/delta);
const auto displacement =
std::max(unscaledSats[i + 1] - unscaledSats[i], UnscaledSat{ 0 });
return std::min(unscaledSats[i] + distance*displacement,
Evaluation { unscaledSats[i + 1] });
};
if (! (scaledSat > scaledSats[0])) {
// s <= sL
return unscaledSats[0];
}
else if (scaledSat < std::min(scaledSats[1], scaledSats[2])) {
// Scaled saturation in interval [sL, sR).
// Map to tabulated saturation in [unscaledSats[0], unscaledSats[1]).
return map(0);
}
else if (scaledSat < scaledSats[2]) {
// Scaled saturation in interval [sR, sU); sR guaranteed to be less
// than sU from previous condition. Map to tabulated saturation in
// [unscaledSats[1], unscaledSats[2]).
return map(1);
}
else {
Scalar delta = scaledSats[2] - scaledSats[1];
if (delta <= 1e-20)
delta = 1.0; // prevent division by zero for (possibly) incorrect input data
return
unscaledSats[1]
+
(scaledSat - scaledSats[1])*((unscaledSats[2] - unscaledSats[1])/delta);
// s >= sU
return unscaledSats[2];
}
}
@ -425,28 +439,35 @@ private:
const PointsContainer& unscaledSats,
const PointsContainer& scaledSats)
{
if (unscaledSats[1] >= unscaledSats[2])
return unscaledToScaledSatTwoPoint_(unscaledSat, unscaledSats, scaledSats);
using ScaledSat = std::remove_cv_t<std::remove_reference_t<decltype(scaledSats[0])>>;
if (unscaledSat < unscaledSats[1]) {
Scalar delta = unscaledSats[1] - unscaledSats[0];
if (delta <= 1e-20)
delta = 1.0; // prevent division by zero for (possibly) incorrect input data
auto map = [&unscaledSat, &unscaledSats, &scaledSats](const std::size_t i)
{
const auto distance = (unscaledSat - unscaledSats[i])
/ (unscaledSats[i + 1] - unscaledSats[i]);
return
scaledSats[0]
+
(unscaledSat - unscaledSats[0])*((scaledSats[1] - scaledSats[0])/delta);
const auto displacement =
std::max(scaledSats[i + 1] - scaledSats[i], ScaledSat{ 0 });
return std::min(scaledSats[i] + distance*displacement,
Evaluation { scaledSats[i + 1] });
};
if (! (unscaledSat > unscaledSats[0])) {
return scaledSats[0];
}
else if (unscaledSat < unscaledSats[1]) {
// Tabulated saturation in interval [unscaledSats[0], unscaledSats[1]).
// Map to scaled saturation in [sL, sR).
return map(0);
}
else if (unscaledSat < unscaledSats[2]) {
// Tabulated saturation in interval [unscaledSats[1], unscaledSats[2]).
// Map to scaled saturation in [sR, sU).
return map(1);
}
else {
Scalar delta = unscaledSats[2] - unscaledSats[1];
if (delta <= 1e-20)
delta = 1.0; // prevent division by zero for (possibly) incorrect input data
return
scaledSats[1]
+
(unscaledSat - unscaledSats[1])*((scaledSats[2] - scaledSats[1])/delta);
return scaledSats[2];
}
}