Loading...
Searching...
No Matches
Tabulated1DFunction.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_TABULATED_1D_FUNCTION_HPP
#define OPM_TABULATED_1D_FUNCTION_HPP
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/material/densead/Math.hpp>
#include <algorithm>
#include <cassert>
#include <iosfwd>
#include <stdexcept>
#include <vector>
namespace Opm {
struct SegmentIndex {
size_t value;
};
template <class Scalar>
class Tabulated1DFunction
{
public:
Tabulated1DFunction()
{}
template <class ScalarArrayX, class ScalarArrayY>
Tabulated1DFunction(size_t nSamples,
const ScalarArrayX& x,
const ScalarArrayY& y,
bool sortInputs = true)
{ this->setXYArrays(nSamples, x, y, sortInputs); }
template <class ScalarContainer>
Tabulated1DFunction(const ScalarContainer& x,
const ScalarContainer& y,
bool sortInputs = true)
{ this->setXYContainers(x, y, sortInputs); }
template <class PointContainer>
explicit Tabulated1DFunction(const PointContainer& points,
bool sortInputs = true)
{ this->setContainerOfTuples(points, sortInputs); }
template <class ScalarArrayX, class ScalarArrayY>
void setXYArrays(size_t nSamples,
const ScalarArrayX& x,
const ScalarArrayY& y,
bool sortInputs = true)
{
assert(nSamples > 1);
resizeArrays_(nSamples);
for (size_t i = 0; i < nSamples; ++i) {
xValues_[i] = x[i];
yValues_[i] = y[i];
}
if (sortInputs)
sortInput_();
else if (xValues_[0] > xValues_[numSamples() - 1])
reverseSamplingPoints_();
}
template <class ScalarContainerX, class ScalarContainerY>
void setXYContainers(const ScalarContainerX& x,
const ScalarContainerY& y,
bool sortInputs = true)
{
assert(x.size() == y.size());
resizeArrays_(x.size());
if (x.size() > 0) {
std::copy(x.begin(), x.end(), xValues_.begin());
std::copy(y.begin(), y.end(), yValues_.begin());
if (sortInputs)
sortInput_();
else if (xValues_[0] > xValues_[numSamples() - 1])
reverseSamplingPoints_();
}
}
template <class PointArray>
void setArrayOfPoints(size_t nSamples,
const PointArray& points,
bool sortInputs = true)
{
// a linear function with less than two sampling points? what an incredible
// bad idea!
assert(nSamples > 1);
resizeArrays_(nSamples);
for (size_t i = 0; i < nSamples; ++i) {
xValues_[i] = points[i][0];
yValues_[i] = points[i][1];
}
if (sortInputs)
sortInput_();
else if (xValues_[0] > xValues_[numSamples() - 1])
reverseSamplingPoints_();
}
template <class XYContainer>
void setContainerOfTuples(const XYContainer& points,
bool sortInputs = true)
{
// a linear function with less than two sampling points? what an incredible
// bad idea!
assert(points.size() > 1);
resizeArrays_(points.size());
typename XYContainer::const_iterator it = points.begin();
typename XYContainer::const_iterator endIt = points.end();
for (unsigned i = 0; it != endIt; ++i, ++it) {
xValues_[i] = std::get<0>(*it);
yValues_[i] = std::get<1>(*it);
}
if (sortInputs)
sortInput_();
else if (xValues_[0] > xValues_[numSamples() - 1])
reverseSamplingPoints_();
}
size_t numSamples() const
{ return xValues_.size(); }
Scalar xMin() const
{ return xValues_[0]; }
Scalar xMax() const
{ return xValues_[numSamples() - 1]; }
Scalar xAt(size_t i) const
{ return xValues_[i]; }
const std::vector<Scalar>& xValues() const
{ return xValues_; }
const std::vector<Scalar>& yValues() const
{ return yValues_; }
Scalar valueAt(size_t i) const
{ return yValues_[i]; }
template <class Evaluation>
bool applies(const Evaluation& x) const
{ return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
template <class Evaluation>
Evaluation eval(const Evaluation& x, bool extrapolate = false) const
{
SegmentIndex segIdx = findSegmentIndex(x, extrapolate);
return eval(x, segIdx);
}
template <class Evaluation>
Evaluation eval(const Evaluation& x, SegmentIndex segIdxIn) const
{
size_t segIdx = segIdxIn.value;
Scalar x0 = xValues_[segIdx];
Scalar x1 = xValues_[segIdx + 1];
Scalar y0 = yValues_[segIdx];
Scalar y1 = yValues_[segIdx + 1];
return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
}
template <class Evaluation>
Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
{
size_t segIdx = findSegmentIndex(x, extrapolate).value;
return evalDerivative_(x, segIdx);
}
template <class Evaluation>
Evaluation evalSecondDerivative(const Evaluation&, bool = false) const
{ return 0.0; }
template <class Evaluation>
Evaluation evalThirdDerivative(const Evaluation&, bool = false) const
{ return 0.0; }
int monotonic(Scalar x0, Scalar x1,
[[maybe_unused]] bool extrapolate = false) const
{
assert(x0 != x1);
// make sure that x0 is smaller than x1
if (x0 > x1)
std::swap(x0, x1);
assert(x0 < x1);
int r = 3;
if (x0 < xMin()) {
assert(extrapolate);
x0 = xMin();
};
size_t i = findSegmentIndex(x0, extrapolate).value;
if (xValues_[i + 1] >= x1) {
// interval is fully contained within a single function
// segment
updateMonotonicity_(i, r);
return r;
}
// the first segment overlaps with the specified interval
// partially
updateMonotonicity_(i, r);
++ i;
// make sure that the segments which are completly in the
// interval [x0, x1] all exhibit the same monotonicity.
size_t iEnd = findSegmentIndex(x1, extrapolate).value;
for (; i < iEnd - 1; ++i) {
updateMonotonicity_(i, r);
if (!r)
return 0;
}
// if the user asked for a part of the function which is
// extrapolated, we need to check the slope at the function's
// endpoint
if (x1 > xMax()) {
assert(extrapolate);
Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2);
if (m < 0)
return (r < 0 || r==3) ? -1 : 0;
else if (m > 0)
return r > 0 ? 1 : 0;
return r;
}
// check for the last segment
updateMonotonicity_(iEnd, r);
return r;
}
int monotonic() const
{ return monotonic(xMin(), xMax()); }
void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const;
bool operator==(const Tabulated1DFunction<Scalar>& data) const {
return xValues_ == data.xValues_ &&
yValues_ == data.yValues_;
}
template <class Evaluation>
SegmentIndex findSegmentIndex(const Evaluation& x, bool extrapolate = false) const
{
if (!isfinite(x)) {
throw std::runtime_error("We can not search for extrapolation/interpolation "
"segment in an 1D table for non-finite value " +
std::to_string(getValue(x)) + " .");
}
if (!extrapolate && !applies(x))
throw std::logic_error("Trying to evaluate a tabulated function outside of its range");
// we need at least two sampling points!
if (numSamples() < 2) {
throw std::logic_error("We need at least two sampling points to "
"do interpolation/extrapolation, "
"and the table only contains " +
std::to_string(numSamples()) +
" sampling points");
}
if (x <= xValues_[1])
return SegmentIndex{0};
else if (x >= xValues_[xValues_.size() - 2])
return SegmentIndex{xValues_.size() - 2};
else {
// bisection
size_t lowerIdx = 1;
size_t upperIdx = xValues_.size() - 2;
while (lowerIdx + 1 < upperIdx) {
size_t pivotIdx = (lowerIdx + upperIdx) / 2;
if (x < xValues_[pivotIdx])
upperIdx = pivotIdx;
else
lowerIdx = pivotIdx;
}
if (xValues_[lowerIdx] > x || x > xValues_[lowerIdx + 1]) {
std::string msg = "Problematic interpolation/extrapolation "
"segment is found for the input value " +
std::to_string(Opm::getValue(x)) +
"\nthe lower index of the found segment is " +
std::to_string(lowerIdx) +
", the size of the table is " +
std::to_string(numSamples()) +
",\nand the end values of the found segment are " +
std::to_string(xValues_[lowerIdx]) +
" and " +
std::to_string(xValues_[lowerIdx + 1]) +
", respectively.\n";
msg += "Outputting the problematic table for more information "
"(with *** marking the found segment):";
for (size_t i = 0; i < numSamples(); ++i) {
if (i % 10 == 0)
msg += "\n";
if (i == lowerIdx)
msg += " ***";
msg += " " + std::to_string(xValues_[i]);
if (i == lowerIdx + 1)
msg += " ***";
}
msg += "\n";
OpmLog::debug(msg);
throw std::runtime_error(msg);
}
return SegmentIndex{lowerIdx};
}
}
private:
template <class Evaluation>
Evaluation evalDerivative_(const Evaluation& x, size_t segIdx) const
{
Scalar x0 = xValues_[segIdx];
Scalar x1 = xValues_[segIdx + 1];
Scalar y0 = yValues_[segIdx];
Scalar y1 = yValues_[segIdx + 1];
Evaluation ret = blank(x);
ret = (y1 - y0)/(x1 - x0);
return ret;
}
// returns the monotonicity of a segment
//
// The return value have the following meaning:
//
// 3: function is constant within interval [x0, x1]
// 1: function is monotonously increasing in the specified interval
// 0: function is not monotonic in the specified interval
// -1: function is monotonously decreasing in the specified interval
int updateMonotonicity_(size_t i, int& r) const
{
if (yValues_[i] < yValues_[i + 1]) {
// monotonically increasing?
if (r == 3 || r == 1)
r = 1;
else
r = 0;
return 1;
}
else if (yValues_[i] > yValues_[i + 1]) {
// monotonically decreasing?
if (r == 3 || r == -1)
r = -1;
else
r = 0;
return -1;
}
return 3;
}
struct ComparatorX_
{
explicit ComparatorX_(const std::vector<Scalar>& x)
: x_(x)
{}
bool operator ()(size_t idxA, size_t idxB) const
{ return x_.at(idxA) < x_.at(idxB); }
const std::vector<Scalar>& x_;
};
void sortInput_()
{
size_t n = numSamples();
// create a vector containing 0...n-1
std::vector<unsigned> idxVector(n);
for (unsigned i = 0; i < n; ++i)
idxVector[i] = i;
// sort the indices according to the x values of the sample
// points
ComparatorX_ cmp(xValues_);
std::sort(idxVector.begin(), idxVector.end(), cmp);
// reorder the sample points
std::vector<Scalar> tmpX(n), tmpY(n);
for (size_t i = 0; i < idxVector.size(); ++ i) {
tmpX[i] = xValues_[idxVector[i]];
tmpY[i] = yValues_[idxVector[i]];
}
xValues_ = tmpX;
yValues_ = tmpY;
}
void reverseSamplingPoints_()
{
// reverse the arrays
size_t n = numSamples();
for (size_t i = 0; i <= (n - 1)/2; ++i) {
std::swap(xValues_[i], xValues_[n - i - 1]);
std::swap(yValues_[i], yValues_[n - i - 1]);
}
}
void resizeArrays_(size_t nSamples)
{
xValues_.resize(nSamples);
yValues_.resize(nSamples);
}
std::vector<Scalar> xValues_;
std::vector<Scalar> yValues_;
};
} // namespace Opm
#endif
Math.hpp
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Opm::Tabulated1DFunction::eval
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition Tabulated1DFunction.hpp:265
Opm::Tabulated1DFunction::printCSV
void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream &os) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition Tabulated1DFunction.cpp:33
Opm::Tabulated1DFunction::numSamples
size_t numSamples() const
Returns the number of sampling points.
Definition Tabulated1DFunction.hpp:215
Opm::Tabulated1DFunction::evalDerivative
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition Tabulated1DFunction.hpp:296
Opm::Tabulated1DFunction::setArrayOfPoints
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition Tabulated1DFunction.hpp:156
Opm::Tabulated1DFunction::xMax
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition Tabulated1DFunction.hpp:227
Opm::Tabulated1DFunction::Tabulated1DFunction
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:86
Opm::Tabulated1DFunction::xMin
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition Tabulated1DFunction.hpp:221
Opm::Tabulated1DFunction::evalThirdDerivative
Evaluation evalThirdDerivative(const Evaluation &, bool=false) const
Evaluate the function's third derivative at a given position.
Definition Tabulated1DFunction.hpp:335
Opm::Tabulated1DFunction::evalSecondDerivative
Evaluation evalSecondDerivative(const Evaluation &, bool=false) const
Evaluate the function's second derivative at a given position.
Definition Tabulated1DFunction.hpp:317
Opm::Tabulated1DFunction::Tabulated1DFunction
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:70
Opm::Tabulated1DFunction::setContainerOfTuples
void setContainerOfTuples(const XYContainer &points, bool sortInputs=true)
Set the sampling points of the piecewise linear function using a STL-compatible container of tuple-li...
Definition Tabulated1DFunction.hpp:191
Opm::Tabulated1DFunction::xAt
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition Tabulated1DFunction.hpp:233
Opm::Tabulated1DFunction::Tabulated1DFunction
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:99
Opm::Tabulated1DFunction::Tabulated1DFunction
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:58
Opm::Tabulated1DFunction::monotonic
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval.
Definition Tabulated1DFunction.hpp:411
Opm::Tabulated1DFunction::setXYArrays
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition Tabulated1DFunction.hpp:109
Opm::Tabulated1DFunction::setXYContainers
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition Tabulated1DFunction.hpp:134
Opm::Tabulated1DFunction::monotonic
int monotonic(Scalar x0, Scalar x1, bool extrapolate=false) const
Returns 1 if the function is monotonically increasing, -1 if the function is mononously decreasing an...
Definition Tabulated1DFunction.hpp:346
Opm::Tabulated1DFunction::applies
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition Tabulated1DFunction.hpp:252
Opm::Tabulated1DFunction::valueAt
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition Tabulated1DFunction.hpp:245
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
Opm::SegmentIndex
Definition Tabulated1DFunction.hpp:41