My Project
Loading...
Searching...
No Matches
Tabulated1DFunction.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef OPM_TABULATED_1D_FUNCTION_HPP
28#define OPM_TABULATED_1D_FUNCTION_HPP
29
30#include <opm/common/OpmLog/OpmLog.hpp>
32
33#include <algorithm>
34#include <cassert>
35#include <iosfwd>
36#include <stdexcept>
37#include <vector>
38
39namespace Opm {
40
42 size_t value;
43};
44
49template <class Scalar>
51{
52public:
60
68 template <class ScalarArrayX, class ScalarArrayY>
69 Tabulated1DFunction(size_t nSamples,
70 const ScalarArrayX& x,
71 const ScalarArrayY& y,
72 bool sortInputs = true)
73 { this->setXYArrays(nSamples, x, y, sortInputs); }
74
83 template <class ScalarContainer>
84 Tabulated1DFunction(const ScalarContainer& x,
85 const ScalarContainer& y,
86 bool sortInputs = true)
87 { this->setXYContainers(x, y, sortInputs); }
88
95 template <class PointContainer>
96 Tabulated1DFunction(const PointContainer& points,
97 bool sortInputs = true)
98 { this->setContainerOfTuples(points, sortInputs); }
99
105 template <class ScalarArrayX, class ScalarArrayY>
106 void setXYArrays(size_t nSamples,
107 const ScalarArrayX& x,
108 const ScalarArrayY& y,
109 bool sortInputs = true)
110 {
111 assert(nSamples > 1);
112
113 resizeArrays_(nSamples);
114 for (size_t i = 0; i < nSamples; ++i) {
115 xValues_[i] = x[i];
116 yValues_[i] = y[i];
117 }
118
119 if (sortInputs)
120 sortInput_();
121 else if (xValues_[0] > xValues_[numSamples() - 1])
122 reverseSamplingPoints_();
123 }
124
130 template <class ScalarContainerX, class ScalarContainerY>
131 void setXYContainers(const ScalarContainerX& x,
132 const ScalarContainerY& y,
133 bool sortInputs = true)
134 {
135 assert(x.size() == y.size());
136
137 resizeArrays_(x.size());
138 if (x.size() > 0) {
139 std::copy(x.begin(), x.end(), xValues_.begin());
140 std::copy(y.begin(), y.end(), yValues_.begin());
141
142 if (sortInputs)
143 sortInput_();
144 else if (xValues_[0] > xValues_[numSamples() - 1])
145 reverseSamplingPoints_();
146 }
147 }
148
152 template <class PointArray>
153 void setArrayOfPoints(size_t nSamples,
154 const PointArray& points,
155 bool sortInputs = true)
156 {
157 // a linear function with less than two sampling points? what an incredible
158 // bad idea!
159 assert(nSamples > 1);
160
161 resizeArrays_(nSamples);
162 for (size_t i = 0; i < nSamples; ++i) {
163 xValues_[i] = points[i][0];
164 yValues_[i] = points[i][1];
165 }
166
167 if (sortInputs)
168 sortInput_();
169 else if (xValues_[0] > xValues_[numSamples() - 1])
170 reverseSamplingPoints_();
171 }
172
187 template <class XYContainer>
188 void setContainerOfTuples(const XYContainer& points,
189 bool sortInputs = true)
190 {
191 // a linear function with less than two sampling points? what an incredible
192 // bad idea!
193 assert(points.size() > 1);
194
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);
201 }
202
203 if (sortInputs)
204 sortInput_();
205 else if (xValues_[0] > xValues_[numSamples() - 1])
206 reverseSamplingPoints_();
207 }
208
212 size_t numSamples() const
213 { return xValues_.size(); }
214
218 Scalar xMin() const
219 { return xValues_[0]; }
220
224 Scalar xMax() const
225 { return xValues_[numSamples() - 1]; }
226
230 Scalar xAt(size_t i) const
231 { return xValues_[i]; }
232
233 const std::vector<Scalar>& xValues() const
234 { return xValues_; }
235
236 const std::vector<Scalar>& yValues() const
237 { return yValues_; }
238
242 Scalar valueAt(size_t i) const
243 { return yValues_[i]; }
244
248 template <class Evaluation>
249 bool applies(const Evaluation& x) const
250 { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
251
261 template <class Evaluation>
262 Evaluation eval(const Evaluation& x, bool extrapolate = false) const
263 {
264 SegmentIndex segIdx = findSegmentIndex(x, extrapolate);
265 return eval(x, segIdx);
266 }
267
268 template <class Evaluation>
269 Evaluation eval(const Evaluation& x, SegmentIndex segIdxIn) const
270 {
271 size_t segIdx = segIdxIn.value;
272 Scalar x0 = xValues_[segIdx];
273 Scalar x1 = xValues_[segIdx + 1];
274
275 Scalar y0 = yValues_[segIdx];
276 Scalar y1 = yValues_[segIdx + 1];
277
278 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
279 }
280
292 template <class Evaluation>
293 Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
294 {
295 size_t segIdx = findSegmentIndex(x, extrapolate).value;
296 return evalDerivative_(x, segIdx);
297 }
298
313 template <class Evaluation>
314 Evaluation evalSecondDerivative(const Evaluation&, bool = false) const
315 { return 0.0; }
316
331 template <class Evaluation>
332 Evaluation evalThirdDerivative(const Evaluation&, bool = false) const
333 { return 0.0; }
334
343 int monotonic(Scalar x0, Scalar x1,
344 [[maybe_unused]] bool extrapolate = false) const
345 {
346 assert(x0 != x1);
347
348 // make sure that x0 is smaller than x1
349 if (x0 > x1)
350 std::swap(x0, x1);
351
352 assert(x0 < x1);
353
354 int r = 3;
355 if (x0 < xMin()) {
356 assert(extrapolate);
357
358 x0 = xMin();
359 };
360
361 size_t i = findSegmentIndex(x0, extrapolate).value;
362 if (xValues_[i + 1] >= x1) {
363 // interval is fully contained within a single function
364 // segment
365 updateMonotonicity_(i, r);
366 return r;
367 }
368
369 // the first segment overlaps with the specified interval
370 // partially
371 updateMonotonicity_(i, r);
372 ++ i;
373
374 // make sure that the segments which are completly in the
375 // interval [x0, x1] all exhibit the same monotonicity.
376 size_t iEnd = findSegmentIndex(x1, extrapolate).value;
377 for (; i < iEnd - 1; ++i) {
378 updateMonotonicity_(i, r);
379 if (!r)
380 return 0;
381 }
382
383 // if the user asked for a part of the function which is
384 // extrapolated, we need to check the slope at the function's
385 // endpoint
386 if (x1 > xMax()) {
387 assert(extrapolate);
388
389 Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2);
390 if (m < 0)
391 return (r < 0 || r==3)?-1:0;
392 else if (m > 0)
393 return (r > 0 || r==3)?1:0;
394
395 return r;
396 }
397
398 // check for the last segment
399 updateMonotonicity_(iEnd, r);
400
401 return r;
402 }
403
408 int monotonic() const
409 { return monotonic(xMin(), xMax()); }
410
427 void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const;
428
429 bool operator==(const Tabulated1DFunction<Scalar>& data) const {
430 return xValues_ == data.xValues_ &&
431 yValues_ == data.yValues_;
432 }
433
434 template <class Evaluation>
435 SegmentIndex findSegmentIndex(const Evaluation& x, bool extrapolate = false) const
436 {
437 if (!isfinite(x)) {
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)) + " .");
441 }
442
443 if (!extrapolate && !applies(x))
444 throw std::logic_error("Trying to evaluate a tabulated function outside of its range");
445
446 // we need at least two sampling points!
447 if (numSamples() < 2) {
448 throw std::logic_error("We need at least two sampling points to "
449 "do interpolation/extrapolation, "
450 "and the table only contains " +
451 std::to_string(numSamples()) +
452 " sampling points");
453 }
454
455 if (x <= xValues_[1])
456 return SegmentIndex{0};
457 else if (x >= xValues_[xValues_.size() - 2])
458 return SegmentIndex{xValues_.size() - 2};
459 else {
460 // bisection
461 size_t lowerIdx = 1;
462 size_t upperIdx = xValues_.size() - 2;
463 while (lowerIdx + 1 < upperIdx) {
464 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
465 if (x < xValues_[pivotIdx])
466 upperIdx = pivotIdx;
467 else
468 lowerIdx = pivotIdx;
469 }
470
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 " +
478 std::to_string(numSamples()) +
479 ",\nand the end values of the found segment are " +
480 std::to_string(xValues_[lowerIdx]) +
481 " and " +
482 std::to_string(xValues_[lowerIdx + 1]) +
483 ", respectively.\n";
484 msg += "Outputting the problematic table for more information "
485 "(with *** marking the found segment):";
486 for (size_t i = 0; i < numSamples(); ++i) {
487 if (i % 10 == 0)
488 msg += "\n";
489 if (i == lowerIdx)
490 msg += " ***";
491 msg += " " + std::to_string(xValues_[i]);
492 if (i == lowerIdx + 1)
493 msg += " ***";
494 }
495 msg += "\n";
496 OpmLog::debug(msg);
497 throw std::runtime_error(msg);
498 }
499 return SegmentIndex{lowerIdx};
500 }
501 }
502
503private:
504 template <class Evaluation>
505 Evaluation evalDerivative_(const Evaluation& x, size_t segIdx) const
506 {
507
508 Scalar x0 = xValues_[segIdx];
509 Scalar x1 = xValues_[segIdx + 1];
510
511 Scalar y0 = yValues_[segIdx];
512 Scalar y1 = yValues_[segIdx + 1];
513
514 Evaluation ret = blank(x);
515 ret = (y1 - y0)/(x1 - x0);
516 return ret;
517 }
518
519 // returns the monotonicity of a segment
520 //
521 // The return value have the following meaning:
522 //
523 // 3: function is constant within interval [x0, x1]
524 // 1: function is monotonously increasing in the specified interval
525 // 0: function is not monotonic in the specified interval
526 // -1: function is monotonously decreasing in the specified interval
527 int updateMonotonicity_(size_t i, int& r) const
528 {
529 if (yValues_[i] < yValues_[i + 1]) {
530 // monotonically increasing?
531 if (r == 3 || r == 1)
532 r = 1;
533 else
534 r = 0;
535 return 1;
536 }
537 else if (yValues_[i] > yValues_[i + 1]) {
538 // monotonically decreasing?
539 if (r == 3 || r == -1)
540 r = -1;
541 else
542 r = 0;
543 return -1;
544 }
545
546 return 3;
547 }
548
552 struct ComparatorX_
553 {
554 ComparatorX_(const std::vector<Scalar>& x)
555 : x_(x)
556 {}
557
558 bool operator ()(size_t idxA, size_t idxB) const
559 { return x_.at(idxA) < x_.at(idxB); }
560
561 const std::vector<Scalar>& x_;
562 };
563
567 void sortInput_()
568 {
569 size_t n = numSamples();
570
571 // create a vector containing 0...n-1
572 std::vector<unsigned> idxVector(n);
573 for (unsigned i = 0; i < n; ++i)
574 idxVector[i] = i;
575
576 // sort the indices according to the x values of the sample
577 // points
578 ComparatorX_ cmp(xValues_);
579 std::sort(idxVector.begin(), idxVector.end(), cmp);
580
581 // reorder the sample points
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]];
586 }
587 xValues_ = tmpX;
588 yValues_ = tmpY;
589 }
590
595 void reverseSamplingPoints_()
596 {
597 // reverse the arrays
598 size_t n = numSamples();
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]);
602 }
603 }
604
608 void resizeArrays_(size_t nSamples)
609 {
610 xValues_.resize(nSamples);
611 yValues_.resize(nSamples);
612 }
613
614 std::vector<Scalar> xValues_;
615 std::vector<Scalar> yValues_;
616};
617
618} // namespace Opm
619
620#endif
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