Merge pull request #328 from atgeirr/interpolation_2d_miscibility_fix

Interpolation 2d miscibility fix
This commit is contained in:
Atgeirr Flø Rasmussen 2019-04-02 08:20:02 +02:00 committed by GitHub
commit 77c963167f
5 changed files with 88 additions and 65 deletions

View File

@ -42,6 +42,9 @@
namespace Opm {
/*!
* \brief Implements a scalar function that depends on two variables and which is sampled
* uniformly in the X direction, but non-uniformly on the Y axis-
@ -56,7 +59,23 @@ class UniformXTabulated2DFunction
typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
public:
UniformXTabulated2DFunction()
/*!
* \brief Indicates how interpolation will be performed.
*
* Normal interpolation is done by interpolating vertically
* between lines of sample points, whereas LeftExtreme or
* RightExtreme implies guided interpolation, where
* interpolation is done parallel to a guide line. With
* LeftExtreme the lowest Y values will be used for the guide,
* and the guide line slope extends unchanged to infinity. With
* RightExtreme, the highest Y values are used, and the slope
* decreases linearly down to 0 (normal interpolation) for y <= 0.
*/
enum class InterpolationPolicy { LeftExtreme, RightExtreme, Vertical };
explicit UniformXTabulated2DFunction(const InterpolationPolicy iGuide)
: interpGuide_(iGuide)
{ }
/*!
@ -281,11 +300,43 @@ public:
// table ...
unsigned i = xSegmentIndex(x, extrapolate);
const Evaluation& alpha = xToAlpha(x, i);
// The 'shift' is used to shift the points used to interpolate within
// the (i) and (i+1) sets of sample points, so that when approaching
// the boundary of the domain given by the samples, one gets the same
// value as one would get by interpolating along the boundary curve
// itself.
Evaluation shift = 0.0;
if (interpGuide_ == InterpolationPolicy::Vertical) {
// Shift is zero, no need to reset it.
} else {
// find upper and lower y value
if (interpGuide_ == InterpolationPolicy::LeftExtreme) {
// The domain is above the boundary curve, up to y = infinity.
// The shift is therefore the same for all values of y.
shift = yPos_[i+1] - yPos_[i];
} else {
assert(interpGuide_ == InterpolationPolicy::RightExtreme);
// The domain is below the boundary curve, down to y = 0.
// The shift is therefore no longer the the same for all
// values of y, since at y = 0 the shift must be zero.
// The shift is computed by linear interpolation between
// the maximal value at the domain boundary curve, and zero.
shift = yPos_[i+1] - yPos_[i];
auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
if (yEnd > 0.) {
shift = shift * y / yEnd;
} else {
shift = 0.;
}
}
}
auto yLower = y - alpha*shift;
auto yUpper = y + (1-alpha)*shift;
unsigned j1 = ySegmentIndex(y, i, extrapolate);
unsigned j2 = ySegmentIndex(y, i + 1, extrapolate);
const Evaluation& beta1 = yToBeta(y, i, j1);
const Evaluation& beta2 = yToBeta(y, i + 1, j2);
unsigned j1 = ySegmentIndex(yLower, i, extrapolate);
unsigned j2 = ySegmentIndex(yUpper, i + 1, extrapolate);
const Evaluation& beta1 = yToBeta(yLower, i, j1);
const Evaluation& beta2 = yToBeta(yUpper, i + 1, j2);
// evaluate the two function values for the same y value ...
const Evaluation& s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
@ -310,12 +361,14 @@ public:
{
if (xPos_.empty() || xPos_.back() < nextX) {
xPos_.push_back(nextX);
samples_.resize(xPos_.size());
yPos_.push_back(-1e100);
samples_.push_back({});
return xPos_.size() - 1;
}
else if (xPos_.front() > nextX) {
// this is slow, but so what?
xPos_.insert(xPos_.begin(), nextX);
yPos_.insert(yPos_.begin(), -1e100);
samples_.insert(samples_.begin(), std::vector<SamplePoint>());
return 0;
}
@ -326,20 +379,25 @@ public:
/*!
* \brief Append a sample point.
*
* Returns the i index of that line.
* Returns the i index of the new point within its line.
*/
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
{
assert(0 <= i && i < numX());
Scalar x = iToX(i);
if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
samples_[i].push_back(SamplePoint(x, y, value));
if (interpGuide_ == InterpolationPolicy::RightExtreme) {
yPos_[i] = y;
}
return samples_[i].size() - 1;
}
else if (std::get<1>(samples_[i].front()) > y) {
// slow, but we still don't care...
samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
if (interpGuide_ == InterpolationPolicy::LeftExtreme) {
yPos_[i] = y;
}
return 0;
}
@ -388,6 +446,11 @@ private:
// the position of each vertical line on the x-axis
std::vector<Scalar> xPos_;
// the position on the y-axis of the guide point
std::vector<Scalar> yPos_;
InterpolationPolicy interpGuide_;
};
} // namespace Opm

View File

@ -724,20 +724,10 @@ public:
if (fluidState.saturation(gasPhaseIdx) > 0.0
&& Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, Opm::scalarValue(T), Opm::scalarValue(p)))
{
if (fluidState.saturation(gasPhaseIdx) < 1e-4) {
// here comes the relatively expensive case: first calculate and then
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
const auto& bSat = oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
const auto& bUndersat = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
return alpha*bSat + (1.0 - alpha)*bUndersat;
}
return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
} else {
return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
}
return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
}
const LhsEval Rs(0.0);
@ -749,20 +739,10 @@ public:
if (fluidState.saturation(oilPhaseIdx) > 0.0
&& Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, Opm::scalarValue(T), Opm::scalarValue(p)))
{
if (fluidState.saturation(oilPhaseIdx) < 1e-4) {
// here comes the relatively expensive case: first calculate and then
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
const auto& bSat = gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
const auto& bUndersat = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
return alpha*bSat + (1.0 - alpha)*bUndersat;
}
return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
} else {
return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
}
return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
}
const LhsEval Rv(0.0);
@ -941,20 +921,10 @@ public:
if (fluidState.saturation(gasPhaseIdx) > 0.0
&& Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, Opm::scalarValue(T), Opm::scalarValue(p)))
{
if (fluidState.saturation(gasPhaseIdx) < 1e-4) {
// here comes the relatively expensive case: first calculate and then
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
const auto& muSat = oilPvt_->saturatedViscosity(regionIdx, T, p);
const auto& muUndersat = oilPvt_->viscosity(regionIdx, T, p, Rs);
return alpha*muSat + (1.0 - alpha)*muUndersat;
}
return oilPvt_->saturatedViscosity(regionIdx, T, p);
} else {
return oilPvt_->viscosity(regionIdx, T, p, Rs);
}
return oilPvt_->viscosity(regionIdx, T, p, Rs);
}
const LhsEval Rs(0.0);
@ -967,20 +937,10 @@ public:
if (fluidState.saturation(oilPhaseIdx) > 0.0
&& Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, Opm::scalarValue(T), Opm::scalarValue(p)))
{
if (fluidState.saturation(oilPhaseIdx) < 1e-4) {
// here comes the relatively expensive case: first calculate and then
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
const auto& muSat = gasPvt_->saturatedViscosity(regionIdx, T, p);
const auto& muUndersat = gasPvt_->viscosity(regionIdx, T, p, Rv);
return alpha*muSat + (1.0 - alpha)*muUndersat;
}
return gasPvt_->saturatedViscosity(regionIdx, T, p);
} else {
return gasPvt_->viscosity(regionIdx, T, p, Rv);
}
return gasPvt_->viscosity(regionIdx, T, p, Rv);
}
const LhsEval Rv(0.0);

View File

@ -238,11 +238,11 @@ public:
{
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseOilBTable_.resize(numRegions);
inverseOilBMuTable_.resize(numRegions);
inverseOilBTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
inverseOilBMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
inverseSaturatedOilBTable_.resize(numRegions);
inverseSaturatedOilBMuTable_.resize(numRegions);
oilMuTable_.resize(numRegions);
oilMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
saturatedOilMuTable_.resize(numRegions);
saturatedGasDissolutionFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);

View File

@ -234,11 +234,11 @@ public:
{
oilReferenceDensity_.resize(numRegions);
gasReferenceDensity_.resize(numRegions);
inverseGasB_.resize(numRegions);
inverseGasBMu_.resize(numRegions);
inverseGasB_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseGasBMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
inverseSaturatedGasB_.resize(numRegions);
inverseSaturatedGasBMu_.resize(numRegions);
gasMu_.resize(numRegions);
gasMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
saturatedOilVaporizationFactorTable_.resize(numRegions);
saturationPressure_.resize(numRegions);
}

View File

@ -90,8 +90,7 @@ struct Test
Scalar yMin = -1/2.0;
Scalar yMax = 1/3.0;
unsigned n = 40;
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>(Opm::UniformXTabulated2DFunction<Scalar>::InterpolationPolicy::Vertical);
for (unsigned i = 0; i < m; ++i) {
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
tab->appendXPos(x);
@ -116,7 +115,8 @@ struct Test
Scalar yMin = - 4.0;
Scalar yMax = 5.0;
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>(Opm::UniformXTabulated2DFunction<Scalar>::InterpolationPolicy::Vertical);
for (unsigned i = 0; i < m; ++i) {
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
tab->appendXPos(x);