Merge pull request #439 from bska/vscale-disp-sat-general

Implement Three-Point Vertical Scaling for Non-Wetting Phase Relative Permeability
This commit is contained in:
Bård Skaflestad 2021-01-27 21:27:47 +01:00 committed by GitHub
commit d10786a18c
3 changed files with 134 additions and 40 deletions

View File

@ -124,6 +124,20 @@ public:
bool enableThreePointKrwScaling() const
{ return this->enableThreePointKrwScaling_; }
/*!
* \brief Specify whether three-point relative permeability value
* scaling is enabled for the wetting phase (e.g., KRORW + KRO).
*/
void setEnableThreePointKrnScaling(const bool enable)
{ this->enableThreePointKrnScaling_ = enable; }
/*!
* \brief Whether or not three-point relative permeability value scaling
* is enabled for the non-wetting phase (e.g., KRORW + KRO).
*/
bool enableThreePointKrnScaling() const
{ return this->enableThreePointKrnScaling_; }
/*!
* \brief Specify whether relative permeability scaling is enabled for the non-wetting phase.
*/
@ -220,15 +234,20 @@ public:
// check if we are supposed to scale the Y axis of the capillary pressure
if (twoPhaseSystemType == EclOilWaterSystem) {
this->setEnableThreePointKrwScaling(hasKR("WR"));
this->setEnableThreePointKrnScaling(hasKR("ORW"));
this->enableKrnScaling_ = hasKR("O");
this->enableKrnScaling_ = hasKR("O") || this->enableThreePointKrnScaling();
this->enableKrwScaling_ = hasKR("W") || this->enableThreePointKrwScaling();
this->enablePcScaling_ = hasPC("W") || fp.has_double("SWATINIT");
}
else {
assert(twoPhaseSystemType == EclGasOilSystem);
this->enableKrnScaling_ = hasKR("G");
this->enableKrwScaling_ = hasKR("O");
this->setEnableThreePointKrwScaling(hasKR("ORG"));
this->setEnableThreePointKrnScaling(hasKR("GR"));
this->enableKrnScaling_ = hasKR("G") || this->enableThreePointKrnScaling();
this->enableKrwScaling_ = hasKR("O") || this->enableThreePointKrwScaling();
this->enablePcScaling_ = hasPC("G");
}
@ -262,6 +281,9 @@ private:
// Employ three-point vertical scaling (e.g., KRWR and KRW).
bool enableThreePointKrwScaling_{false};
// Employ three-point vertical scaling (e.g., KRORW and KRO).
bool enableThreePointKrnScaling_{false};
};
} // namespace Opm

View File

@ -250,7 +250,7 @@ public:
{
const Evaluation SwUnscaled = scaledToUnscaledSatKrn(params, SwScaled);
const Evaluation krnUnscaled = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), SwUnscaled);
return unscaledToScaledKrn_(params, krnUnscaled);
return unscaledToScaledKrn_(SwScaled, params, krnUnscaled);
}
template <class Evaluation>
@ -581,14 +581,68 @@ private:
* \brief Scale the non-wetting phase relative permeability of a phase according to the given parameters
*/
template <class Evaluation>
static Evaluation unscaledToScaledKrn_(const Params& params, const Evaluation& unscaledKrn)
static Evaluation unscaledToScaledKrn_(const Evaluation& SwScaled,
const Params& params,
const Evaluation& unscaledKrn)
{
if (!params.config().enableKrnScaling())
const auto& cfg = params.config();
if (! cfg.enableKrnScaling())
return unscaledKrn;
//TODO: three point krn y-scaling
Scalar alpha = params.scaledPoints().maxKrn()/params.unscaledPoints().maxKrn();
return unscaledKrn*alpha;
const auto& scaled = params.scaledPoints();
const auto& unscaled = params.unscaledPoints();
if (! cfg.enableThreePointKrnScaling()) {
// Simple case: Run uses pure vertical scaling of non-wetting
// phase's relative permeability (e.g., KRG)
const Scalar alpha = scaled.maxKrn() / unscaled.maxKrn();
return unscaledKrn * alpha;
}
// Otherwise, run uses three-point vertical scaling (e.g., keywords KRGR and KRG)
const auto fdisp = unscaled.krnr();
const auto fmax = unscaled.maxKrn();
const auto sl = scaled.saturationKrnPoints()[0];
const auto sr = std::max(scaled.saturationKrnPoints()[1], sl);
const auto fr = scaled.krnr();
const auto fm = scaled.maxKrn();
// Note logic here. Krn is a decreasing function of Sw (dKrn/dSw <=
// 0) so the roles of left and right intervals are reversed viz
// unscaledToScaledKrw_().
if (! (SwScaled < sr)) {
// Pure vertical scaling in right-hand interval ([SR, SWU])
return unscaledKrn * (fr / fdisp);
}
else if (fmax > fdisp) {
// s \in [SWL, SR), SWL < SR; normal case: Kr(Swl) > Kr(Sr).
//
// Linear function between (sr,fr) and (sl,fm) in terms of
// function value 'unscaledKrn'. This usually alters the shape
// of the relative permeability function in this interval (e.g.,
// roughly quadratic to linear).
const auto t = (unscaledKrn - fdisp) / (fmax - fdisp);
return fr + t*(fm - fr);
}
else if (sr > sl) {
// s \in [SWL, SR), SWL < SR; special case: Kr(Swl) == Kr(Sr).
//
// Linear function between (sr,fr) and (sl,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 = (sr - SwScaled) / (sr - sl);
return fr + t*(fm - fr);
}
else {
// sl == sr (pure scaling). Almost arbitrarily pick 'fm'.
return fm;
}
}
template <class Evaluation>

View File

@ -48,6 +48,7 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableColumn.hpp>
#include <algorithm>
#include <cassert>
@ -703,6 +704,9 @@ private:
// the situation for the gas phase is complicated that all saturations are
// shifted by the connate water saturation.
const Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl;
const auto tolcrit = eclState.runspec().saturationFunctionControls()
.minimumRelpermMobilityThreshold();
const auto& tableManager = eclState.getTableManager();
switch (eclState.runspec().saturationFunctionControls().family()) {
@ -711,12 +715,10 @@ private:
const TableContainer& sgofTables = tableManager.getSgofTables();
const TableContainer& slgofTables = tableManager.getSlgofTables();
if (!sgofTables.empty())
readGasOilEffectiveParametersSgof_(effParams,
Swco,
readGasOilEffectiveParametersSgof_(effParams, Swco, tolcrit,
sgofTables.getTable<SgofTable>(satRegionIdx));
else if (!slgofTables.empty())
readGasOilEffectiveParametersSlgof_(effParams,
Swco,
readGasOilEffectiveParametersSlgof_(effParams, Swco, tolcrit,
slgofTables.getTable<SlgofTable>(satRegionIdx));
break;
}
@ -727,16 +729,11 @@ private:
if (!hasWater) {
// oil and gas case
const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>( satRegionIdx );
readGasOilEffectiveParametersFamily2_(effParams,
Swco,
sof2Table,
sgfnTable);
} else {
readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable);
}
else {
const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
readGasOilEffectiveParametersFamily2_(effParams,
Swco,
sof3Table,
sgfnTable);
readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable);
}
break;
}
@ -747,7 +744,8 @@ private:
}
void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
Scalar Swco,
const Scalar Swco,
const double tolcrit,
const Opm::SgofTable& sgofTable)
{
// convert the saturations of the SGOF keyword from gas to oil saturations
@ -756,14 +754,15 @@ private:
SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx);
}
effParams.setKrwSamples(SoSamples, sgofTable.getColumn("KROG").vectorCopy());
effParams.setKrnSamples(SoSamples, sgofTable.getColumn("KRG").vectorCopy());
effParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KROG")));
effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KRG")));
effParams.setPcnwSamples(SoSamples, sgofTable.getColumn("PCOG").vectorCopy());
effParams.finalize();
}
void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
Scalar Swco,
const Scalar Swco,
const double tolcrit,
const Opm::SlgofTable& slgofTable)
{
// convert the saturations of the SLGOF keyword from "liquid" to oil saturations
@ -772,14 +771,15 @@ private:
SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco;
}
effParams.setKrwSamples(SoSamples, slgofTable.getColumn("KROG").vectorCopy());
effParams.setKrnSamples(SoSamples, slgofTable.getColumn("KRG").vectorCopy());
effParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KROG")));
effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KRG")));
effParams.setPcnwSamples(SoSamples, slgofTable.getColumn("PCOG").vectorCopy());
effParams.finalize();
}
void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
Scalar Swco,
const Scalar Swco,
const double tolcrit,
const Opm::Sof3Table& sof3Table,
const Opm::SgfnTable& sgfnTable)
{
@ -790,14 +790,15 @@ private:
SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
}
effParams.setKrwSamples(SoColumn, sof3Table.getColumn("KROG").vectorCopy());
effParams.setKrnSamples(SoSamples, sgfnTable.getColumn("KRG").vectorCopy());
effParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROG")));
effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
effParams.finalize();
}
void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
Scalar Swco,
const Scalar Swco,
const double tolcrit,
const Opm::Sof2Table& sof2Table,
const Opm::SgfnTable& sgfnTable)
{
@ -808,8 +809,8 @@ private:
SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
}
effParams.setKrwSamples(SoColumn, sof2Table.getColumn("KRO").vectorCopy());
effParams.setKrnSamples(SoSamples, sgfnTable.getColumn("KRG").vectorCopy());
effParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO")));
effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
effParams.finalize();
}
@ -825,6 +826,9 @@ private:
dest[satRegionIdx] = std::make_shared<OilWaterEffectiveTwoPhaseParams>();
const auto tolcrit = eclState.runspec().saturationFunctionControls()
.minimumRelpermMobilityThreshold();
const auto& tableManager = eclState.getTableManager();
auto& effParams = *dest[satRegionIdx];
@ -832,10 +836,10 @@ private:
case SatFuncControls::KeywordFamily::Family_I:
{
const auto& swofTable = tableManager.getSwofTables().getTable<SwofTable>(satRegionIdx);
std::vector<double> SwColumn = swofTable.getColumn("SW").vectorCopy();
const std::vector<double> SwColumn = swofTable.getColumn("SW").vectorCopy();
effParams.setKrwSamples(SwColumn, swofTable.getColumn("KRW").vectorCopy());
effParams.setKrnSamples(SwColumn, swofTable.getColumn("KROW").vectorCopy());
effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KRW")));
effParams.setKrnSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KROW")));
effParams.setPcnwSamples(SwColumn, swofTable.getColumn("PCOW").vectorCopy());
effParams.finalize();
break;
@ -845,15 +849,15 @@ private:
{
const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
const std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
// convert the saturations of the SOF3 keyword from oil to water saturations
std::vector<double> SwSamples(sof3Table.numRows());
for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx);
effParams.setKrwSamples(SwColumn, swfnTable.getColumn("KRW").vectorCopy());
effParams.setKrnSamples(SwSamples, sof3Table.getColumn("KROW").vectorCopy());
effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROW")));
effParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
effParams.finalize();
break;
@ -985,6 +989,20 @@ private:
}
}
// Relative permeability values not strictly greater than 'tolcrit' treated as zero.
std::vector<double> normalizeKrValues_(const double tolcrit,
const TableColumn& krValues) const
{
auto kr = krValues.vectorCopy();
std::transform(kr.begin(), kr.end(), kr.begin(),
[tolcrit](const double kri)
{
return (kri > tolcrit) ? kri : 0.0;
});
return kr;
}
bool enableEndPointScaling_;
std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;