Fixed numerical problems with bubble/dew point pressures

This commit is contained in:
babrodtk 2017-03-23 11:34:22 +01:00
parent 943f022277
commit d0b2eb9606
2 changed files with 39 additions and 25 deletions

View File

@ -529,16 +529,14 @@ public:
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());
// If the derivative is "zero" Newton will not converge,
// so simply return our initial guess.
if (std::abs(Toolbox::scalarValue(fPrime)) < 1.0e-30) {
return pSat;
}
const Evaluation& delta = f/fPrime;
pSat -= delta;
if (pSat < 0.0) {
@ -556,15 +554,18 @@ public:
}
std::stringstream errlog;
errlog << "Could find the oil saturation pressure for Rs = " << Rs;
OpmLog::problem("liveoil saturationpressure", errlog.str());
errlog << "Finding saturation pressure did not converge:"
<< " pSat = " << pSat
<< ", Rs = " << Rs;
OpmLog::problem("liveoil psat", errlog.str());
OPM_THROW_NOLOG(NumericalProblem, errlog.str());
}
private:
void updateSaturationPressure_(unsigned regionIdx)
{
auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
typedef std::pair<Scalar, Scalar> Pair;
const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
// create the function representing saturation pressure depending of the mass
// fraction in gas
@ -579,9 +580,15 @@ private:
/*temperature=*/Scalar(1e30),
pSat);
std::pair<Scalar, Scalar> val(Rs, pSat);
Pair val(Rs, pSat);
pSatSamplePoints.push_back(val);
}
//Prune duplicate Rs values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}

View File

@ -277,7 +277,7 @@ public:
{
auto& invGasB = inverseGasB_[regionIdx];
auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
Scalar T = 273.15 + 15.56; // [K]
@ -557,16 +557,14 @@ public:
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());
// If the derivative is "zero" Newton will not converge,
// so simply return our initial guess.
if (std::abs(Toolbox::scalarValue(fPrime)) < 1.0e-30) {
return pSat;
}
const Evaluation& delta = f/fPrime;
pSat -= delta;
if (pSat < 0.0) {
@ -584,15 +582,18 @@ public:
}
std::stringstream errlog;
errlog << "Could find the gas saturation pressure for Rv = " << Rv;
OpmLog::problem("wetgas saturationpressure", errlog.str());
errlog << "Finding saturation pressure did not converge:"
<< " pSat = " << pSat
<< ", Rv = " << Rv;
OpmLog::problem("wetgas psat", errlog.str());
OPM_THROW_NOLOG(NumericalProblem, errlog.str());
}
private:
void updateSaturationPressure_(unsigned regionIdx)
{
auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
typedef std::pair<Scalar, Scalar> Pair;
const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
// create the taublated function representing saturation pressure depending of
// Rv
@ -605,9 +606,15 @@ private:
Scalar pSat = oilVaporizationFac.xMin() + i*delta;
Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
std::pair<Scalar, Scalar> val(Rv, pSat);
Pair val(Rv, pSat);
pSatSamplePoints.push_back(val);
}
//Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
pSatSamplePoints.erase(last, pSatSamplePoints.end());
saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
}