mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
4be1685604
The oil-water capillary pressure is a non-increasing function of water saturation so we need to ensure that we use the capillary pressure value at the minimum tabulated water saturation. Add special purpose code to enforce this rule. This adds a certain amount of computational overhead because we now compute function values at two saturation points instead of just one for the case of pure vertical scaling. Most of the time those extra function values will just be subsequently discarded. Note that this is a bit of hack, because it relies on the fact that the current implementation assigns TableEndPoints::disp = TableEndPoints::low in the case of two-point horizontal scaling. We may wish to make that rule more explicit.
2583 lines
82 KiB
C++
2583 lines
82 KiB
C++
/*
|
|
Copyright 2018 Equinor ASA.
|
|
Copyright 2017 Statoil ASA.
|
|
|
|
This file is part of the Open Porous Media Project (OPM).
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#if HAVE_CONFIG_H
|
|
#include <config.h>
|
|
#endif
|
|
|
|
#include <opm/utility/ECLEndPointScaling.hpp>
|
|
|
|
#include <opm/utility/ECLGraph.hpp>
|
|
#include <opm/utility/ECLPhaseIndex.hpp>
|
|
#include <opm/utility/ECLResultData.hpp>
|
|
|
|
#include <opm/utility/ECLUnitHandling.hpp>
|
|
#include <opm/utility/imported/Units.hpp>
|
|
|
|
#include <algorithm>
|
|
#include <cassert>
|
|
#include <cmath>
|
|
#include <exception>
|
|
#include <iterator>
|
|
#include <limits>
|
|
#include <numeric>
|
|
#include <stdexcept>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include <ert/ecl/ecl_kw_magic.h>
|
|
|
|
namespace {
|
|
std::vector<Opm::SatFunc::EPSEvalInterface::TableEndPoints>
|
|
unscaledTwoPt(const std::vector<double>& min,
|
|
const std::vector<double>& max)
|
|
{
|
|
assert ((min.size() == max.size()) && "Internal Logic Error");
|
|
assert ((! min.empty()) && "Internal Logic Error");
|
|
|
|
using TEP = Opm::SatFunc::EPSEvalInterface::TableEndPoints;
|
|
|
|
auto tep = std::vector<TEP>{};
|
|
tep.reserve(min.size());
|
|
|
|
for (auto n = min.size(), i = 0*n; i < n; ++i) {
|
|
const auto smin = min[i];
|
|
|
|
// Ignore 'sdisp' in the two-point scaling.
|
|
tep.push_back(TEP{ smin, smin, max[i] });
|
|
}
|
|
|
|
return tep;
|
|
}
|
|
|
|
std::vector<Opm::SatFunc::EPSEvalInterface::TableEndPoints>
|
|
unscaledThreePt(const std::vector<double>& min ,
|
|
const std::vector<double>& disp,
|
|
const std::vector<double>& max )
|
|
{
|
|
assert ((min.size() == max .size()) && "Internal Logic Error");
|
|
assert ((min.size() == disp.size()) && "Internal Logic Error");
|
|
assert ((! min.empty()) && "Internal Logic Error");
|
|
|
|
using TEP = Opm::SatFunc::EPSEvalInterface::TableEndPoints;
|
|
|
|
auto tep = std::vector<TEP>{};
|
|
tep.reserve(min.size());
|
|
|
|
for (auto n = min.size(), i = 0*n; i < n; ++i) {
|
|
tep.push_back(TEP{ min[i], disp[i], max[i] });
|
|
}
|
|
|
|
return tep;
|
|
}
|
|
|
|
bool haveKeywordData(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::string& vector)
|
|
{
|
|
const auto& grids = G.activeGrids();
|
|
|
|
return std::any_of(std::begin(grids), std::end(grids),
|
|
[&init, &vector](const std::string& gridID)
|
|
{
|
|
return init.haveKeywordData(vector, gridID);
|
|
});
|
|
}
|
|
|
|
template <class CvrtVal>
|
|
std::vector<double>
|
|
gridDefaultedVector(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::string& vector,
|
|
const std::vector<double>& dflt,
|
|
CvrtVal&& cvrt,
|
|
const std::vector<double>& fallback
|
|
= std::vector<double>{})
|
|
{
|
|
auto ret = std::vector<double>(G.numCells());
|
|
|
|
const auto sentinel = 1.0e20;
|
|
|
|
assert (! dflt.empty() && "Internal Error");
|
|
|
|
auto cellID = std::vector<double>::size_type{0};
|
|
for (const auto& gridID : G.activeGrids()) {
|
|
const auto nc = G.numCells(gridID);
|
|
|
|
const auto& snum = init.haveKeywordData("SATNUM", gridID)
|
|
? G.rawLinearisedCellData<int>(init, "SATNUM", gridID)
|
|
: std::vector<int>(nc, 1);
|
|
|
|
const auto& val = init.haveKeywordData(vector, gridID)
|
|
? G.rawLinearisedCellData<double>(init, vector, gridID)
|
|
: std::vector<double>(nc, -1.0e21);
|
|
|
|
for (auto c = 0*nc; c < nc; ++c, ++cellID) {
|
|
const auto fb = fallback.empty()
|
|
? -sentinel : fallback[c];
|
|
|
|
auto& r = ret[cellID];
|
|
|
|
if (std::abs(val[c]) < sentinel) {
|
|
r = cvrt(val[c]);
|
|
}
|
|
else if (std::abs(fb) < sentinel) {
|
|
r = cvrt(fb);
|
|
}
|
|
else {
|
|
r = dflt[snum[c] - 1];
|
|
}
|
|
}
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
std::vector<double>
|
|
sgCrit(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SGCR", tep.crit.gas,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
std::vector<double>
|
|
sogCrit(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SOGCR", tep.crit.oil_in_gas,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
std::vector<double>
|
|
sowCrit(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SOWCR", tep.crit.oil_in_water,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
std::vector<double>
|
|
swCrit(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SWCR", tep.crit.water,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
std::vector<double>
|
|
sgMax(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SGU", tep.smax.gas,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
std::vector<double>
|
|
soMax(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
auto smax = std::vector<double>(G.numCells());
|
|
|
|
const auto sgl = ::Opm::SatFunc::scaledConnateGas (G, init, tep);
|
|
const auto swl = ::Opm::SatFunc::scaledConnateWater(G, init, tep);
|
|
|
|
std::transform(std::begin(sgl), std::end(sgl), std::begin(swl),
|
|
std::begin(smax),
|
|
[](const double sg, const double sw) -> double
|
|
{
|
|
return 1.0 - (sg + sw);
|
|
});
|
|
|
|
return smax;
|
|
}
|
|
|
|
std::vector<double>
|
|
swMax(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SWU", tep.smax.water,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
double defaultedScaledSaturation(const double s, const double dflt)
|
|
{
|
|
// Use input scaled saturation ('s') if not defaulted (|s| < 1e20),
|
|
// default scaled saturation otherwise. The sentinel value 1e20 is
|
|
// the common value used to represent unset/defaulted values in ECL
|
|
// result sets.
|
|
|
|
return (std::abs(s) < 1.0e+20) ? s : dflt;
|
|
}
|
|
|
|
bool
|
|
haveScaledRelPermAtCritSat(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::EPSOptions& opt)
|
|
{
|
|
switch (opt.thisPh) {
|
|
case ::Opm::ECLPhaseIndex::Aqua:
|
|
return haveKeywordData(G, init, "KRWR");
|
|
|
|
case ::Opm::ECLPhaseIndex::Liquid:
|
|
return (opt.subSys == ::Opm::SatFunc::CreateEPS::SubSystem::OilGas)
|
|
? haveKeywordData(G, init, "KRORG")
|
|
: haveKeywordData(G, init, "KRORW");
|
|
|
|
case ::Opm::ECLPhaseIndex::Vapour:
|
|
return haveKeywordData(G, init, "KRGR");
|
|
}
|
|
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Class Opm::TwoPointScaling::Impl
|
|
// ---------------------------------------------------------------------
|
|
|
|
class Opm::SatFunc::TwoPointScaling::Impl
|
|
{
|
|
public:
|
|
Impl(std::vector<double> smin,
|
|
std::vector<double> smax)
|
|
: smin_(std::move(smin)), smax_(std::move(smax))
|
|
{
|
|
if (this->smin_.size() != this->smax_.size()) {
|
|
throw std::invalid_argument {
|
|
"Size Mismatch Between Minimum and "
|
|
"Maximum Saturation Arrays"
|
|
};
|
|
}
|
|
}
|
|
|
|
std::vector<double>
|
|
eval(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const;
|
|
|
|
std::vector<double>
|
|
reverse(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const;
|
|
|
|
private:
|
|
std::vector<double> smin_;
|
|
std::vector<double> smax_;
|
|
|
|
double sMin(const std::vector<int>::size_type cell,
|
|
const TableEndPoints& tep) const
|
|
{
|
|
// Use model's scaled connate saturation if defined, otherwise fall
|
|
// back to table's unscaled connate saturation.
|
|
|
|
return defaultedScaledSaturation(this->smin_[cell], tep.low);
|
|
}
|
|
|
|
double sMax(const std::vector<int>::size_type cell,
|
|
const TableEndPoints& tep) const
|
|
{
|
|
// Use model's scaled maximum saturation if defined, otherwise fall
|
|
// back to table's unscaled maximum saturation.
|
|
|
|
return defaultedScaledSaturation(this->smax_[cell], tep.high);
|
|
}
|
|
};
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::TwoPointScaling::
|
|
Impl::eval(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
const auto srng = tep.high - tep.low;
|
|
|
|
auto effsat = std::vector<double>{};
|
|
effsat.reserve(sp.size());
|
|
|
|
for (const auto& eval_pt : sp) {
|
|
const auto cell = eval_pt.cell;
|
|
|
|
const auto sLO = this->sMin(cell, tep);
|
|
const auto sHI = this->sMax(cell, tep);
|
|
|
|
effsat.push_back(0.0);
|
|
auto& s_eff = effsat.back();
|
|
|
|
if (! (eval_pt.sat > sLO)) {
|
|
// s <= sLO
|
|
s_eff = tep.low;
|
|
}
|
|
else if (! (eval_pt.sat < sHI)) {
|
|
// s >= sHI
|
|
s_eff = tep.high;
|
|
}
|
|
else {
|
|
// s \in (sLO, sHI)
|
|
const auto scaled_sat =
|
|
((eval_pt.sat - sLO) / (sHI - sLO)) * srng;
|
|
|
|
s_eff = tep.low + scaled_sat;
|
|
}
|
|
}
|
|
|
|
return effsat;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::TwoPointScaling::
|
|
Impl::reverse(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
const auto srng = tep.high - tep.low;
|
|
|
|
auto unscaledsat = std::vector<double>{};
|
|
unscaledsat.reserve(sp.size());
|
|
|
|
for (const auto& eval_pt : sp) {
|
|
const auto cell = eval_pt.cell;
|
|
|
|
const auto sLO = this->sMin(cell, tep);
|
|
const auto sHI = this->sMax(cell, tep);
|
|
|
|
unscaledsat.push_back(0.0);
|
|
auto& s_unsc = unscaledsat.back();
|
|
|
|
if (! (eval_pt.sat > tep.low)) {
|
|
// s <= minimum tabulated saturation.
|
|
// Map to Minimum Input Saturation in cell (sLO).
|
|
s_unsc = sLO;
|
|
}
|
|
else if (! (eval_pt.sat < tep.high)) {
|
|
// s >= maximum tabulated saturation.
|
|
// Map to Maximum Input Saturation in cell (sHI).
|
|
s_unsc = sHI;
|
|
}
|
|
else {
|
|
// s in tabulated interval (tep.low, tep.high)
|
|
// Map to Input Saturation in (sLO, sHI)
|
|
const auto t =
|
|
(eval_pt.sat - tep.low) / srng;
|
|
|
|
s_unsc = sLO + t*(sHI - sLO);
|
|
}
|
|
}
|
|
|
|
return unscaledsat;
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Class Opm::SatFunc::PureVerticalScaling::Impl
|
|
// ---------------------------------------------------------------------
|
|
|
|
class Opm::SatFunc::PureVerticalScaling::Impl
|
|
{
|
|
public:
|
|
explicit Impl(std::vector<double> fmax)
|
|
: fmax_(std::move(fmax))
|
|
, inf_ (std::numeric_limits<double>::has_infinity
|
|
? std::numeric_limits<double>::infinity()
|
|
: std::numeric_limits<double>::max())
|
|
{}
|
|
|
|
std::vector<double>
|
|
vertScale(const FunctionValues& f,
|
|
const SaturationPoints& sp,
|
|
std::vector<double> val) const;
|
|
|
|
private:
|
|
std::vector<double> fmax_;
|
|
|
|
double inf_;
|
|
|
|
std::vector<double>
|
|
zeroMaxVal(const SaturationPoints& sp) const;
|
|
|
|
std::vector<double>
|
|
nonZeroMaxVal(const SaturationPoints& sp,
|
|
const double maxVal,
|
|
std::vector<double> val) const;
|
|
};
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::PureVerticalScaling::Impl::
|
|
vertScale(const FunctionValues& f,
|
|
const SaturationPoints& sp,
|
|
std::vector<double> val) const
|
|
{
|
|
assert (sp.size() == val.size() && "Internal Error in Vertical Scaling");
|
|
|
|
const auto maxVal = f.max.val;
|
|
|
|
if (! (std::abs(maxVal) > 0.0)) {
|
|
return this->zeroMaxVal(sp);
|
|
}
|
|
|
|
return this->nonZeroMaxVal(sp, maxVal, std::move(val));
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::PureVerticalScaling::Impl::
|
|
zeroMaxVal(const SaturationPoints& sp) const
|
|
{
|
|
auto ret = std::vector<double>(sp.size(), 0.0);
|
|
|
|
for (auto n = sp.size(), i = 0*n; i < n; ++i) {
|
|
const auto fmax = this->fmax_[ sp[i].cell ];
|
|
|
|
ret[i] = (std::abs(fmax) > 0.0)
|
|
? (std::signbit(fmax) ? -this->inf_ : this->inf_)
|
|
: 0.0;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::PureVerticalScaling::Impl::
|
|
nonZeroMaxVal(const SaturationPoints& sp,
|
|
const double maxVal,
|
|
std::vector<double> val) const
|
|
{
|
|
auto ret = std::move(val);
|
|
|
|
for (auto n = sp.size(), i = 0*n; i < n; ++i) {
|
|
ret[i] *= this->fmax_[ sp[i].cell ] / maxVal;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Class Opm::SatFunc::CritSatVerticalScaling::Impl
|
|
// ---------------------------------------------------------------------
|
|
|
|
class Opm::SatFunc::CritSatVerticalScaling::Impl
|
|
{
|
|
public:
|
|
explicit Impl(std::vector<double> sdisp,
|
|
std::vector<double> fdisp,
|
|
std::vector<double> smax,
|
|
std::vector<double> fmax)
|
|
: sdisp_(std::move(sdisp))
|
|
, fdisp_(std::move(fdisp))
|
|
, smax_ (std::move(smax))
|
|
, fmax_ (std::move(fmax))
|
|
{}
|
|
|
|
std::vector<double>
|
|
vertScale(const FunctionValues& f,
|
|
const SaturationPoints& sp,
|
|
std::vector<double> val) const;
|
|
|
|
private:
|
|
std::vector<double> sdisp_;
|
|
std::vector<double> fdisp_;
|
|
std::vector<double> smax_;
|
|
std::vector<double> fmax_;
|
|
};
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::CritSatVerticalScaling::Impl::
|
|
vertScale(const FunctionValues& f,
|
|
const SaturationPoints& sp,
|
|
std::vector<double> val) const
|
|
{
|
|
assert ((sp.size() == val.size()) && "Internal Error in Vertical Scaling");
|
|
|
|
auto ret = std::move(val);
|
|
|
|
const auto fdisp = f.disp.val;
|
|
const auto fmax = std::max(f.disp.val, f.max.val);
|
|
const auto sepfv = fmax > fdisp;
|
|
|
|
for (auto n = sp.size(), i = 0*n; i < n; ++i) {
|
|
auto& y = ret[i];
|
|
|
|
const auto c = sp[i].cell;
|
|
const auto s = sp[i].sat;
|
|
const auto sr = std::min(this->sdisp_[c], this->smax_[c]);
|
|
const auto fr = this->fdisp_[c];
|
|
const auto sm = this->smax_ [c];
|
|
const auto fm = this->fmax_ [c];
|
|
|
|
if (s < sr) {
|
|
// s < sr: Pure vertical scaling in left interval.
|
|
y *= fr / fdisp;
|
|
}
|
|
else if (sepfv) {
|
|
// 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 'y'. This usually alters the shape of the
|
|
// relative permeability function in this interval (e.g.,
|
|
// roughly quadratic to linear).
|
|
const auto t = (y - fdisp) / (fmax - fdisp);
|
|
|
|
y = fr + t*(fm - fr);
|
|
}
|
|
else if (s < sm) {
|
|
// s \in [sr, sm), sm > sr; special case: Kr(Smax) == Kr(Sr).
|
|
//
|
|
// Use linear function between (sr,fr) and (sm,fm) in terms of
|
|
// saturation value 's'. This usually alters the shape of the
|
|
// relative permeability function in this interval (e.g.,
|
|
// roughly quadratic to linear).
|
|
const auto t = (s - sr) / (sm - sr);
|
|
|
|
y = fr + t*(fm - fr);
|
|
}
|
|
else {
|
|
// sm == sr (pure scaling). Almost arbitrarily pick fmax_[c].
|
|
y = fm;
|
|
}
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Class Opm::ThreePointScaling::Impl
|
|
// ---------------------------------------------------------------------
|
|
|
|
class Opm::SatFunc::ThreePointScaling::Impl
|
|
{
|
|
public:
|
|
Impl(std::vector<double> smin ,
|
|
std::vector<double> sdisp,
|
|
std::vector<double> smax )
|
|
: smin_(std::move(smin)), sdisp_(std::move(sdisp)), smax_(std::move(smax))
|
|
{
|
|
if ((this->sdisp_.size() != this->smin_.size()) ||
|
|
(this->sdisp_.size() != this->smax_.size()))
|
|
{
|
|
throw std::invalid_argument {
|
|
"Size Mismatch Between Minimum, Displacing "
|
|
"and Maximum Saturation Arrays"
|
|
};
|
|
}
|
|
}
|
|
|
|
std::vector<double>
|
|
eval(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const;
|
|
|
|
std::vector<double>
|
|
reverse(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const;
|
|
|
|
private:
|
|
std::vector<double> smin_;
|
|
std::vector<double> sdisp_;
|
|
std::vector<double> smax_;
|
|
|
|
double sMin(const std::vector<int>::size_type cell,
|
|
const TableEndPoints& tep) const
|
|
{
|
|
// Use model's scaled connate saturation if defined, otherwise fall
|
|
// back to table's unscaled connate saturation.
|
|
|
|
return defaultedScaledSaturation(this->smin_[cell], tep.low);
|
|
}
|
|
|
|
double sDisp(const std::vector<int>::size_type cell,
|
|
const TableEndPoints& tep) const
|
|
{
|
|
// Use model's scaled displacing saturation if defined, otherwise
|
|
// fall back to table's unscaled displacing saturation.
|
|
|
|
return defaultedScaledSaturation(this->sdisp_[cell], tep.disp);
|
|
}
|
|
|
|
double sMax(const std::vector<int>::size_type cell,
|
|
const TableEndPoints& tep) const
|
|
{
|
|
// Use model's scaled maximum saturation if defined, otherwise fall
|
|
// back to table's unscaled maximum saturation.
|
|
|
|
return defaultedScaledSaturation(this->smax_[cell], tep.high);
|
|
}
|
|
};
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::ThreePointScaling::
|
|
Impl::eval(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
auto effsat = std::vector<double>{};
|
|
effsat.reserve(sp.size());
|
|
|
|
for (const auto& eval_pt : sp) {
|
|
const auto cell = eval_pt.cell;
|
|
|
|
const auto sLO = this->sMin (cell, tep);
|
|
const auto sR = this->sDisp(cell, tep);
|
|
const auto sHI = this->sMax (cell, tep);
|
|
|
|
effsat.push_back(0.0);
|
|
auto& s_eff = effsat.back();
|
|
|
|
if (! (eval_pt.sat > sLO)) {
|
|
// s <= sLO
|
|
s_eff = tep.low;
|
|
}
|
|
else if (eval_pt.sat < std::min(sR, sHI)) {
|
|
// s in scaled interval [sLO, sR)
|
|
// Map to tabulated saturation in [tep.low, tep.disp)
|
|
const auto t = (eval_pt.sat - sLO) / (sR - sLO);
|
|
|
|
s_eff = tep.low + t*(tep.disp - tep.low);
|
|
}
|
|
else if (eval_pt.sat < sHI) {
|
|
// s in scaled interval [sR, sHI)
|
|
// Map to tabulated saturation in [tep.disp, tep.high)
|
|
assert (sHI > sR);
|
|
|
|
const auto t = (eval_pt.sat - sR) / (sHI - sR);
|
|
|
|
s_eff = tep.disp + t*(tep.high - tep.disp);
|
|
}
|
|
else {
|
|
// s >= sHI
|
|
s_eff = tep.high;
|
|
}
|
|
}
|
|
|
|
return effsat;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::ThreePointScaling::
|
|
Impl::reverse(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
auto unscaledsat = std::vector<double>{};
|
|
unscaledsat.reserve(sp.size());
|
|
|
|
for (const auto& eval_pt : sp) {
|
|
const auto cell = eval_pt.cell;
|
|
|
|
const auto sLO = this->sMin (cell, tep);
|
|
const auto sR = this->sDisp(cell, tep);
|
|
const auto sHI = this->sMax (cell, tep);
|
|
|
|
unscaledsat.push_back(0.0);
|
|
auto& s_unsc = unscaledsat.back();
|
|
|
|
if (! (eval_pt.sat > tep.low)) {
|
|
// s <= minimum tabulated saturation.
|
|
// Map to Minimum Input Saturation in cell (sLO).
|
|
s_unsc = sLO;
|
|
}
|
|
else if (eval_pt.sat < tep.disp) {
|
|
// s in tabulated interval [tep.low, tep.disp)
|
|
// Map to Input Saturation in (sLO, sR)
|
|
const auto t =
|
|
(eval_pt.sat - tep.low)
|
|
/ (tep.disp - tep.low);
|
|
|
|
s_unsc = std::min(sLO + t*(sR - sLO), sHI);
|
|
}
|
|
else if (eval_pt.sat < tep.high) {
|
|
// s in tabulated interval [tep.disp, tep.high)
|
|
// Map to Input Saturation in [sR, sHI)
|
|
assert (tep.high > tep.disp);
|
|
|
|
const auto t =
|
|
(eval_pt.sat - tep.disp)
|
|
/ (tep.high - tep.disp);
|
|
|
|
s_unsc = std::min(sR + t*std::max(sHI - sR, 0.0), sHI);
|
|
}
|
|
else {
|
|
// s >= maximum tabulated saturation.
|
|
//
|
|
// Map to Maximum Input Saturation in cell (sHI) if maximum
|
|
// tabulated saturation is strictly greater than critical
|
|
// displacing saturation--otherwise map to critical displacing
|
|
// saturation.
|
|
//
|
|
// Needed to handle cases in which \code tep.disp==tep.high
|
|
// \endcode but scaled versions of these might differ, i.e. when
|
|
// sR < sHI, but the corresponding saturation points in the
|
|
// underlying input table coincide.
|
|
s_unsc = (tep.high > tep.disp) ? sHI : sR;
|
|
}
|
|
}
|
|
|
|
return unscaledsat;
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// EPS factory functions for two-point and three-point scaling options
|
|
// ---------------------------------------------------------------------
|
|
|
|
namespace Create {
|
|
using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions;
|
|
using RTEP = ::Opm::SatFunc::CreateEPS::RawTableEndPoints;
|
|
using TEP = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints;
|
|
|
|
namespace TwoPoint {
|
|
using EPS = ::Opm::SatFunc::TwoPointScaling;
|
|
using EPSPtr = std::unique_ptr<EPS>;
|
|
|
|
struct Kr {
|
|
static EPSPtr
|
|
G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
OG(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
};
|
|
|
|
struct Pc {
|
|
static EPSPtr
|
|
GO(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
};
|
|
|
|
static EPSPtr
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const RTEP& tep);
|
|
|
|
static std::vector<TEP>
|
|
unscaledEndPoints(const RTEP& ep, const EPSOpt& opt);
|
|
} // namespace TwoPoint
|
|
|
|
namespace ThreePoint {
|
|
using EPS = ::Opm::SatFunc::ThreePointScaling;
|
|
using EPSPtr = std::unique_ptr<EPS>;
|
|
|
|
struct Kr {
|
|
static EPSPtr
|
|
G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
OG(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static EPSPtr
|
|
W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
};
|
|
|
|
static EPSPtr
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const RTEP& tep);
|
|
|
|
static std::vector<TEP>
|
|
unscaledEndPoints(const RTEP& ep, const EPSOpt& opt);
|
|
} // namespace ThreePoint
|
|
|
|
namespace PureVertical {
|
|
using Scaling = ::Opm::SatFunc::PureVerticalScaling;
|
|
using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions;
|
|
using FValVec = ::Opm::SatFunc::CreateEPS::Vertical::FuncValVector;
|
|
using ScalPtr = std::unique_ptr<Scaling>;
|
|
|
|
struct Kr {
|
|
static ScalPtr
|
|
G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt);
|
|
|
|
static ScalPtr
|
|
O(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt);
|
|
|
|
static ScalPtr
|
|
W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt);
|
|
};
|
|
|
|
struct Pc {
|
|
static ScalPtr
|
|
GO(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt);
|
|
|
|
static ScalPtr
|
|
OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt);
|
|
};
|
|
|
|
static ScalPtr
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const FValVec& fvals);
|
|
} // namespace PureVertical
|
|
|
|
namespace CritSatVertical {
|
|
using Scaling = ::Opm::SatFunc::CritSatVerticalScaling;
|
|
using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions;
|
|
using FValVec = ::Opm::SatFunc::CreateEPS::Vertical::FuncValVector;
|
|
using ScalPtr = std::unique_ptr<Scaling>;
|
|
|
|
namespace CritDispSat {
|
|
struct KrG {
|
|
static std::vector<double>
|
|
twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil);
|
|
|
|
static std::vector<double>
|
|
alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil);
|
|
};
|
|
|
|
struct KrGO {
|
|
static std::vector<double>
|
|
twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static std::vector<double>
|
|
alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
};
|
|
|
|
struct KrOW {
|
|
static std::vector<double>
|
|
twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
|
|
static std::vector<double>
|
|
alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep);
|
|
};
|
|
|
|
struct KrW {
|
|
static std::vector<double>
|
|
twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil);
|
|
|
|
static std::vector<double>
|
|
alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil);
|
|
};
|
|
|
|
std::vector<double>
|
|
transformedCritSat(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& t,
|
|
const std::vector<double>& left,
|
|
const std::vector<double>& right);
|
|
} // namespace CritDispSat
|
|
|
|
struct Kr {
|
|
static ScalPtr
|
|
G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval);
|
|
|
|
static ScalPtr
|
|
GO(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval);
|
|
|
|
static ScalPtr
|
|
OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval);
|
|
|
|
static ScalPtr
|
|
W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval);
|
|
};
|
|
|
|
static ScalPtr
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const RTEP& tep,
|
|
const FValVec& fvals);
|
|
} // namespace CritSatVertical
|
|
} // namespace Create
|
|
|
|
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
// Implementation of Create::TwoPoint::scalingFunction()
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::Kr::G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sgcr = sgCrit(G, init, tep);
|
|
auto sgu = sgMax (G, init, tep);
|
|
|
|
if ((sgcr.size() != sgu.size()) ||
|
|
(sgcr.size() != G.numCells()))
|
|
{
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Gas End-Point "
|
|
"Specifications (SGCR and/or SGU)"
|
|
};
|
|
}
|
|
|
|
return EPSPtr {
|
|
new EPS { std::move(sgcr), std::move(sgu) }
|
|
};
|
|
}
|
|
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::Kr::OG(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sogcr = sogCrit(G, init, tep);
|
|
|
|
if (sogcr.size() != G.numCells()) {
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Critical Oil "
|
|
"Saturation in Oil/Gas System"
|
|
};
|
|
}
|
|
|
|
auto smax = soMax(G, init, tep);
|
|
|
|
return EPSPtr {
|
|
new EPS { std::move(sogcr), std::move(smax) }
|
|
};
|
|
}
|
|
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::Kr::OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sowcr = sowCrit(G, init, tep);
|
|
|
|
if (sowcr.size() != G.numCells()) {
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Critical Oil "
|
|
"Saturation in Oil/Water System"
|
|
};
|
|
}
|
|
|
|
auto smax = soMax(G, init, tep);
|
|
|
|
return EPSPtr {
|
|
new EPS { std::move(sowcr), std::move(smax) }
|
|
};
|
|
}
|
|
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::Kr::W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto swcr = swCrit(G, init, tep);
|
|
auto swu = swMax (G, init, tep);
|
|
|
|
if (swcr.empty() || swu.empty()) {
|
|
throw std::invalid_argument {
|
|
"Missing Water End-Point Specifications (SWCR and/or SWU)"
|
|
};
|
|
}
|
|
|
|
return EPSPtr {
|
|
new EPS { std::move(swcr), std::move(swu) }
|
|
};
|
|
}
|
|
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
// Use dedicated scaled Sg_conn for Pc if defined in at least one
|
|
// subgrid. Use SGL otherwise. Same default value (i.e., table's
|
|
// connate gas saturation) for both vectors.
|
|
const auto sgl = ::Opm::SatFunc::scaledConnateGas(G, init, tep);
|
|
|
|
auto sglpc = haveKeywordData(G, init, "SGLPC")
|
|
? gridDefaultedVector(G, init, "SGLPC", tep.conn.gas,
|
|
[](const double s) { return s; }, sgl)
|
|
: sgl;
|
|
|
|
auto sgu = sgMax(G, init, tep);
|
|
|
|
if ((sgl.size() != sgu.size()) ||
|
|
(sgl.size() != G.numCells()))
|
|
{
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Connate or Maximum Gas "
|
|
"Saturation in Pcgo EPS"
|
|
};
|
|
}
|
|
|
|
return EPSPtr {
|
|
new EPS { std::move(sglpc), std::move(sgu) }
|
|
};
|
|
}
|
|
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
// Use dedicated scaled Sw_conn for Pc if defined in at least one
|
|
// subgrid. Use SWL otherwise. Same default value (i.e., table's
|
|
// connate water saturation) for both vectors.
|
|
const auto swl = ::Opm::SatFunc::scaledConnateWater(G, init, tep);
|
|
|
|
auto swlpc = haveKeywordData(G, init, "SWLPC")
|
|
? gridDefaultedVector(G, init, "SWLPC", tep.conn.water,
|
|
[](const double s) { return s; }, swl)
|
|
: swl;
|
|
|
|
auto swu = swMax(G, init, tep);
|
|
|
|
if ((swl.size() != swu.size()) ||
|
|
(swl.size() != G.numCells()))
|
|
{
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Connate or Maximum Water "
|
|
"Saturation in Pcow EPS"
|
|
};
|
|
}
|
|
|
|
return EPSPtr {
|
|
new EPS { std::move(swlpc), std::move(swu) }
|
|
};
|
|
}
|
|
|
|
Create::TwoPoint::EPSPtr
|
|
Create::TwoPoint::
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const RTEP& tep)
|
|
{
|
|
using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory;
|
|
using SSys = ::Opm::SatFunc::CreateEPS::SubSystem;
|
|
using PhIdx = ::Opm::ECLPhaseIndex;
|
|
|
|
assert (((! opt.use3PtScaling) || (opt.curve == FCat::CapPress))
|
|
&& "Internal Error Selecting EPS Family");
|
|
|
|
if (opt.curve == FCat::Relperm) {
|
|
if (opt.subSys == SSys::OilWater) {
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create an EPS for Gas Relperm "
|
|
"in an Oil/Water System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return Create::TwoPoint::Kr::W(G, init, tep);
|
|
}
|
|
|
|
return Create::TwoPoint::Kr::OW(G, init, tep);
|
|
}
|
|
|
|
if (opt.subSys == SSys::OilGas) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create an EPS for Water Relperm "
|
|
"in an Oil/Gas System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return Create::TwoPoint::Kr::G(G, init, tep);
|
|
}
|
|
|
|
return Create::TwoPoint::Kr::OG(G, init, tep);
|
|
}
|
|
}
|
|
|
|
if (opt.curve == FCat::CapPress) {
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
throw std::invalid_argument {
|
|
"Creating Capillary Pressure EPS as a Function "
|
|
"of Oil Saturation is not Supported"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return Create::TwoPoint::Pc::GO(G, init, tep);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return Create::TwoPoint::Pc::OW(G, init, tep);
|
|
}
|
|
}
|
|
|
|
// Invalid.
|
|
return EPSPtr{};
|
|
}
|
|
|
|
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
// Implementation of Create::TwoPoint::unscaledEndPoints()
|
|
std::vector<Create::TEP>
|
|
Create::TwoPoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt)
|
|
{
|
|
using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory;
|
|
using SSys = ::Opm::SatFunc::CreateEPS::SubSystem;
|
|
using PhIdx = ::Opm::ECLPhaseIndex;
|
|
|
|
assert (((opt.curve == FCat::CapPress) ||
|
|
(! opt.use3PtScaling)) && "Internal Logic Error");
|
|
|
|
if (opt.curve == FCat::CapPress) {
|
|
// Left node is connate saturation, right node is max saturation.
|
|
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
throw std::invalid_argument {
|
|
"No Capillary Pressure Function for Oil"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return unscaledTwoPt(ep.conn.water, ep.smax.water);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return unscaledTwoPt(ep.conn.gas, ep.smax.gas);
|
|
}
|
|
}
|
|
|
|
if (opt.curve == FCat::Relperm) {
|
|
// Left node is critical saturation, right node is max saturation.
|
|
|
|
if (opt.subSys == SSys::OilGas) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Void Request for Unscaled Water Saturation "
|
|
"End-Points in Oil-Gas System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
return unscaledTwoPt(ep.crit.oil_in_gas, ep.smax.oil);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return unscaledTwoPt(ep.crit.gas, ep.smax.gas);
|
|
}
|
|
}
|
|
|
|
if (opt.subSys == SSys::OilWater) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return unscaledTwoPt(ep.crit.water, ep.smax.water);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
return unscaledTwoPt(ep.crit.oil_in_water, ep.smax.oil);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Void Request for Unscaled Gas Saturation "
|
|
"End-Points in Oil-Water System"
|
|
};
|
|
}
|
|
}
|
|
}
|
|
|
|
// Invalid.
|
|
return {};
|
|
}
|
|
|
|
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
// Implementation of Create::ThreePoint::scalingFunction()
|
|
|
|
Create::ThreePoint::EPSPtr
|
|
Create::ThreePoint::Kr::G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sgcr = sgCrit(G, init, tep);
|
|
auto sgu = sgMax (G, init, tep);
|
|
|
|
if ((sgcr.size() != sgu.size()) ||
|
|
(sgcr.size() != G.numCells()))
|
|
{
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Gas End-Point "
|
|
"Specifications (SGCR and/or SGU)"
|
|
};
|
|
}
|
|
|
|
auto sdisp = CritSatVertical::CritDispSat::
|
|
KrG::alternateMethod(G, init, tep, true);
|
|
|
|
return EPSPtr {
|
|
new EPS {
|
|
std::move(sgcr), std::move(sdisp), std::move(sgu)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::ThreePoint::EPSPtr
|
|
Create::ThreePoint::Kr::OG(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sogcr = sogCrit(G, init, tep);
|
|
|
|
if (sogcr.size() != G.numCells()) {
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Critical Oil "
|
|
"Saturation in Oil/Gas System"
|
|
};
|
|
}
|
|
|
|
auto sdisp = CritSatVertical::CritDispSat::
|
|
KrGO::alternateMethod(G, init, tep);
|
|
|
|
auto smax = soMax(G, init, tep);
|
|
|
|
return EPSPtr {
|
|
new EPS {
|
|
std::move(sogcr), std::move(sdisp), std::move(smax)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::ThreePoint::EPSPtr
|
|
Create::ThreePoint::Kr::OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sowcr = sowCrit(G, init, tep);
|
|
|
|
if (sowcr.size() != G.numCells()) {
|
|
throw std::invalid_argument {
|
|
"Missing or Mismatching Critical Oil "
|
|
"Saturation in Oil/Water System"
|
|
};
|
|
}
|
|
|
|
auto sdisp = CritSatVertical::CritDispSat::
|
|
KrOW::alternateMethod(G, init, tep);
|
|
|
|
auto smax = soMax(G, init, tep);
|
|
|
|
return EPSPtr {
|
|
new EPS {
|
|
std::move(sowcr), std::move(sdisp), std::move(smax)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::ThreePoint::EPSPtr
|
|
Create::ThreePoint::Kr::W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto swcr = swCrit(G, init, tep);
|
|
auto swu = swMax (G, init, tep);
|
|
|
|
if ((swcr.size() != G.numCells()) ||
|
|
(swcr.size() != swu.size()))
|
|
{
|
|
throw std::invalid_argument {
|
|
"Missing Water End-Point Specifications (SWCR and/or SWU)"
|
|
};
|
|
}
|
|
|
|
auto sdisp = CritSatVertical::CritDispSat::
|
|
KrW::alternateMethod(G, init, tep, true);
|
|
|
|
return EPSPtr {
|
|
new EPS {
|
|
std::move(swcr), std::move(sdisp), std::move(swu)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::ThreePoint::EPSPtr
|
|
Create::ThreePoint::
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const RTEP& tep)
|
|
{
|
|
#if !defined(NDEBUG)
|
|
using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory;
|
|
#endif // !defined(NDEBUG)
|
|
|
|
using SSys = ::Opm::SatFunc::CreateEPS::SubSystem;
|
|
using PhIdx = ::Opm::ECLPhaseIndex;
|
|
|
|
assert ((opt.use3PtScaling && (opt.curve == FCat::Relperm))
|
|
&& "Internal Error Selecting EPS Family");
|
|
|
|
if (opt.subSys == SSys::OilWater) {
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create a Three-Point EPS for "
|
|
"Gas Relperm in an Oil/Water System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return Create::ThreePoint::Kr::W(G, init, tep);
|
|
}
|
|
|
|
return Create::ThreePoint::Kr::OW(G, init, tep);
|
|
}
|
|
|
|
if (opt.subSys == SSys::OilGas) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create a Three-Point EPS for "
|
|
"Water Relperm in an Oil/Gas System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return Create::ThreePoint::Kr::G(G, init, tep);
|
|
}
|
|
|
|
return Create::ThreePoint::Kr::OG(G, init, tep);
|
|
}
|
|
|
|
// Invalid.
|
|
return EPSPtr{};
|
|
}
|
|
|
|
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
// Implementation of Create::ThreePoint::unscaledEndPoints()
|
|
std::vector<Create::TEP>
|
|
Create::ThreePoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt)
|
|
{
|
|
#if !defined(NDEBUG)
|
|
using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory;
|
|
#endif // !defined(NDEBUG)
|
|
|
|
using SSys = ::Opm::SatFunc::CreateEPS::SubSystem;
|
|
using PhIdx = ::Opm::ECLPhaseIndex;
|
|
|
|
assert ((opt.use3PtScaling && (opt.curve == FCat::Relperm))
|
|
&& "Internal Error Selecting EPS Family");
|
|
|
|
auto sdisp = [](const std::vector<double>& s1,
|
|
const std::vector<double>& s2)
|
|
-> std::vector<double>
|
|
{
|
|
auto sr = std::vector<double>(s1.size(), 1.0);
|
|
|
|
for (auto n = s1.size(), i = 0*n; i < n; ++i) {
|
|
sr[i] -= s1[i] + s2[i];
|
|
}
|
|
|
|
return sr;
|
|
};
|
|
|
|
// Left node is critical saturation, middle node is displacing critical
|
|
// saturation, and right node is maximum saturation.
|
|
|
|
if (opt.subSys == SSys::OilGas) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Void Request for Unscaled Water Saturation "
|
|
"End-Points in Oil-Gas System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
return unscaledThreePt(ep.crit.oil_in_gas,
|
|
sdisp(ep.crit.gas, ep.conn.water),
|
|
ep.smax.oil);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return unscaledThreePt(ep.crit.gas,
|
|
sdisp(ep.crit.oil_in_gas, ep.conn.water),
|
|
ep.smax.gas);
|
|
}
|
|
}
|
|
|
|
if (opt.subSys == SSys::OilWater) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return unscaledThreePt(ep.crit.water,
|
|
sdisp(ep.crit.oil_in_water, ep.conn.gas),
|
|
ep.smax.water);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
return unscaledThreePt(ep.crit.oil_in_water,
|
|
sdisp(ep.crit.water, ep.conn.gas),
|
|
ep.smax.oil);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Void Request for Unscaled Gas Saturation "
|
|
"End-Points in Oil-Water System"
|
|
};
|
|
}
|
|
}
|
|
|
|
// Invalid.
|
|
return {};
|
|
}
|
|
|
|
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
// Implementation of Create::PureVertical::scalingFunction()
|
|
|
|
namespace {
|
|
Create::PureVertical::ScalPtr
|
|
pureVerticalRelpermScaling(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt,
|
|
const std::string& vector)
|
|
{
|
|
auto scaledMaxKr =
|
|
gridDefaultedVector(G, init, vector, dflt,
|
|
[](const double kr) { return kr; });
|
|
|
|
return Create::PureVertical::ScalPtr {
|
|
new ::Opm::SatFunc::PureVerticalScaling(std::move(scaledMaxKr))
|
|
};
|
|
}
|
|
|
|
Create::PureVertical::ScalPtr
|
|
pureVerticalCapPressScaling(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt,
|
|
const std::string& vector)
|
|
{
|
|
const auto& ih = init.keywordData<int>(INTEHEAD_KW);
|
|
|
|
const auto pscale = ::Opm::ECLUnits::
|
|
createUnitSystem(ih[ INTEHEAD_UNIT_INDEX ])->pressure();
|
|
|
|
auto scaledMaxPc =
|
|
gridDefaultedVector(G, init, vector, dflt,
|
|
[pscale](const double pc)
|
|
{
|
|
return ::ImportedOpm::unit::convert::from(pc, pscale);
|
|
});
|
|
|
|
return Create::PureVertical::ScalPtr {
|
|
new ::Opm::SatFunc::PureVerticalScaling(std::move(scaledMaxPc))
|
|
};
|
|
}
|
|
|
|
} // Anonymous
|
|
|
|
Create::PureVertical::ScalPtr
|
|
Create::PureVertical::Kr::G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt)
|
|
{
|
|
return pureVerticalRelpermScaling(G, init, dflt, "KRG");
|
|
}
|
|
|
|
Create::PureVertical::ScalPtr
|
|
Create::PureVertical::Kr::O(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt)
|
|
{
|
|
return pureVerticalRelpermScaling(G, init, dflt, "KRO");
|
|
}
|
|
|
|
Create::PureVertical::ScalPtr
|
|
Create::PureVertical::Kr::W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt)
|
|
{
|
|
return pureVerticalRelpermScaling(G, init, dflt, "KRW");
|
|
}
|
|
|
|
Create::PureVertical::ScalPtr
|
|
Create::PureVertical::Pc::GO(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt)
|
|
{
|
|
return pureVerticalCapPressScaling(G, init, dflt, "PCG");
|
|
}
|
|
|
|
Create::PureVertical::ScalPtr
|
|
Create::PureVertical::Pc::OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& dflt)
|
|
{
|
|
return pureVerticalCapPressScaling(G, init, dflt, "PCW");
|
|
}
|
|
|
|
Create::PureVertical::ScalPtr
|
|
Create::PureVertical::
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const FValVec& fvals)
|
|
{
|
|
using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory;
|
|
using SSys = ::Opm::SatFunc::CreateEPS::SubSystem;
|
|
using PhIdx = ::Opm::ECLPhaseIndex;
|
|
using FVal = ::Opm::SatFunc::VerticalScalingInterface::FunctionValues;
|
|
|
|
auto dfltVals = std::vector<double>(fvals.size(), 0.0);
|
|
std::transform(std::begin(fvals), std::end(fvals),
|
|
std::begin(dfltVals),
|
|
[](const FVal& fv)
|
|
{
|
|
return fv.max.val;
|
|
});
|
|
|
|
if (opt.curve == FCat::Relperm) {
|
|
if (opt.subSys == SSys::OilGas) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create Vertical Scaling for "
|
|
"Water Relperm in an Oil/Gas System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return Create::PureVertical::Kr::G(G, init, dfltVals);
|
|
}
|
|
|
|
return Create::PureVertical::Kr::O(G, init, dfltVals);
|
|
}
|
|
|
|
if (opt.subSys == SSys::OilWater) {
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create Vertical Scaling for "
|
|
"Gas Relperm in an Oil/Water System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return Create::PureVertical::Kr::W(G, init, dfltVals);
|
|
}
|
|
|
|
return Create::PureVertical::Kr::O(G, init, dfltVals);
|
|
}
|
|
}
|
|
|
|
if (opt.curve == FCat::CapPress) {
|
|
if (opt.thisPh == PhIdx::Liquid) {
|
|
throw std::invalid_argument {
|
|
"Creating Capillary Pressure Vertical Scaling "
|
|
"as a Function of Oil Saturation is not Supported"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return Create::PureVertical::Pc::GO(G, init, dfltVals);
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return Create::PureVertical::Pc::OW(G, init, dfltVals);
|
|
}
|
|
}
|
|
|
|
// Invalid.
|
|
return {};
|
|
}
|
|
|
|
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
// Implementation of Create::CritSatVertical::scalingFunction()
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrG::twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil)
|
|
{
|
|
const auto& sgcr = tep.crit.gas;
|
|
const auto& sgu = tep.smax.gas;
|
|
|
|
auto t = std::vector<double>(sgcr.size(), 0.0);
|
|
if (activeOil) {
|
|
// G/O or G/O/W system
|
|
for (auto n = sgcr.size(), i = 0*n; i < n; ++i) {
|
|
const auto sr = 1.0 - (tep.crit.oil_in_gas[i] +
|
|
tep.conn.water [i]);
|
|
|
|
t[i] = (sr - sgcr[i]) / (sgu[i] - sgcr[i]);
|
|
}
|
|
}
|
|
else {
|
|
// G/W system.
|
|
for (auto n = sgcr.size(), i = 0*n; i < n; ++i) {
|
|
const auto sr = 1.0 - tep.crit.water[i];
|
|
|
|
t[i] = (sr - sgcr[i]) / (sgu[i] - sgcr[i]);
|
|
}
|
|
}
|
|
|
|
return transformedCritSat(G, init, t,
|
|
sgCrit(G, init, tep),
|
|
sgMax (G, init, tep));
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrG::alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil)
|
|
{
|
|
auto sdisp = std::vector<double>(G.numCells(), 0.0);
|
|
|
|
if (activeOil) {
|
|
// G/O or G/O/W system.
|
|
const auto sogcr = sogCrit(G, init, tep);
|
|
const auto swl = ::Opm::SatFunc::scaledConnateWater(G, init, tep);
|
|
|
|
std::transform(std::begin(sogcr), std::end(sogcr), std::begin(swl),
|
|
std::begin(sdisp),
|
|
[](const double so, const double sw) -> double
|
|
{
|
|
return 1.0 - (so + sw);
|
|
});
|
|
}
|
|
else {
|
|
// G/W system.
|
|
const auto swcr = swCrit(G, init, tep);
|
|
|
|
std::transform(std::begin(swcr), std::end(swcr),
|
|
std::begin(sdisp),
|
|
[](const double sw) -> double
|
|
{
|
|
return 1.0 - sw;
|
|
});
|
|
}
|
|
|
|
return sdisp;
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrGO::twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
const auto& sogcr = tep.crit.oil_in_gas;
|
|
|
|
auto t = std::vector<double>(sogcr.size(), 0.0);
|
|
{
|
|
const auto& sgco = tep.conn.gas;
|
|
const auto& swco = tep.conn.water;
|
|
const auto& sgcr = tep.crit.gas;
|
|
|
|
for (auto n = sgcr.size(), i = 0*n; i < n; ++i) {
|
|
const auto sr = 1.0 - (sgcr[i] + swco[i]);
|
|
const auto smax = 1.0 - (sgco[i] + swco[i]); // >= sr
|
|
|
|
t[i] = (sr - sogcr[i]) / (smax - sogcr[i]);
|
|
}
|
|
}
|
|
|
|
return transformedCritSat(G, init, t,
|
|
sogCrit(G, init, tep),
|
|
soMax (G, init, tep));
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrGO::alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sdisp = std::vector<double>(G.numCells(), 0.0);
|
|
|
|
const auto sgcr = sgCrit(G, init, tep);
|
|
const auto swl = ::Opm::SatFunc::scaledConnateWater(G, init, tep);
|
|
|
|
std::transform(std::begin(sgcr), std::end(sgcr), std::begin(swl),
|
|
std::begin(sdisp),
|
|
[](const double sg, const double sw) -> double
|
|
{
|
|
return 1.0 - (sg + sw);
|
|
});
|
|
|
|
return sdisp;
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrOW::twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
const auto& sowcr = tep.crit.oil_in_water;
|
|
|
|
auto t = std::vector<double>(sowcr.size(), 0.0);
|
|
{
|
|
const auto& sgco = tep.conn.gas;
|
|
const auto& swco = tep.conn.water;
|
|
const auto& swcr = tep.crit.water;
|
|
|
|
for (auto n = swcr.size(), i = 0*n; i < n; ++i) {
|
|
const auto sr = 1.0 - (swcr[i] + sgco[i]);
|
|
const auto smax = 1.0 - (swco[i] + sgco[i]); // >= sr
|
|
|
|
t[i] = (sr - sowcr[i]) / (smax - sowcr[i]);
|
|
}
|
|
}
|
|
|
|
return transformedCritSat(G, init, t,
|
|
sowCrit(G, init, tep),
|
|
soMax (G, init, tep));
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrOW::alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep)
|
|
{
|
|
auto sdisp = std::vector<double>(G.numCells(), 0.0);
|
|
|
|
const auto swcr = swCrit(G, init, tep);
|
|
const auto sgl = ::Opm::SatFunc::scaledConnateGas(G, init, tep);
|
|
|
|
std::transform(std::begin(swcr), std::end(swcr), std::begin(sgl),
|
|
std::begin(sdisp),
|
|
[](const double sw, const double sg) -> double
|
|
{
|
|
return 1.0 - (sg + sw);
|
|
});
|
|
|
|
return sdisp;
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrW::twoPointMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil)
|
|
{
|
|
const auto& swcr = tep.crit.water;
|
|
const auto& swu = tep.smax.water;
|
|
|
|
auto t = std::vector<double>(swcr.size(), 0.0);
|
|
if (activeOil) {
|
|
// G/O or G/O/W system
|
|
for (auto n = swcr.size(), i = 0*n; i < n; ++i) {
|
|
const auto sr = 1.0 - (tep.crit.oil_in_water[i] +
|
|
tep.conn.gas [i]);
|
|
|
|
t[i] = (sr - swcr[i]) / (swu[i] - swcr[i]);
|
|
}
|
|
}
|
|
else {
|
|
// G/W system.
|
|
for (auto n = swcr.size(), i = 0*n; i < n; ++i) {
|
|
const auto sr = 1.0 - tep.crit.gas[i];
|
|
|
|
t[i] = (sr - swcr[i]) / (swu[i] - swcr[i]);
|
|
}
|
|
}
|
|
|
|
return transformedCritSat(G, init, t,
|
|
swCrit(G, init, tep),
|
|
swMax (G, init, tep));
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
KrW::alternateMethod(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const RTEP& tep,
|
|
const bool activeOil)
|
|
{
|
|
auto sdisp = std::vector<double>(G.numCells(), 0.0);
|
|
|
|
if (activeOil) {
|
|
// G/O or G/O/W system.
|
|
const auto sowcr = sowCrit(G, init, tep);
|
|
const auto sgl = ::Opm::SatFunc::scaledConnateGas(G, init, tep);
|
|
|
|
std::transform(std::begin(sowcr), std::end(sowcr), std::begin(sgl),
|
|
std::begin(sdisp),
|
|
[](const double so, const double sg) -> double
|
|
{
|
|
return 1.0 - (so + sg);
|
|
});
|
|
}
|
|
else {
|
|
// G/W system.
|
|
const auto sgcr = sgCrit(G, init, tep);
|
|
|
|
std::transform(std::begin(sgcr), std::end(sgcr),
|
|
std::begin(sdisp),
|
|
[](const double sg) -> double
|
|
{
|
|
return 1.0 - sg;
|
|
});
|
|
}
|
|
|
|
return sdisp;
|
|
}
|
|
|
|
std::vector<double>
|
|
Create::CritSatVertical::CritDispSat::
|
|
transformedCritSat(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const std::vector<double>& t,
|
|
const std::vector<double>& left,
|
|
const std::vector<double>& right)
|
|
{
|
|
auto sdisp = std::vector<double>(G.numCells(), 0.0);
|
|
|
|
auto cellID = std::vector<double>::size_type{0};
|
|
for (const auto& gridID : G.activeGrids()) {
|
|
const auto nc = G.numCells(gridID);
|
|
|
|
const auto& snum = init.haveKeywordData("SATNUM", gridID)
|
|
? G.rawLinearisedCellData<int>(init, "SATNUM", gridID)
|
|
: std::vector<int>(nc, 1);
|
|
|
|
for (auto c = 0*nc; c < nc; ++c, ++cellID) {
|
|
const auto x = t[snum[c] - 1];
|
|
|
|
sdisp[cellID] = (1.0 - x)*left[cellID] + x*right[cellID];
|
|
}
|
|
}
|
|
|
|
return sdisp;
|
|
}
|
|
|
|
namespace {
|
|
std::vector<double>
|
|
critDispSat(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::EPSOptions& opt,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& rtep)
|
|
{
|
|
namespace CDS = Create::CritSatVertical::CritDispSat;
|
|
|
|
if (opt.curve != ::Opm::SatFunc::CreateEPS::FunctionCategory::Relperm) {
|
|
return {};
|
|
}
|
|
|
|
const auto& ih = init.keywordData<int>(INTEHEAD_KW);
|
|
const auto activeOil =
|
|
(ih[INTEHEAD_PHASE_INDEX] & (1u << 0u)) != 0;
|
|
|
|
if (opt.subSys == ::Opm::SatFunc::CreateEPS::SubSystem::OilGas) {
|
|
if (opt.thisPh == ::Opm::ECLPhaseIndex::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Cannot request Critical Scaled Saturation "
|
|
"for water in Gas/Oil system"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == ::Opm::ECLPhaseIndex::Liquid) {
|
|
return !opt.use3PtScaling
|
|
? CDS::KrGO::twoPointMethod (G, init, rtep)
|
|
: CDS::KrGO::alternateMethod(G, init, rtep);
|
|
}
|
|
|
|
return !opt.use3PtScaling
|
|
? CDS::KrG::twoPointMethod (G, init, rtep, activeOil)
|
|
: CDS::KrG::alternateMethod(G, init, rtep, activeOil);
|
|
}
|
|
|
|
if (opt.subSys == ::Opm::SatFunc::CreateEPS::SubSystem::OilWater) {
|
|
if (opt.thisPh == ::Opm::ECLPhaseIndex::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Cannot request Critical Scaled Saturation "
|
|
"for gas in Oil/Water system"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == ::Opm::ECLPhaseIndex::Liquid) {
|
|
return !opt.use3PtScaling
|
|
? CDS::KrOW::twoPointMethod (G, init, rtep)
|
|
: CDS::KrOW::alternateMethod(G, init, rtep);
|
|
}
|
|
|
|
return !opt.use3PtScaling
|
|
? CDS::KrW::twoPointMethod (G, init, rtep, activeOil)
|
|
: CDS::KrW::alternateMethod(G, init, rtep, activeOil);
|
|
}
|
|
|
|
// Invalid
|
|
return {};
|
|
}
|
|
|
|
std::vector<double>
|
|
maximumSat(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const ::Opm::SatFunc::CreateEPS::EPSOptions& opt,
|
|
const ::Opm::SatFunc::CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
switch (opt.thisPh) {
|
|
case ::Opm::ECLPhaseIndex::Aqua: return swMax(G, init, tep);
|
|
case ::Opm::ECLPhaseIndex::Liquid: return soMax(G, init, tep);
|
|
case ::Opm::ECLPhaseIndex::Vapour: return sgMax(G, init, tep);
|
|
}
|
|
|
|
throw std::invalid_argument {
|
|
"Unsupported Phase Index"
|
|
};
|
|
}
|
|
}
|
|
|
|
Create::CritSatVertical::ScalPtr
|
|
Create::CritSatVertical::Kr::G(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval)
|
|
{
|
|
using FVal = ::Opm::SatFunc::VerticalScalingInterface::FunctionValues;
|
|
|
|
auto dflt_fdisp = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fdisp),
|
|
[](const FVal& fv) { return fv.disp.val;});
|
|
|
|
auto fdisp =
|
|
gridDefaultedVector(G, init, "KRGR", dflt_fdisp,
|
|
[](const double kr) { return kr; });
|
|
|
|
auto dflt_fmax = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fmax),
|
|
[](const FVal& fv) { return fv.max.val; });
|
|
|
|
auto fmax =
|
|
gridDefaultedVector(G, init, "KRG", dflt_fmax,
|
|
[](const double kr) { return kr; });
|
|
|
|
return ScalPtr {
|
|
new ::Opm::SatFunc::CritSatVerticalScaling {
|
|
std::move(sdisp), std::move(fdisp),
|
|
std::move(smax) , std::move(fmax)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::CritSatVertical::ScalPtr
|
|
Create::CritSatVertical::Kr::GO(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval)
|
|
{
|
|
using FVal = ::Opm::SatFunc::VerticalScalingInterface::FunctionValues;
|
|
|
|
auto dflt_fdisp = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fdisp),
|
|
[](const FVal& fv) { return fv.disp.val;});
|
|
|
|
auto fdisp =
|
|
gridDefaultedVector(G, init, "KRORG", dflt_fdisp,
|
|
[](const double kr) { return kr; });
|
|
|
|
auto dflt_fmax = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fmax),
|
|
[](const FVal& fv) { return fv.max.val; });
|
|
|
|
auto fmax =
|
|
gridDefaultedVector(G, init, "KRO", dflt_fmax,
|
|
[](const double kr) { return kr; });
|
|
|
|
return ScalPtr {
|
|
new ::Opm::SatFunc::CritSatVerticalScaling {
|
|
std::move(sdisp), std::move(fdisp),
|
|
std::move(smax) , std::move(fmax)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::CritSatVertical::ScalPtr
|
|
Create::CritSatVertical::Kr::OW(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval)
|
|
{
|
|
using FVal = ::Opm::SatFunc::VerticalScalingInterface::FunctionValues;
|
|
|
|
auto dflt_fdisp = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fdisp),
|
|
[](const FVal& fv) { return fv.disp.val;});
|
|
|
|
auto fdisp =
|
|
gridDefaultedVector(G, init, "KRORW", dflt_fdisp,
|
|
[](const double kr) { return kr; });
|
|
|
|
auto dflt_fmax = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fmax),
|
|
[](const FVal& fv) { return fv.max.val; });
|
|
|
|
auto fmax =
|
|
gridDefaultedVector(G, init, "KRO", dflt_fmax,
|
|
[](const double kr) { return kr; });
|
|
|
|
return ScalPtr {
|
|
new ::Opm::SatFunc::CritSatVerticalScaling {
|
|
std::move(sdisp), std::move(fdisp),
|
|
std::move(smax) , std::move(fmax)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::CritSatVertical::ScalPtr
|
|
Create::CritSatVertical::Kr::W(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
std::vector<double>&& sdisp,
|
|
std::vector<double>&& smax,
|
|
const FValVec& fval)
|
|
{
|
|
using FVal = ::Opm::SatFunc::VerticalScalingInterface::FunctionValues;
|
|
|
|
auto dflt_fdisp = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fdisp),
|
|
[](const FVal& fv) { return fv.disp.val;});
|
|
|
|
auto fdisp =
|
|
gridDefaultedVector(G, init, "KRWR", dflt_fdisp,
|
|
[](const double kr) { return kr; });
|
|
|
|
auto dflt_fmax = std::vector<double>(fval.size(), 0.0);
|
|
std::transform(std::begin(fval), std::end(fval),
|
|
std::begin(dflt_fmax),
|
|
[](const FVal& fv) { return fv.max.val; });
|
|
|
|
auto fmax =
|
|
gridDefaultedVector(G, init, "KRW", dflt_fmax,
|
|
[](const double kr) { return kr; });
|
|
|
|
return ScalPtr {
|
|
new ::Opm::SatFunc::CritSatVerticalScaling {
|
|
std::move(sdisp), std::move(fdisp),
|
|
std::move(smax) , std::move(fmax)
|
|
}
|
|
};
|
|
}
|
|
|
|
Create::CritSatVertical::ScalPtr
|
|
Create::CritSatVertical::
|
|
scalingFunction(const ::Opm::ECLGraph& G,
|
|
const ::Opm::ECLInitFileData& init,
|
|
const EPSOpt& opt,
|
|
const RTEP& tep,
|
|
const FValVec& fvals)
|
|
{
|
|
using SSys = ::Opm::SatFunc::CreateEPS::SubSystem;
|
|
using PhIdx = ::Opm::ECLPhaseIndex;
|
|
|
|
auto scr = critDispSat(G, init, opt, tep);
|
|
auto smax = maximumSat (G, init, opt, tep);
|
|
|
|
if (opt.subSys == SSys::OilWater) {
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create Critical Saturation Vertical "
|
|
"Scaling for Gas Relperm in an Oil/Water System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
return Create::CritSatVertical::
|
|
Kr::W(G, init, std::move(scr), std::move(smax), fvals);
|
|
}
|
|
|
|
return Create::CritSatVertical::
|
|
Kr::OW(G, init, std::move(scr), std::move(smax), fvals);
|
|
}
|
|
|
|
if (opt.subSys == SSys::OilGas) {
|
|
if (opt.thisPh == PhIdx::Aqua) {
|
|
throw std::invalid_argument {
|
|
"Cannot Create Critical Saturation Vertical "
|
|
"Scaling for Water Relperm in an Oil/Gas System"
|
|
};
|
|
}
|
|
|
|
if (opt.thisPh == PhIdx::Vapour) {
|
|
return Create::CritSatVertical::
|
|
Kr::G(G, init, std::move(scr), std::move(smax), fvals);
|
|
}
|
|
|
|
return Create::CritSatVertical::
|
|
Kr::GO(G, init, std::move(scr), std::move(smax), fvals);
|
|
}
|
|
|
|
// Invalid.
|
|
return {};
|
|
}
|
|
|
|
// #####################################################################
|
|
// =====================================================================
|
|
// Public Interface Below Separator
|
|
// =====================================================================
|
|
// #####################################################################
|
|
|
|
// Class Opm::SatFunc::EPSEvalInterface
|
|
Opm::SatFunc::EPSEvalInterface::~EPSEvalInterface()
|
|
{}
|
|
|
|
// ---------------------------------------------------------------------
|
|
|
|
// Class Opm::SatFunc::VerticalScalingInterface
|
|
Opm::SatFunc::VerticalScalingInterface::~VerticalScalingInterface()
|
|
{}
|
|
|
|
// ---------------------------------------------------------------------
|
|
|
|
// Class Opm::SatFunc::TwoPointScaling
|
|
Opm::SatFunc::TwoPointScaling::
|
|
TwoPointScaling(std::vector<double> smin,
|
|
std::vector<double> smax)
|
|
: pImpl_(new Impl(std::move(smin), std::move(smax)))
|
|
{}
|
|
|
|
Opm::SatFunc::TwoPointScaling::~TwoPointScaling()
|
|
{}
|
|
|
|
Opm::SatFunc::TwoPointScaling::
|
|
TwoPointScaling(const TwoPointScaling& rhs)
|
|
: pImpl_(new Impl(*rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::TwoPointScaling::
|
|
TwoPointScaling(TwoPointScaling&& rhs)
|
|
: pImpl_(std::move(rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::TwoPointScaling&
|
|
Opm::SatFunc::TwoPointScaling::operator=(const TwoPointScaling& rhs)
|
|
{
|
|
this->pImpl_.reset(new Impl(*rhs.pImpl_));
|
|
|
|
return *this;
|
|
}
|
|
|
|
Opm::SatFunc::TwoPointScaling&
|
|
Opm::SatFunc::TwoPointScaling::operator=(TwoPointScaling&& rhs)
|
|
{
|
|
this->pImpl_ = std::move(rhs.pImpl_);
|
|
|
|
return *this;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::TwoPointScaling::eval(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
return this->pImpl_->eval(tep, sp);
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::TwoPointScaling::reverse(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
return this->pImpl_->reverse(tep, sp);
|
|
}
|
|
|
|
std::unique_ptr<Opm::SatFunc::EPSEvalInterface>
|
|
Opm::SatFunc::TwoPointScaling::clone() const
|
|
{
|
|
return std::unique_ptr<TwoPointScaling>(new TwoPointScaling(*this));
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
|
|
// Class Opm::SatFunc::PureVerticalScaling
|
|
|
|
Opm::SatFunc::PureVerticalScaling::
|
|
PureVerticalScaling(std::vector<double> fmax)
|
|
: pImpl_(new Impl(std::move(fmax)))
|
|
{}
|
|
|
|
Opm::SatFunc::PureVerticalScaling::~PureVerticalScaling()
|
|
{}
|
|
|
|
Opm::SatFunc::PureVerticalScaling::
|
|
PureVerticalScaling(const PureVerticalScaling& rhs)
|
|
: pImpl_(new Impl(*rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::PureVerticalScaling::
|
|
PureVerticalScaling(PureVerticalScaling&& rhs)
|
|
: pImpl_(std::move(rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::PureVerticalScaling&
|
|
Opm::SatFunc::PureVerticalScaling::operator=(const PureVerticalScaling& rhs)
|
|
{
|
|
this->pImpl_.reset(new Impl(*rhs.pImpl_));
|
|
|
|
return *this;
|
|
}
|
|
|
|
Opm::SatFunc::PureVerticalScaling&
|
|
Opm::SatFunc::PureVerticalScaling::operator=(PureVerticalScaling&& rhs)
|
|
{
|
|
this->pImpl_ = std::move(rhs.pImpl_);
|
|
|
|
return *this;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::PureVerticalScaling::
|
|
vertScale(const FunctionValues& f,
|
|
const SaturationPoints& sp,
|
|
const std::vector<double>& val) const
|
|
{
|
|
return this->pImpl_->vertScale(f, sp, val);
|
|
}
|
|
|
|
std::unique_ptr<Opm::SatFunc::VerticalScalingInterface>
|
|
Opm::SatFunc::PureVerticalScaling::clone() const
|
|
{
|
|
return std::unique_ptr<PureVerticalScaling>(new PureVerticalScaling(*this));
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
|
|
// Class Opm::SatFunc::ThreePointScaling
|
|
Opm::SatFunc::ThreePointScaling::
|
|
ThreePointScaling(std::vector<double> smin,
|
|
std::vector<double> sdisp,
|
|
std::vector<double> smax)
|
|
: pImpl_(new Impl(std::move(smin), std::move(sdisp), std::move(smax)))
|
|
{}
|
|
|
|
Opm::SatFunc::ThreePointScaling::~ThreePointScaling()
|
|
{}
|
|
|
|
Opm::SatFunc::ThreePointScaling::
|
|
ThreePointScaling(const ThreePointScaling& rhs)
|
|
: pImpl_(new Impl(*rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::ThreePointScaling::ThreePointScaling(ThreePointScaling&& rhs)
|
|
: pImpl_(std::move(rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::ThreePointScaling&
|
|
Opm::SatFunc::ThreePointScaling::operator=(const ThreePointScaling& rhs)
|
|
{
|
|
this->pImpl_.reset(new Impl(*rhs.pImpl_));
|
|
|
|
return *this;
|
|
}
|
|
|
|
Opm::SatFunc::ThreePointScaling&
|
|
Opm::SatFunc::ThreePointScaling::operator=(ThreePointScaling&& rhs)
|
|
{
|
|
this->pImpl_ = std::move(rhs.pImpl_);
|
|
|
|
return *this;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::ThreePointScaling::eval(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
return this->pImpl_->eval(tep, sp);
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::ThreePointScaling::reverse(const TableEndPoints& tep,
|
|
const SaturationPoints& sp) const
|
|
{
|
|
return this->pImpl_->reverse(tep, sp);
|
|
}
|
|
|
|
std::unique_ptr<Opm::SatFunc::EPSEvalInterface>
|
|
Opm::SatFunc::ThreePointScaling::clone() const
|
|
{
|
|
return std::unique_ptr<ThreePointScaling>(new ThreePointScaling(*this));
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
|
|
// Class Opm::SatFunc::CritSatVerticalScaling
|
|
Opm::SatFunc::CritSatVerticalScaling::
|
|
CritSatVerticalScaling(std::vector<double> sdisp,
|
|
std::vector<double> fdisp,
|
|
std::vector<double> smax,
|
|
std::vector<double> fmax)
|
|
: pImpl_(new Impl(std::move(sdisp), std::move(fdisp),
|
|
std::move(smax) , std::move(fmax)))
|
|
{}
|
|
|
|
Opm::SatFunc::CritSatVerticalScaling::~CritSatVerticalScaling()
|
|
{}
|
|
|
|
Opm::SatFunc::CritSatVerticalScaling::
|
|
CritSatVerticalScaling(const CritSatVerticalScaling& rhs)
|
|
: pImpl_(new Impl(*rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::CritSatVerticalScaling::
|
|
CritSatVerticalScaling(CritSatVerticalScaling&& rhs)
|
|
: pImpl_(std::move(rhs.pImpl_))
|
|
{}
|
|
|
|
Opm::SatFunc::CritSatVerticalScaling&
|
|
Opm::SatFunc::CritSatVerticalScaling::
|
|
operator=(const CritSatVerticalScaling& rhs)
|
|
{
|
|
this->pImpl_.reset(new Impl(*rhs.pImpl_));
|
|
|
|
return *this;
|
|
}
|
|
|
|
Opm::SatFunc::CritSatVerticalScaling&
|
|
Opm::SatFunc::CritSatVerticalScaling::
|
|
operator=(CritSatVerticalScaling&& rhs)
|
|
{
|
|
this->pImpl_ = std::move(rhs.pImpl_);
|
|
|
|
return *this;
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::CritSatVerticalScaling::
|
|
vertScale(const FunctionValues& f,
|
|
const SaturationPoints& sp,
|
|
const std::vector<double>& val) const
|
|
{
|
|
return this->pImpl_->vertScale(f, sp, val);
|
|
}
|
|
|
|
std::unique_ptr<Opm::SatFunc::VerticalScalingInterface>
|
|
Opm::SatFunc::CritSatVerticalScaling::clone() const
|
|
{
|
|
return std::unique_ptr<CritSatVerticalScaling> {
|
|
new CritSatVerticalScaling(*this)
|
|
};
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Factory function Opm::SatFunc::CreateEPS::Horizontal::fromECLOutput()
|
|
|
|
std::unique_ptr<Opm::SatFunc::EPSEvalInterface>
|
|
Opm::SatFunc::CreateEPS::Horizontal::
|
|
fromECLOutput(const ECLGraph& G,
|
|
const ECLInitFileData& init,
|
|
const EPSOptions& opt,
|
|
const RawTableEndPoints& tep)
|
|
{
|
|
if ((opt.curve == FunctionCategory::CapPress) ||
|
|
(! opt.use3PtScaling))
|
|
{
|
|
return Create::TwoPoint::scalingFunction(G, init, opt, tep);
|
|
}
|
|
|
|
if ((opt.curve == FunctionCategory::Relperm) && opt.use3PtScaling)
|
|
{
|
|
return Create::ThreePoint::scalingFunction(G, init, opt, tep);
|
|
}
|
|
|
|
// Invalid
|
|
return std::unique_ptr<Opm::SatFunc::EPSEvalInterface>{};
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Factory function Opm::SatFunc::CreateEPS::Horizontal::unscaledEndPoints()
|
|
|
|
std::vector<Opm::SatFunc::EPSEvalInterface::TableEndPoints>
|
|
Opm::SatFunc::CreateEPS::Horizontal::
|
|
unscaledEndPoints(const EPSOptions& opt,
|
|
const RawTableEndPoints& ep)
|
|
{
|
|
if ((opt.curve == FunctionCategory::CapPress) ||
|
|
(! opt.use3PtScaling))
|
|
{
|
|
return Create::TwoPoint::unscaledEndPoints(ep, opt);
|
|
}
|
|
|
|
if ((opt.curve == FunctionCategory::Relperm) && opt.use3PtScaling)
|
|
{
|
|
return Create::ThreePoint::unscaledEndPoints(ep, opt);
|
|
}
|
|
|
|
// Invalid
|
|
return {};
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Factory function Opm::SatFunc::CreateEPS::Vertical::fromECLOutput()
|
|
|
|
std::unique_ptr<Opm::SatFunc::VerticalScalingInterface>
|
|
Opm::SatFunc::CreateEPS::Vertical::
|
|
fromECLOutput(const ECLGraph& G,
|
|
const ECLInitFileData& init,
|
|
const EPSOptions& opt,
|
|
const RawTableEndPoints& tep,
|
|
const FuncValVector& fvals)
|
|
{
|
|
const auto haveScaleCRS = haveScaledRelPermAtCritSat(G, init, opt);
|
|
|
|
if ((opt.curve == FunctionCategory::CapPress) || (! haveScaleCRS))
|
|
{
|
|
return Create::PureVertical::
|
|
scalingFunction(G, init, opt, fvals);
|
|
}
|
|
|
|
if ((opt.curve == FunctionCategory::Relperm) && haveScaleCRS)
|
|
{
|
|
return Create::CritSatVertical::
|
|
scalingFunction(G, init, opt, tep, fvals);
|
|
}
|
|
|
|
// Invalid
|
|
return {};
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Factory function Opm::SatFunc::CreateEPS::Vertical::unscaledFunctionValues()
|
|
|
|
std::vector<Opm::SatFunc::VerticalScalingInterface::FunctionValues>
|
|
Opm::SatFunc::CreateEPS::Vertical::
|
|
unscaledFunctionValues(const ECLGraph& G,
|
|
const ECLInitFileData& init,
|
|
const RawTableEndPoints& ep,
|
|
const EPSOptions& opt,
|
|
const SatFuncEvaluator& evalSF)
|
|
{
|
|
auto ret = std::vector<VerticalScalingInterface::FunctionValues>{};
|
|
|
|
const auto haveScaleCRS = haveScaledRelPermAtCritSat(G, init, opt);
|
|
|
|
if ((opt.curve == FunctionCategory::CapPress) || (! haveScaleCRS)) {
|
|
auto opt_cpy = opt;
|
|
opt_cpy.use3PtScaling = false;
|
|
|
|
const auto uep =
|
|
Create::TwoPoint::unscaledEndPoints(ep, opt_cpy);
|
|
|
|
ret.resize(uep.size());
|
|
|
|
for (auto n = uep.size(), i = 0*n; i < n; ++i) {
|
|
ret[i].disp.sat = uep[i].disp;
|
|
ret[i].disp.val = evalSF(static_cast<int>(i), ret[i].disp.sat);
|
|
|
|
ret[i].max.sat = uep[i].high;
|
|
ret[i].max.val = evalSF(static_cast<int>(i), ret[i].max.sat);
|
|
}
|
|
}
|
|
else {
|
|
auto opt_cpy = opt;
|
|
opt_cpy.use3PtScaling = true;
|
|
|
|
const auto uep =
|
|
Create::ThreePoint::unscaledEndPoints(ep, opt_cpy);
|
|
|
|
ret.resize(uep.size());
|
|
|
|
for (auto n = uep.size(), i = 0*n; i < n; ++i) {
|
|
ret[i].disp.sat = uep[i].disp;
|
|
ret[i].disp.val = evalSF(static_cast<int>(i), ret[i].disp.sat);
|
|
|
|
ret[i].max.sat = uep[i].high;
|
|
ret[i].max.val = evalSF(static_cast<int>(i), ret[i].max.sat);
|
|
}
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Factory functions Opm::SatFunc::scaledConnate*()
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::scaledConnateGas(const ECLGraph& G,
|
|
const ECLInitFileData& init,
|
|
const CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SGL", tep.conn.gas,
|
|
[](const double s) { return s; });
|
|
}
|
|
|
|
std::vector<double>
|
|
Opm::SatFunc::scaledConnateWater(const ECLGraph& G,
|
|
const ECLInitFileData& init,
|
|
const CreateEPS::RawTableEndPoints& tep)
|
|
{
|
|
return gridDefaultedVector(G, init, "SWL", tep.conn.water,
|
|
[](const double s) { return s; });
|
|
}
|