Loading...
Searching...
No Matches
UniformXTabulated2DFunction.hpp
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
#define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
#include <opm/common/Exceptions.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <cassert>
#include <cmath>
#include <iosfwd>
#include <limits>
#include <tuple>
#include <type_traits>
#include <vector>
namespace Opm {
template <class Scalar>
class UniformXTabulated2DFunction
{
public:
typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
enum InterpolationPolicy {
LeftExtreme,
RightExtreme,
Vertical
};
explicit UniformXTabulated2DFunction(const InterpolationPolicy interpolationGuide = Vertical)
: interpolationGuide_(interpolationGuide)
{ }
UniformXTabulated2DFunction(const std::vector<Scalar>& xPos,
const std::vector<Scalar>& yPos,
const std::vector<std::vector<SamplePoint>>& samples,
InterpolationPolicy interpolationGuide)
: samples_(samples)
, xPos_(xPos)
, yPos_(yPos)
, interpolationGuide_(interpolationGuide)
{ }
Scalar xMin() const
{ return xPos_.front(); }
Scalar xMax() const
{ return xPos_.back(); }
Scalar xAt(size_t i) const
{ return xPos_[i]; }
Scalar yAt(size_t i, size_t j) const
{ return std::get<1>(samples_[i][j]); }
Scalar valueAt(size_t i, size_t j) const
{ return std::get<2>(samples_[i][j]); }
size_t numX() const
{ return xPos_.size(); }
Scalar yMin(unsigned i) const
{ return std::get<1>(samples_.at(i).front()); }
Scalar yMax(unsigned i) const
{ return std::get<1>(samples_.at(i).back()); }
size_t numY(unsigned i) const
{ return samples_.at(i).size(); }
Scalar iToX(unsigned i) const
{
assert(i < numX());
return xPos_.at(i);
}
const std::vector<std::vector<SamplePoint>>& samples() const
{
return samples_;
}
const std::vector<Scalar>& xPos() const
{
return xPos_;
}
const std::vector<Scalar>& yPos() const
{
return yPos_;
}
InterpolationPolicy interpolationGuide() const
{
return interpolationGuide_;
}
Scalar jToY(unsigned i, unsigned j) const
{
assert(i < numX());
assert(size_t(j) < samples_[i].size());
return std::get<1>(samples_.at(i).at(j));
}
template <class Evaluation>
unsigned xSegmentIndex(const Evaluation& x,
[[maybe_unused]] bool extrapolate = false) const
{
assert(extrapolate || (xMin() <= x && x <= xMax()));
// we need at least two sampling points!
assert(xPos_.size() >= 2);
if (x <= xPos_[1])
return 0;
else if (x >= xPos_[xPos_.size() - 2])
return xPos_.size() - 2;
else {
assert(xPos_.size() >= 3);
// bisection
unsigned lowerIdx = 1;
unsigned upperIdx = xPos_.size() - 2;
while (lowerIdx + 1 < upperIdx) {
unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
if (x < xPos_[pivotIdx])
upperIdx = pivotIdx;
else
lowerIdx = pivotIdx;
}
return lowerIdx;
}
}
template <class Evaluation>
Evaluation xToAlpha(const Evaluation& x, unsigned segmentIdx) const
{
Scalar x1 = xPos_[segmentIdx];
Scalar x2 = xPos_[segmentIdx + 1];
return (x - x1)/(x2 - x1);
}
template <class Evaluation>
unsigned ySegmentIndex(const Evaluation& y, unsigned xSampleIdx,
[[maybe_unused]] bool extrapolate = false) const
{
assert(xSampleIdx < numX());
const auto& colSamplePoints = samples_.at(xSampleIdx);
assert(colSamplePoints.size() >= 2);
assert(extrapolate || (yMin(xSampleIdx) <= y && y <= yMax(xSampleIdx)));
if (y <= std::get<1>(colSamplePoints[1]))
return 0;
else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
return colSamplePoints.size() - 2;
else {
assert(colSamplePoints.size() >= 3);
// bisection
unsigned lowerIdx = 1;
unsigned upperIdx = colSamplePoints.size() - 2;
while (lowerIdx + 1 < upperIdx) {
unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
if (y < std::get<1>(colSamplePoints[pivotIdx]))
upperIdx = pivotIdx;
else
lowerIdx = pivotIdx;
}
return lowerIdx;
}
}
template <class Evaluation>
Evaluation yToBeta(const Evaluation& y, unsigned xSampleIdx, unsigned ySegmentIdx) const
{
assert(xSampleIdx < numX());
assert(ySegmentIdx < numY(xSampleIdx) - 1);
const auto& colSamplePoints = samples_.at(xSampleIdx);
Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
return (y - y1)/(y2 - y1);
}
template <class Evaluation>
bool applies(const Evaluation& x, const Evaluation& y) const
{
if (x < xMin() || xMax() < x)
return false;
unsigned i = xSegmentIndex(x, /*extrapolate=*/false);
Scalar alpha = xToAlpha(decay<Scalar>(x), i);
const auto& col1SamplePoints = samples_.at(i);
const auto& col2SamplePoints = samples_.at(i + 1);
Scalar minY =
alpha*std::get<1>(col1SamplePoints.front()) +
(1 - alpha)*std::get<1>(col2SamplePoints.front());
Scalar maxY =
alpha*std::get<1>(col1SamplePoints.back()) +
(1 - alpha)*std::get<1>(col2SamplePoints.back());
return minY <= y && y <= maxY;
}
template <class Evaluation>
Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
{
Evaluation alpha, beta1, beta2;
unsigned i, j1, j2;
findPoints(i, j1, j2, alpha, beta1, beta2, x, y, extrapolate);
return eval(i, j1, j2, alpha, beta1, beta2);
}
template <class Evaluation>
void findPoints(unsigned& i,
unsigned& j1,
unsigned& j2,
Evaluation& alpha,
Evaluation& beta1,
Evaluation& beta2,
const Evaluation& x,
const Evaluation& y,
bool extrapolate) const
{
#ifndef NDEBUG
if (!extrapolate && !applies(x, y)) {
if constexpr (std::is_floating_point_v<Evaluation>) {
throw NumericalProblem("Attempt to get undefined table value (" +
std::to_string(x) + ", " +
std::to_string(y) + ")");
} else {
throw NumericalProblem("Attempt to get undefined table value (" +
std::to_string(x.value()) + ", " +
std::to_string(y.value()) + ")");
}
};
#endif
// bi-linear interpolation: first, calculate the x and y indices in the lookup
// table ...
i = xSegmentIndex(x, extrapolate);
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 (interpolationGuide_ == InterpolationPolicy::Vertical) {
// Shift is zero, no need to reset it.
} else {
// find upper and lower y value
if (interpolationGuide_ == 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(interpolationGuide_ == 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;
j1 = ySegmentIndex(yLower, i, extrapolate);
j2 = ySegmentIndex(yUpper, i + 1, extrapolate);
beta1 = yToBeta(yLower, i, j1);
beta2 = yToBeta(yUpper, i + 1, j2);
}
template <class Evaluation>
Evaluation eval(const unsigned& i, const unsigned& j1, const unsigned& j2, const Evaluation& alpha,const Evaluation& beta1,const Evaluation& beta2) const
{
// evaluate the two function values for the same y value ...
const Evaluation& s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
const Evaluation& s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
Valgrind::CheckDefined(s1);
Valgrind::CheckDefined(s2);
// ... and combine them using the x position
const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
Valgrind::CheckDefined(result);
return result;
}
size_t appendXPos(Scalar nextX)
{
if (xPos_.empty() || xPos_.back() < nextX) {
xPos_.push_back(nextX);
yPos_.push_back(std::numeric_limits<Scalar>::lowest() / 2);
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(), std::numeric_limits<Scalar>::lowest() / 2);
samples_.insert(samples_.begin(), std::vector<SamplePoint>());
return 0;
}
throw std::invalid_argument("Sampling points should be specified either monotonically "
"ascending or descending.");
}
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
{
assert(i < numX());
Scalar x = iToX(i);
if (samples_[i].empty()) {
samples_[i].emplace_back(x, y, value);
yPos_[i] = y;
return 0;
}
else if (std::get<1>(samples_[i].back()) < y) {
samples_[i].emplace_back(x, y, value);
if (interpolationGuide_ == 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].emplace(samples_[i].begin(), x, y, value);
if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
yPos_[i] = y;
}
return 0;
}
throw std::invalid_argument("Sampling points must be specified in either monotonically "
"ascending or descending order.");
}
void print(std::ostream& os) const;
bool operator==(const UniformXTabulated2DFunction<Scalar>& data) const {
return this->xPos() == data.xPos() &&
this->yPos() == data.yPos() &&
this->samples() == data.samples() &&
this->interpolationGuide() == data.interpolationGuide();
}
private:
// the vector which contains the values of the sample points
// f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
// instead!
std::vector<std::vector<SamplePoint> > samples_;
// 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 interpolationGuide_;
};
} // namespace Opm
#endif
Exceptions.hpp
Provides the OPM specific exception classes.
MathToolbox.hpp
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Valgrind.hpp
Some templates to wrap the valgrind client request macros.
Opm::NumericalProblem
Definition Exceptions.hpp:40
Opm::UniformXTabulated2DFunction
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition UniformXTabulated2DFunction.hpp:55
Opm::UniformXTabulated2DFunction::xMax
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:98
Opm::UniformXTabulated2DFunction::appendSamplePoint
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
Append a sample point.
Definition UniformXTabulated2DFunction.hpp:441
Opm::UniformXTabulated2DFunction::xToAlpha
Evaluation xToAlpha(const Evaluation &x, unsigned segmentIdx) const
Return the relative position of an x value in an intervall.
Definition UniformXTabulated2DFunction.hpp:225
Opm::UniformXTabulated2DFunction::yToBeta
Evaluation yToBeta(const Evaluation &y, unsigned xSampleIdx, unsigned ySegmentIdx) const
Return the relative position of an y value in an interval.
Definition UniformXTabulated2DFunction.hpp:274
Opm::UniformXTabulated2DFunction::xAt
Scalar xAt(size_t i) const
Returns the value of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:104
Opm::UniformXTabulated2DFunction::jToY
Scalar jToY(unsigned i, unsigned j) const
Return the position on the y-axis of the j-th interval.
Definition UniformXTabulated2DFunction.hpp:176
Opm::UniformXTabulated2DFunction::InterpolationPolicy
InterpolationPolicy
Indicates how interpolation will be performed.
Definition UniformXTabulated2DFunction.hpp:69
Opm::UniformXTabulated2DFunction::valueAt
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition UniformXTabulated2DFunction.hpp:116
Opm::UniformXTabulated2DFunction::ySegmentIndex
unsigned ySegmentIndex(const Evaluation &y, unsigned xSampleIdx, bool extrapolate=false) const
Return the interval index of a given position on the y-axis.
Definition UniformXTabulated2DFunction.hpp:236
Opm::UniformXTabulated2DFunction::xMin
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:92
Opm::UniformXTabulated2DFunction::print
void print(std::ostream &os) const
Print the table for debugging purposes.
Definition UniformXTabulated2DFunction.cpp:32
Opm::UniformXTabulated2DFunction::yMin
Scalar yMin(unsigned i) const
Returns the minimum of the Y coordinate of the sampling points for a given column.
Definition UniformXTabulated2DFunction.hpp:128
Opm::UniformXTabulated2DFunction::xSegmentIndex
unsigned xSegmentIndex(const Evaluation &x, bool extrapolate=false) const
Return the interval index of a given position on the x-axis.
Definition UniformXTabulated2DFunction.hpp:188
Opm::UniformXTabulated2DFunction::numY
size_t numY(unsigned i) const
Returns the number of sampling points in Y direction a given column.
Definition UniformXTabulated2DFunction.hpp:140
Opm::UniformXTabulated2DFunction::numX
size_t numX() const
Returns the number of sampling points in X direction.
Definition UniformXTabulated2DFunction.hpp:122
Opm::UniformXTabulated2DFunction::yMax
Scalar yMax(unsigned i) const
Returns the maximum of the Y coordinate of the sampling points for a given column.
Definition UniformXTabulated2DFunction.hpp:134
Opm::UniformXTabulated2DFunction::iToX
Scalar iToX(unsigned i) const
Return the position on the x-axis of the i-th interval.
Definition UniformXTabulated2DFunction.hpp:146
Opm::UniformXTabulated2DFunction::yAt
Scalar yAt(size_t i, size_t j) const
Returns the value of the Y coordinate of a sampling point.
Definition UniformXTabulated2DFunction.hpp:110
Opm::UniformXTabulated2DFunction::appendXPos
size_t appendXPos(Scalar nextX)
Set the x-position of a vertical line.
Definition UniformXTabulated2DFunction.hpp:417
Opm::UniformXTabulated2DFunction::eval
Evaluation eval(const Evaluation &x, const Evaluation &y, bool extrapolate=false) const
Evaluate the function at a given (x,y) position.
Definition UniformXTabulated2DFunction.hpp:319
Opm::UniformXTabulated2DFunction::applies
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true iff a coordinate lies in the tabulated range.
Definition UniformXTabulated2DFunction.hpp:291
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30