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

LiveOil, WetGas: improve the saturationPressure() methods

Browse Source 
we now ask the tabulation class for the derivative directly and we do
not return negative saturation pressures anymore. (the latter is done
via putting the Newton method on "probation" if the saturation
pressure becomes negative. The first time it does this offense we try
to help it out of its mess, but if it happens again, we abort it and
return 0 as the saturation pressure.)
This commit is contained in:
Andreas Lauser 2017-03-09 17:14:28 +01:00
parent 777f8f0bc8
commit 00659e5847
3 changed files with 47 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -482,10 +482,10 @@ public:
*/
template <class Evaluation>
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
const Evaluation& /*temperature*/,
const Evaluation& pressure,
const Evaluation& oilSaturation,
Scalar maxOilSaturation) const
const Evaluation& /*temperature*/,
const Evaluation& pressure,
const Evaluation& oilSaturation,
Scalar maxOilSaturation) const
{
typedef typename Opm::MathToolbox<Evaluation> Toolbox;
Evaluation tmp =
@ -510,32 +510,44 @@ public:
*/
template <class Evaluation>
Evaluation saturationPressure(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& Rs) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
const auto& RsTable = saturatedGasDissolutionFactorTable_[regionIdx];
static constexpr Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
// use the saturation pressure function to get a pretty good initial value
Evaluation pSat = saturationPressure_[regionIdx].eval(Rs, /*extrapolate=*/true);
Evaluation eps = pSat*1e-11;
// Newton method to do the remaining work. If the initial
// value is good, this should only take two to three
// iterations...
bool onProbation = false;
for (int i = 0; i < 20; ++i) {
const Evaluation& f = saturatedGasDissolutionFactor(regionIdx, temperature, pSat) - Rs;
const Evaluation& fPrime = ((saturatedGasDissolutionFactor(regionIdx, temperature, pSat + eps) - Rs) - f)/eps;
const Evaluation& f = RsTable.eval(pSat, /*extrapolate=*/true) - Rs;
const Evaluation& fPrime = RsTable.evalDerivative(pSat, /*extrapolate=*/true);
const Evaluation& delta = f/fPrime;
pSat -= delta;
Scalar absDelta = std::abs(Toolbox::scalarValue(delta));
if (absDelta < Toolbox::scalarValue(pSat) * 1e-10 || absDelta < 1e-4)
if (pSat < 0.0) {
// if the pressure is lower than 0 Pascals, we set it back to 0. if this
// happens twice, we give up and just return 0 Pa...
if (onProbation)
return 0.0;
onProbation = true;
pSat = 0.0;
}
if (std::abs(Toolbox::scalarValue(delta)) < std::abs(Toolbox::scalarValue(pSat))*eps)
return pSat;
}
std::stringstream errlog;
errlog << "Could find the oil saturation pressure for X_o^G = " << Rs;
errlog << "Could find the oil saturation pressure for Rs = " << Rs;
OpmLog::problem("liveoil saturationpressure", errlog.str());
OPM_THROW_NOLOG(NumericalProblem, errlog.str());
}
@ -547,7 +559,7 @@ private:
// create the function representing saturation pressure depending of the mass
// fraction in gas
size_t n = gasDissolutionFac.numSamples()*5;
size_t n = gasDissolutionFac.numSamples();
Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/(n + 1);
SamplingPoints pSatSamplePoints;

4
opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp
View File

@ -219,8 +219,8 @@ public:
*/
template <class Evaluation>
Evaluation saturationPressure(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& Rs) const
const Evaluation& temperature,
const Evaluation& Rs) const
{ OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.saturationPressure(regionIdx, temperature, Rs)); return 0; }
void setApproach(OilPvtApproach appr)

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

@ -538,31 +538,44 @@ public:
*/
template <class Evaluation>
Evaluation saturationPressure(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& Rv) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
static constexpr Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
// use the tabulated saturation pressure function to get a pretty good initial value
Evaluation pSat = saturationPressure_[regionIdx].eval(Rv, /*extrapolate=*/true);
const Evaluation& eps = pSat*1e-11;
// Newton method to do the remaining work. If the initial
// value is good, this should only take two to three
// iterations...
bool onProbation = false;
for (unsigned i = 0; i < 20; ++i) {
const Evaluation& f = saturatedOilVaporizationFactor(regionIdx, temperature, pSat) - Rv;
const Evaluation& fPrime = ((saturatedOilVaporizationFactor(regionIdx, temperature, pSat + eps) - Rv) - f)/eps;
const Evaluation& f = RvTable.eval(pSat, /*extrapolate=*/true) - Rv;
const Evaluation& fPrime = RvTable.evalDerivative(pSat, /*extrapolate=*/true);
const Evaluation& delta = f/fPrime;
pSat -= delta;
if (std::abs(Toolbox::scalarValue(delta)) < std::abs(Toolbox::scalarValue(pSat)) * 1e-10)
if (pSat < 0.0) {
// if the pressure is lower than 0 Pascals, we set it back to 0. if this
// happens twice, we give up and just return 0 Pa...
if (onProbation)
return 0.0;
onProbation = true;
pSat = 0.0;
}
if (std::abs(Toolbox::scalarValue(delta)) < std::abs(Toolbox::scalarValue(pSat))*eps)
return pSat;
}
std::stringstream errlog;
errlog << "Could find the gas saturation pressure for X_g^O = " << Rv;
errlog << "Could find the gas saturation pressure for Rv = " << Rv;
OpmLog::problem("wetgas saturationpressure", errlog.str());
OPM_THROW_NOLOG(NumericalProblem, errlog.str());
}
@ -574,14 +587,14 @@ private:
// create the taublated function representing saturation pressure depending of
// Rv
size_t n = oilVaporizationFac.numSamples()*5;
size_t n = oilVaporizationFac.numSamples();
Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/(n + 1);
SamplingPoints pSatSamplePoints;
Scalar Rv = 0;
for (size_t i = 0; i <= n; ++ i) {
Scalar pSat = oilVaporizationFac.xMin() + i*delta;
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e100), pSat);
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
std::pair<Scalar, Scalar> val(Rv, pSat);
pSatSamplePoints.push_back(val);