Merge pull request #436 from bska/vscale-disp-sat
Implement KRWR Vertical Scaling
This commit is contained in:
Bård Skaflestad
GitHub
commit
292ee5174a
@ -110,6 +110,20 @@ public:
|
||||
bool enableKrwScaling() const
|
||||
{ return enableKrwScaling_; }
|
||||
|
||||
/*!
|
||||
* \brief Specify whether three-point relative permeability value
|
||||
* scaling is enabled for the wetting phase (KRWR + KRW).
|
||||
*/
|
||||
void setEnableThreePointKrwScaling(const bool enable)
|
||||
{ this->enableThreePointKrwScaling_ = enable; }
|
||||
|
||||
/*!
|
||||
* \brief Whether or not three-point relative permeability value scaling
|
||||
* is enabled for the wetting phase (KRWR + KRW).
|
||||
*/
|
||||
bool enableThreePointKrwScaling() const
|
||||
{ return this->enableThreePointKrwScaling_; }
|
||||
|
||||
/*!
|
||||
* \brief Specify whether relative permeability scaling is enabled for the non-wetting phase.
|
||||
*/
|
||||
@ -161,7 +175,9 @@ public:
|
||||
* This requires that the opm-parser module is available.
|
||||
*/
|
||||
void initFromState(const Opm::EclipseState& eclState,
|
||||
Opm::EclTwoPhaseSystemType twoPhaseSystemType)
|
||||
Opm::EclTwoPhaseSystemType twoPhaseSystemType,
|
||||
const std::string& prefix = "",
|
||||
const std::string& suffix = "")
|
||||
{
|
||||
const auto& endscale = eclState.runspec().endpointScaling();
|
||||
// find out if endpoint scaling is used in the first place
|
||||
@ -192,16 +208,28 @@ public:
|
||||
}
|
||||
|
||||
const auto& fp = eclState.fieldProps();
|
||||
auto hasKR = [&fp, &prefix, &suffix](const std::string& scaling)
|
||||
{
|
||||
return fp.has_double(prefix + "KR" + scaling + suffix);
|
||||
};
|
||||
auto hasPC = [&fp, &prefix](const std::string& scaling)
|
||||
{
|
||||
return fp.has_double(prefix + "PC" + scaling);
|
||||
};
|
||||
|
||||
// check if we are supposed to scale the Y axis of the capillary pressure
|
||||
if (twoPhaseSystemType == EclOilWaterSystem) {
|
||||
enableKrnScaling_ = fp.has_double("KRO");
|
||||
enableKrwScaling_ = fp.has_double("KRW");
|
||||
enablePcScaling_ = fp.has_double("PCW") || fp.has_double("SWATINIT");
|
||||
} else {
|
||||
this->setEnableThreePointKrwScaling(hasKR("WR"));
|
||||
|
||||
this->enableKrnScaling_ = hasKR("O");
|
||||
this->enableKrwScaling_ = hasKR("W") || this->enableThreePointKrwScaling();
|
||||
this->enablePcScaling_ = hasPC("W") || fp.has_double("SWATINIT");
|
||||
}
|
||||
else {
|
||||
assert(twoPhaseSystemType == EclGasOilSystem);
|
||||
enableKrnScaling_ = fp.has_double("KRG");
|
||||
enableKrwScaling_ = fp.has_double("KRO");
|
||||
enablePcScaling_ = fp.has_double("PCG");
|
||||
this->enableKrnScaling_ = hasKR("G");
|
||||
this->enableKrwScaling_ = hasKR("O");
|
||||
this->enablePcScaling_ = hasPC("G");
|
||||
}
|
||||
|
||||
if (enablePcScaling_ && enableLeverettScaling_) {
|
||||
@ -231,6 +259,9 @@ private:
|
||||
bool enableLeverettScaling_;
|
||||
bool enableKrwScaling_;
|
||||
bool enableKrnScaling_;
|
||||
|
||||
// Employ three-point vertical scaling (e.g., KRWR and KRW).
|
||||
bool enableThreePointKrwScaling_{false};
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
@ -99,8 +99,12 @@ public:
|
||||
this->compressed_pcw = try_get( fp, kwPrefix+"PCW");
|
||||
this->compressed_pcg = try_get( fp, kwPrefix+"PCG");
|
||||
this->compressed_krw = try_get( fp, kwPrefix+"KRW");
|
||||
this->compressed_krwr = try_get( fp, kwPrefix+"KRWR");
|
||||
this->compressed_kro = try_get( fp, kwPrefix+"KRO");
|
||||
this->compressed_krorg = try_get( fp, kwPrefix+"KRORG");
|
||||
this->compressed_krorw = try_get( fp, kwPrefix+"KRORW");
|
||||
this->compressed_krg = try_get( fp, kwPrefix+"KRG");
|
||||
this->compressed_krgr = try_get( fp, kwPrefix+"KRGR");
|
||||
|
||||
// _may_ be needed to calculate the Leverett capillary pressure scaling factor
|
||||
if (fp.has_double("PORO"))
|
||||
@ -185,14 +189,29 @@ public:
|
||||
return this->satfunc(this->compressed_krw, active_index);
|
||||
}
|
||||
|
||||
const double * krwr(std::size_t active_index) const {
|
||||
return this->satfunc(this->compressed_krwr, active_index);
|
||||
}
|
||||
|
||||
const double * krg(std::size_t active_index) const {
|
||||
return this->satfunc(this->compressed_krg, active_index);
|
||||
}
|
||||
|
||||
const double * krgr(std::size_t active_index) const {
|
||||
return this->satfunc(this->compressed_krgr, active_index);
|
||||
}
|
||||
|
||||
const double * kro(std::size_t active_index) const {
|
||||
return this->satfunc(this->compressed_kro, active_index);
|
||||
}
|
||||
|
||||
const double * krorg(std::size_t active_index) const {
|
||||
return this->satfunc(this->compressed_krorg, active_index);
|
||||
}
|
||||
|
||||
const double * krorw(std::size_t active_index) const {
|
||||
return this->satfunc(this->compressed_krorw, active_index);
|
||||
}
|
||||
|
||||
private:
|
||||
const double *
|
||||
@ -215,8 +234,12 @@ private:
|
||||
std::vector<double> compressed_pcw;
|
||||
std::vector<double> compressed_pcg;
|
||||
std::vector<double> compressed_krw;
|
||||
std::vector<double> compressed_krwr;
|
||||
std::vector<double> compressed_kro;
|
||||
std::vector<double> compressed_krorg;
|
||||
std::vector<double> compressed_krorw;
|
||||
std::vector<double> compressed_krg;
|
||||
std::vector<double> compressed_krgr;
|
||||
|
||||
std::vector<double> compressed_permx;
|
||||
std::vector<double> compressed_permy;
|
||||
|
||||
@ -83,6 +83,12 @@ struct EclEpsScalingPointsInfo
|
||||
Scalar pcowLeverettFactor;
|
||||
Scalar pcgoLeverettFactor;
|
||||
|
||||
// Scaled relative permeabilities at residual displacing saturation
|
||||
Scalar Krwr; // water
|
||||
Scalar Krgr; // gas
|
||||
Scalar Krorw; // oil in water-oil system
|
||||
Scalar Krorg; // oil in gas-oil system
|
||||
|
||||
// maximum relative permabilities
|
||||
Scalar maxKrw; // maximum relative permability of water
|
||||
Scalar maxKrow; // maximum relative permability of oil in the oil-water system
|
||||
@ -103,6 +109,10 @@ struct EclEpsScalingPointsInfo
|
||||
maxPcgo == data.maxPcgo &&
|
||||
pcowLeverettFactor == data.pcowLeverettFactor &&
|
||||
pcgoLeverettFactor == data.pcgoLeverettFactor &&
|
||||
Krwr == data.Krwr &&
|
||||
Krgr == data.Krgr &&
|
||||
Krorw == data.Krorw &&
|
||||
Krorg == data.Krorg &&
|
||||
maxKrw == data.maxKrw &&
|
||||
maxKrow == data.maxKrow &&
|
||||
maxKrog == data.maxKrog &&
|
||||
@ -123,6 +133,10 @@ struct EclEpsScalingPointsInfo
|
||||
<< " maxPcgo: " << maxPcgo << '\n'
|
||||
<< " pcowLeverettFactor: " << pcowLeverettFactor << '\n'
|
||||
<< " pcgoLeverettFactor: " << pcgoLeverettFactor << '\n'
|
||||
<< " Krwr: " << Krwr << '\n'
|
||||
<< " Krgr: " << Krgr << '\n'
|
||||
<< " Krorw: " << Krorw << '\n'
|
||||
<< " Krorg: " << Krorg << '\n'
|
||||
<< " maxKrw: " << maxKrw << '\n'
|
||||
<< " maxKrg: " << maxKrg << '\n'
|
||||
<< " maxKrow: " << maxKrow << '\n'
|
||||
@ -163,6 +177,11 @@ struct EclEpsScalingPointsInfo
|
||||
this->pcowLeverettFactor = 1.0;
|
||||
this->pcgoLeverettFactor = 1.0;
|
||||
|
||||
this->Krwr = rfunc.krw.r [satRegionIdx];
|
||||
this->Krgr = rfunc.krg.r [satRegionIdx];
|
||||
this->Krorw = rfunc.kro.rw[satRegionIdx];
|
||||
this->Krorg = rfunc.kro.rg[satRegionIdx];
|
||||
|
||||
this->maxKrw = rfunc.krw.max[satRegionIdx];
|
||||
this->maxKrow = rfunc.kro.max[satRegionIdx];
|
||||
this->maxKrog = rfunc.kro.max[satRegionIdx];
|
||||
@ -196,9 +215,14 @@ struct EclEpsScalingPointsInfo
|
||||
update(Sgu, epsProperties.sgu(activeIndex));
|
||||
update(maxPcow, epsProperties.pcw(activeIndex));
|
||||
update(maxPcgo, epsProperties.pcg(activeIndex));
|
||||
|
||||
update(this->Krwr, epsProperties.krwr(activeIndex));
|
||||
update(this->Krgr, epsProperties.krgr(activeIndex));
|
||||
update(this->Krorw, epsProperties.krorw(activeIndex));
|
||||
update(this->Krorg, epsProperties.krorg(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,37 +322,29 @@ 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;
|
||||
else
|
||||
maxPcnwOrLeverettFactor_ = epsInfo.maxPcow;
|
||||
|
||||
Krwr_ = epsInfo.Krwr;
|
||||
Krnr_ = epsInfo.Krorw;
|
||||
|
||||
maxKrw_ = epsInfo.maxKrw;
|
||||
maxKrn_ = epsInfo.maxKrow;
|
||||
}
|
||||
@ -337,38 +353,29 @@ 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;
|
||||
else
|
||||
maxPcnwOrLeverettFactor_ = epsInfo.maxPcgo;
|
||||
|
||||
Krwr_ = epsInfo.Krorg;
|
||||
Krnr_ = epsInfo.Krgr;
|
||||
|
||||
maxKrw_ = epsInfo.maxKrog;
|
||||
maxKrn_ = epsInfo.maxKrg;
|
||||
}
|
||||
@ -383,7 +390,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_; }
|
||||
|
||||
/*!
|
||||
@ -434,6 +441,20 @@ public:
|
||||
Scalar leverettFactor() const
|
||||
{ return maxPcnwOrLeverettFactor_; }
|
||||
|
||||
/*!
|
||||
* \brief Set wetting-phase relative permeability at residual saturation
|
||||
* of non-wetting phase.
|
||||
*/
|
||||
void setKrwr(Scalar value)
|
||||
{ this->Krwr_ = value; }
|
||||
|
||||
/*!
|
||||
* \brief Returns wetting-phase relative permeability at residual
|
||||
* saturation of non-wetting phase.
|
||||
*/
|
||||
Scalar krwr() const
|
||||
{ return this->Krwr_; }
|
||||
|
||||
/*!
|
||||
* \brief Sets the maximum wetting phase relative permeability
|
||||
*/
|
||||
@ -446,6 +467,20 @@ public:
|
||||
Scalar maxKrw() const
|
||||
{ return maxKrw_; }
|
||||
|
||||
/*!
|
||||
* \brief Set non-wetting phase relative permeability at residual
|
||||
* saturation of wetting phase.
|
||||
*/
|
||||
void setKrnr(Scalar value)
|
||||
{ this->Krnr_ = value; }
|
||||
|
||||
/*!
|
||||
* \brief Returns non-wetting phase relative permeability at residual
|
||||
* saturation of wetting phase.
|
||||
*/
|
||||
Scalar krnr() const
|
||||
{ return this->Krnr_; }
|
||||
|
||||
/*!
|
||||
* \brief Sets the maximum wetting phase relative permeability
|
||||
*/
|
||||
@ -466,17 +501,25 @@ 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
|
||||
// Scaled wetting phase relative permeability value at residual
|
||||
// saturation of non-wetting phase.
|
||||
Scalar Krwr_;
|
||||
|
||||
// Maximum non-wetting phase relative permability value
|
||||
Scalar maxKrn_;
|
||||
|
||||
// Scaled non-wetting phase relative permeability value at residual
|
||||
// saturation of wetting phase.
|
||||
Scalar Krnr_;
|
||||
|
||||
// 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_;
|
||||
|
||||
@ -31,6 +31,10 @@
|
||||
|
||||
#include <opm/material/fluidstates/SaturationOverlayFluidState.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <cstddef>
|
||||
#include <type_traits>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup FluidMatrixInteractions
|
||||
@ -221,7 +225,7 @@ public:
|
||||
{
|
||||
const Evaluation& SwUnscaled = scaledToUnscaledSatKrw(params, SwScaled);
|
||||
const Evaluation& krwUnscaled = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), SwUnscaled);
|
||||
return unscaledToScaledKrw_(params, krwUnscaled);
|
||||
return unscaledToScaledKrw_(SwScaled, params, krwUnscaled);
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
@ -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];
|
||||
}
|
||||
}
|
||||
|
||||
@ -487,14 +508,63 @@ private:
|
||||
* \brief Scale the wetting phase relative permeability of a phase according to the given parameters
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation unscaledToScaledKrw_(const Params& params, const Evaluation& unscaledKrw)
|
||||
static Evaluation unscaledToScaledKrw_(const Evaluation& SwScaled,
|
||||
const Params& params,
|
||||
const Evaluation& unscaledKrw)
|
||||
{
|
||||
if (!params.config().enableKrwScaling())
|
||||
const auto& cfg = params.config();
|
||||
|
||||
if (! cfg.enableKrwScaling())
|
||||
return unscaledKrw;
|
||||
|
||||
// TODO: three point krw y-scaling
|
||||
Scalar alpha = params.scaledPoints().maxKrw()/params.unscaledPoints().maxKrw();
|
||||
return unscaledKrw*alpha;
|
||||
const auto& scaled = params.scaledPoints();
|
||||
const auto& unscaled = params.unscaledPoints();
|
||||
|
||||
if (! cfg.enableThreePointKrwScaling()) {
|
||||
// Simple case: Run uses pure vertical scaling of water relperm (keyword KRW)
|
||||
const Scalar alpha = scaled.maxKrw() / unscaled.maxKrw();
|
||||
return unscaledKrw * alpha;
|
||||
}
|
||||
|
||||
// Otherwise, run uses three-point vertical scaling (keywords KRWR and KRW)
|
||||
const auto fdisp = unscaled.krwr();
|
||||
const auto fmax = unscaled.maxKrw();
|
||||
|
||||
const auto sm = scaled.saturationKrwPoints()[2];
|
||||
const auto sr = std::min(scaled.saturationKrwPoints()[1], sm);
|
||||
const auto fr = scaled.krwr();
|
||||
const auto fm = scaled.maxKrw();
|
||||
|
||||
if (! (SwScaled > sr)) {
|
||||
// Pure vertical scaling in left interval ([SWL, SR])
|
||||
return unscaledKrw * (fr / fdisp);
|
||||
}
|
||||
else if (fmax > fdisp) {
|
||||
// s \in [sr, sm), sm > sr; normal case: Kr(Smax) > Kr(Sr).
|
||||
//
|
||||
// Linear function between (sr,fr) and (sm,fm) in terms of
|
||||
// function value 'unscaledKrw'. This usually alters the shape
|
||||
// of the relative permeability function in this interval (e.g.,
|
||||
// roughly quadratic to linear).
|
||||
const auto t = (unscaledKrw - fdisp) / (fmax - fdisp);
|
||||
|
||||
return fr + t*(fm - fr);
|
||||
}
|
||||
else if (sr < sm) {
|
||||
// s \in [sr, sm), sm > sr; special case: Kr(Smax) == Kr(Sr).
|
||||
//
|
||||
// Linear function between (sr,fr) and (sm,fm) in terms of
|
||||
// saturation value 'SwScaled'. This usually alters the shape
|
||||
// of the relative permeability function in this interval (e.g.,
|
||||
// roughly quadratic to linear).
|
||||
const auto t = (SwScaled - sr) / (sm - sr);
|
||||
|
||||
return fr + t*(fm - fr);
|
||||
}
|
||||
else {
|
||||
// sm == sr (pure scaling). Almost arbitrarily pick 'fm'.
|
||||
return fm;
|
||||
}
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
|
||||
Loading…
Reference in New Issue
Block a user