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

Merge pull request #212 from andlaus/satpressure_improvements

Browse Source 
Saturation pressure & tabulation improvements
This commit is contained in:
Andreas Lauser 2017-03-09 21:30:42 +01:00 committed by GitHub
commit fc9885229c
5 changed files with 133 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
opm/material/common/Spline.hpp
View File

@ -30,6 +30,7 @@
#include <opm/material/common/TridiagonalMatrix.hpp>
#include <opm/material/common/PolynomialUtils.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/Unused.hpp>
#include <ostream>
@ -798,7 +799,11 @@ public:
template <class Evaluation>
Evaluation eval(const Evaluation& x, bool extrapolate = false) const
{
assert(extrapolate || applies(x));
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a spline outside of its range");
// handle extrapolation
if (extrapolate) {
@ -815,7 +820,7 @@ public:
}
}
return eval_(x, segmentIdx_(x));
return eval_(x, segmentIdx_(Toolbox::scalarValue(x)));
}
/*!
@ -833,7 +838,13 @@ public:
template <class Evaluation>
Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
{
assert(extrapolate || applies(x));
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate the derivative of a spline outside of its range");
// handle extrapolation
if (extrapolate) {
if (x < xAt(0))
return evalDerivative_(xAt(0), /*segmentIdx=*/0);
@ -841,7 +852,7 @@ public:
return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
}
return evalDerivative_(x, segmentIdx_(x));
return evalDerivative_(x, segmentIdx_(Toolbox::scalarValue(x)));
}
/*!
@ -859,11 +870,15 @@ public:
template <class Evaluation>
Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
{
assert(extrapolate || applies(x));
if (extrapolate)
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate the second derivative of a spline outside of its range");
else if (extrapolate)
return 0.0;
return evalDerivative2_(x, segmentIdx_(x));
return evalDerivative2_(x, segmentIdx_(Toolbox::scalarValue(x)));
}
/*!
@ -881,11 +896,15 @@ public:
template <class Evaluation>
Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
{
assert(extrapolate || applies(x));
if (extrapolate)
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate the third derivative of a spline outside of its range");
else if (extrapolate)
return 0.0;
return evalDerivative3_(x, segmentIdx_(x));
return evalDerivative3_(x, segmentIdx_(Toolbox::scalarValue(x)));
}
/*!
@ -1730,13 +1749,8 @@ protected:
}
// find the segment index for a given x coordinate
template <class Evaluation>
size_t segmentIdx_(const Evaluation& xEval) const
size_t segmentIdx_(Scalar x) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
Scalar x = Toolbox::scalarValue(xEval);
// bisection
size_t iLow = 0;
size_t iHigh = numSamples() - 1;

45
opm/material/common/Tabulated1DFunction.hpp
View File

@ -27,10 +27,10 @@
#ifndef OPM_TABULATED_1D_FUNCTION_HPP
#define OPM_TABULATED_1D_FUNCTION_HPP
#include <opm/material/densead/Math.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/Unused.hpp>
#include <opm/material/densead/Math.hpp>
#include <algorithm>
#include <cassert>
@ -235,7 +235,8 @@ public:
/*!
* \brief Return true iff the given x is in range [x1, xn].
*/
bool applies(Scalar x) const
template <class Evaluation>
bool applies(const Evaluation& x) const
{ return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
/*!
@ -252,15 +253,18 @@ public:
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a tabulated function outside of its range");
size_t segIdx;
if (extrapolate && x < xValues_.front())
segIdx = 0;
else if (extrapolate && x > xValues_.back())
segIdx = numSamples() - 2;
else {
assert(xValues_.front() <= x && x <= xValues_.back());
else
segIdx = findSegmentIndex_(Toolbox::scalarValue(x));
}
Scalar x0 = xValues_[segIdx];
Scalar x1 = xValues_[segIdx + 1];
@ -283,11 +287,18 @@ public:
* cause a failed assertation.
*/
template <class Evaluation>
Evaluation evalDerivative(const Evaluation& x, bool extrapolate OPM_UNUSED = false) const
Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
{
unsigned segIdx = findSegmentIndex_(x);
typedef Opm::MathToolbox<Evaluation> Toolbox;
return evalDerivative(x, segIdx);
if (!extrapolate && !applies(x)) {
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a derivative of a tabulated"
" function outside of its range");
}
unsigned segIdx = findSegmentIndex_(Toolbox::scalarValue(x));
return evalDerivative_(x, segIdx);
}
/*!
@ -305,9 +316,14 @@ public:
* cause a failed assertation.
*/
template <class Evaluation>
Evaluation evalSecondDerivative(const Evaluation OPM_OPTIM_UNUSED& x, bool extrapolate OPM_OPTIM_UNUSED = false) const
Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
{
assert(extrapolate || applies(x));
if (!extrapolate && !applies(x)) {
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a second derivative of a tabulated "
" function outside of its range");
}
return 0.0;
}
@ -326,9 +342,14 @@ public:
* cause a failed assertation.
*/
template <class Evaluation>
Evaluation evalThirdDerivative(const Evaluation OPM_OPTIM_UNUSED& x, bool extrapolate OPM_OPTIM_UNUSED = false) const
Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
{
assert(extrapolate || applies(x));
if (!extrapolate && !applies(x)) {
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a third derivative of a tabulated "
" function outside of its range");
}
return 0.0;
}

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

@ -125,9 +125,8 @@ public:
// update the tables for the formation volume factor and for the gas
// dissolution factor of saturated oil
{
std::vector<double> tmpPressureColumn = saturatedTable.getColumn("P").vectorCopy();
std::vector<double> tmpGasSolubilityColumn = saturatedTable.getColumn("RS").vectorCopy();
std::vector<double> tmpMuColumn = saturatedTable.getColumn("MU").vectorCopy();
const auto& tmpPressureColumn = saturatedTable.getColumn("P");
const auto& tmpGasSolubilityColumn = saturatedTable.getColumn("RS");
invSatOilB.setXYContainers(tmpPressureColumn, invSatOilBArray);
satOilMu.setXYContainers(tmpPressureColumn, satOilMuArray);
@ -483,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 =
@ -511,32 +510,53 @@ 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;
if (!Opm::isfinite(delta)) {
std::stringstream errlog;
errlog << "Obtain delta with non-finite value for Rs = " << Rs << " in the "
<< i << "th iteration during finding saturation pressure iteratively";
OpmLog::problem("wetgas NaN delta value", errlog.str());
OPM_THROW_NOLOG(NumericalProblem, errlog.str());
}
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());
}
@ -548,7 +568,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)

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

@ -92,12 +92,9 @@ public:
auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
{
std::vector<double> pressure = saturatedTable.getColumn("PG").vectorCopy( );
std::vector<double> rv = saturatedTable.getColumn("RV").vectorCopy( );
oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
pressure , rv );
}
oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
saturatedTable.getColumn("PG"),
saturatedTable.getColumn("RV"));
std::vector<Scalar> invSatGasBArray;
std::vector<Scalar> invSatGasBMuArray;
@ -541,31 +538,53 @@ 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;
if (!Opm::isfinite(delta)) {
std::stringstream errlog;
errlog << "Obtain delta with non-finite value for Rv = " << Rv << " in the "
<< i << "th iteration during finding saturation pressure iteratively";
OpmLog::problem("wetgas NaN delta value", errlog.str());
OPM_THROW_NOLOG(NumericalProblem, errlog.str());
}
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());
}
@ -577,14 +596,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);