27#ifndef OPM_TABULATED_1D_FUNCTION_HPP
28#define OPM_TABULATED_1D_FUNCTION_HPP
30#include <opm/common/OpmLog/OpmLog.hpp>
49template <
class Scalar>
68 template <
class ScalarArrayX,
class ScalarArrayY>
70 const ScalarArrayX& x,
71 const ScalarArrayY& y,
72 bool sortInputs =
true)
83 template <
class ScalarContainer>
85 const ScalarContainer& y,
86 bool sortInputs =
true)
95 template <
class Po
intContainer>
97 bool sortInputs =
true)
105 template <
class ScalarArrayX,
class ScalarArrayY>
107 const ScalarArrayX& x,
108 const ScalarArrayY& y,
109 bool sortInputs =
true)
111 assert(nSamples > 1);
113 resizeArrays_(nSamples);
114 for (
size_t i = 0; i < nSamples; ++i) {
121 else if (xValues_[0] > xValues_[
numSamples() - 1])
122 reverseSamplingPoints_();
130 template <
class ScalarContainerX,
class ScalarContainerY>
132 const ScalarContainerY& y,
133 bool sortInputs =
true)
135 assert(x.size() == y.size());
137 resizeArrays_(x.size());
139 std::copy(x.begin(), x.end(), xValues_.begin());
140 std::copy(y.begin(), y.end(), yValues_.begin());
144 else if (xValues_[0] > xValues_[
numSamples() - 1])
145 reverseSamplingPoints_();
152 template <
class Po
intArray>
154 const PointArray& points,
155 bool sortInputs =
true)
159 assert(nSamples > 1);
161 resizeArrays_(nSamples);
162 for (
size_t i = 0; i < nSamples; ++i) {
163 xValues_[i] = points[i][0];
164 yValues_[i] = points[i][1];
169 else if (xValues_[0] > xValues_[
numSamples() - 1])
170 reverseSamplingPoints_();
187 template <
class XYContainer>
189 bool sortInputs =
true)
193 assert(points.size() > 1);
195 resizeArrays_(points.size());
196 typename XYContainer::const_iterator it = points.begin();
197 typename XYContainer::const_iterator endIt = points.end();
198 for (
unsigned i = 0; it != endIt; ++i, ++it) {
199 xValues_[i] = std::get<0>(*it);
200 yValues_[i] = std::get<1>(*it);
205 else if (xValues_[0] > xValues_[
numSamples() - 1])
206 reverseSamplingPoints_();
213 {
return xValues_.size(); }
219 {
return xValues_[0]; }
230 Scalar
xAt(
size_t i)
const
231 {
return xValues_[i]; }
233 const std::vector<Scalar>& xValues()
const
236 const std::vector<Scalar>& yValues()
const
243 {
return yValues_[i]; }
248 template <
class Evaluation>
250 {
return xValues_[0] <= x && x <= xValues_[
numSamples() - 1]; }
261 template <
class Evaluation>
262 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const
265 return eval(x, segIdx);
268 template <
class Evaluation>
271 size_t segIdx = segIdxIn.value;
272 Scalar x0 = xValues_[segIdx];
273 Scalar x1 = xValues_[segIdx + 1];
275 Scalar y0 = yValues_[segIdx];
276 Scalar y1 = yValues_[segIdx + 1];
278 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
292 template <
class Evaluation>
295 size_t segIdx = findSegmentIndex(x, extrapolate).value;
296 return evalDerivative_(x, segIdx);
313 template <
class Evaluation>
331 template <
class Evaluation>
344 [[maybe_unused]]
bool extrapolate =
false)
const
361 size_t i = findSegmentIndex(x0, extrapolate).value;
362 if (xValues_[i + 1] >= x1) {
365 updateMonotonicity_(i, r);
371 updateMonotonicity_(i, r);
376 size_t iEnd = findSegmentIndex(x1, extrapolate).value;
377 for (; i < iEnd - 1; ++i) {
378 updateMonotonicity_(i, r);
391 return (r < 0 || r==3)?-1:0;
393 return (r > 0 || r==3)?1:0;
399 updateMonotonicity_(iEnd, r);
427 void printCSV(Scalar xi0, Scalar xi1,
unsigned k, std::ostream& os)
const;
430 return xValues_ == data.xValues_ &&
431 yValues_ == data.yValues_;
434 template <
class Evaluation>
435 SegmentIndex findSegmentIndex(
const Evaluation& x,
bool extrapolate =
false)
const
438 throw std::runtime_error(
"We can not search for extrapolation/interpolation "
439 "segment in an 1D table for non-finite value " +
440 std::to_string(getValue(x)) +
" .");
443 if (!extrapolate && !
applies(x))
444 throw std::logic_error(
"Trying to evaluate a tabulated function outside of its range");
448 throw std::logic_error(
"We need at least two sampling points to "
449 "do interpolation/extrapolation, "
450 "and the table only contains " +
455 if (x <= xValues_[1])
456 return SegmentIndex{0};
457 else if (x >= xValues_[xValues_.size() - 2])
458 return SegmentIndex{xValues_.size() - 2};
462 size_t upperIdx = xValues_.size() - 2;
463 while (lowerIdx + 1 < upperIdx) {
464 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
465 if (x < xValues_[pivotIdx])
471 if (xValues_[lowerIdx] > x || x > xValues_[lowerIdx + 1]) {
472 std::string msg =
"Problematic interpolation/extrapolation "
473 "segment is found for the input value " +
474 std::to_string(Opm::getValue(x)) +
475 "\nthe lower index of the found segment is " +
476 std::to_string(lowerIdx) +
477 ", the size of the table is " +
479 ",\nand the end values of the found segment are " +
480 std::to_string(xValues_[lowerIdx]) +
482 std::to_string(xValues_[lowerIdx + 1]) +
484 msg +=
"Outputting the problematic table for more information "
485 "(with *** marking the found segment):";
491 msg +=
" " + std::to_string(xValues_[i]);
492 if (i == lowerIdx + 1)
497 throw std::runtime_error(msg);
499 return SegmentIndex{lowerIdx};
504 template <
class Evaluation>
505 Evaluation evalDerivative_(
const Evaluation& x,
size_t segIdx)
const
508 Scalar x0 = xValues_[segIdx];
509 Scalar x1 = xValues_[segIdx + 1];
511 Scalar y0 = yValues_[segIdx];
512 Scalar y1 = yValues_[segIdx + 1];
514 Evaluation ret = blank(x);
515 ret = (y1 - y0)/(x1 - x0);
527 int updateMonotonicity_(
size_t i,
int& r)
const
529 if (yValues_[i] < yValues_[i + 1]) {
531 if (r == 3 || r == 1)
537 else if (yValues_[i] > yValues_[i + 1]) {
539 if (r == 3 || r == -1)
554 ComparatorX_(
const std::vector<Scalar>& x)
558 bool operator ()(
size_t idxA,
size_t idxB)
const
559 {
return x_.at(idxA) < x_.at(idxB); }
561 const std::vector<Scalar>& x_;
572 std::vector<unsigned> idxVector(n);
573 for (
unsigned i = 0; i < n; ++i)
578 ComparatorX_ cmp(xValues_);
579 std::sort(idxVector.begin(), idxVector.end(), cmp);
582 std::vector<Scalar> tmpX(n), tmpY(n);
583 for (
size_t i = 0; i < idxVector.size(); ++ i) {
584 tmpX[i] = xValues_[idxVector[i]];
585 tmpY[i] = yValues_[idxVector[i]];
595 void reverseSamplingPoints_()
599 for (
size_t i = 0; i <= (n - 1)/2; ++i) {
600 std::swap(xValues_[i], xValues_[n - i - 1]);
601 std::swap(yValues_[i], yValues_[n - i - 1]);
608 void resizeArrays_(
size_t nSamples)
610 xValues_.resize(nSamples);
611 yValues_.resize(nSamples);
614 std::vector<Scalar> xValues_;
615 std::vector<Scalar> yValues_;
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Implements a linearly interpolated scalar function that depends on one variable.
Definition Tabulated1DFunction.hpp:51
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition Tabulated1DFunction.hpp:262
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
size_t numSamples() const
Returns the number of sampling points.
Definition Tabulated1DFunction.hpp:212
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition Tabulated1DFunction.hpp:293
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition Tabulated1DFunction.hpp:153
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition Tabulated1DFunction.hpp:224
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:84
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition Tabulated1DFunction.hpp:218
Evaluation evalThirdDerivative(const Evaluation &, bool=false) const
Evaluate the function's third derivative at a given position.
Definition Tabulated1DFunction.hpp:332
Evaluation evalSecondDerivative(const Evaluation &, bool=false) const
Evaluate the function's second derivative at a given position.
Definition Tabulated1DFunction.hpp:314
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:69
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:188
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition Tabulated1DFunction.hpp:230
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:96
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition Tabulated1DFunction.hpp:58
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval.
Definition Tabulated1DFunction.hpp:408
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:106
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition Tabulated1DFunction.hpp:131
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:343
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition Tabulated1DFunction.hpp:249
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition Tabulated1DFunction.hpp:242
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
Definition Tabulated1DFunction.hpp:41