My Project
Loading...
Searching...
No Matches
UniformXTabulated2DFunction.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*/
28#ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
29#define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
30
32
35
36#include <cassert>
37#include <cmath>
38#include <iosfwd>
39#include <limits>
40#include <tuple>
41#include <type_traits>
42#include <vector>
43
44namespace Opm {
53template <class Scalar>
55{
56public:
57 typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
58
70 LeftExtreme,
71 RightExtreme,
72 Vertical
73 };
74
75 explicit UniformXTabulated2DFunction(const InterpolationPolicy interpolationGuide = Vertical)
76 : interpolationGuide_(interpolationGuide)
77 { }
78
79 UniformXTabulated2DFunction(const std::vector<Scalar>& xPos,
80 const std::vector<Scalar>& yPos,
81 const std::vector<std::vector<SamplePoint>>& samples,
82 InterpolationPolicy interpolationGuide)
83 : samples_(samples)
84 , xPos_(xPos)
85 , yPos_(yPos)
86 , interpolationGuide_(interpolationGuide)
87 { }
88
92 Scalar xMin() const
93 { return xPos_.front(); }
94
98 Scalar xMax() const
99 { return xPos_.back(); }
100
104 Scalar xAt(size_t i) const
105 { return xPos_[i]; }
106
110 Scalar yAt(size_t i, size_t j) const
111 { return std::get<1>(samples_[i][j]); }
112
116 Scalar valueAt(size_t i, size_t j) const
117 { return std::get<2>(samples_[i][j]); }
118
122 size_t numX() const
123 { return xPos_.size(); }
124
128 Scalar yMin(unsigned i) const
129 { return std::get<1>(samples_.at(i).front()); }
130
134 Scalar yMax(unsigned i) const
135 { return std::get<1>(samples_.at(i).back()); }
136
140 size_t numY(unsigned i) const
141 { return samples_.at(i).size(); }
142
146 Scalar iToX(unsigned i) const
147 {
148 assert(i < numX());
149
150 return xPos_.at(i);
151 }
152
153 const std::vector<std::vector<SamplePoint>>& samples() const
154 {
155 return samples_;
156 }
157
158 const std::vector<Scalar>& xPos() const
159 {
160 return xPos_;
161 }
162
163 const std::vector<Scalar>& yPos() const
164 {
165 return yPos_;
166 }
167
168 InterpolationPolicy interpolationGuide() const
169 {
170 return interpolationGuide_;
171 }
172
176 Scalar jToY(unsigned i, unsigned j) const
177 {
178 assert(i < numX());
179 assert(size_t(j) < samples_[i].size());
180
181 return std::get<1>(samples_.at(i).at(j));
182 }
183
187 template <class Evaluation>
188 unsigned xSegmentIndex(const Evaluation& x,
189 [[maybe_unused]] bool extrapolate = false) const
190 {
191 assert(extrapolate || (xMin() <= x && x <= xMax()));
192
193 // we need at least two sampling points!
194 assert(xPos_.size() >= 2);
195
196 if (x <= xPos_[1])
197 return 0;
198 else if (x >= xPos_[xPos_.size() - 2])
199 return xPos_.size() - 2;
200 else {
201 assert(xPos_.size() >= 3);
202
203 // bisection
204 unsigned lowerIdx = 1;
205 unsigned upperIdx = xPos_.size() - 2;
206 while (lowerIdx + 1 < upperIdx) {
207 unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
208 if (x < xPos_[pivotIdx])
209 upperIdx = pivotIdx;
210 else
211 lowerIdx = pivotIdx;
212 }
213
214 return lowerIdx;
215 }
216 }
217
224 template <class Evaluation>
225 Evaluation xToAlpha(const Evaluation& x, unsigned segmentIdx) const
226 {
227 Scalar x1 = xPos_[segmentIdx];
228 Scalar x2 = xPos_[segmentIdx + 1];
229 return (x - x1)/(x2 - x1);
230 }
231
235 template <class Evaluation>
236 unsigned ySegmentIndex(const Evaluation& y, unsigned xSampleIdx,
237 [[maybe_unused]] bool extrapolate = false) const
238 {
239 assert(xSampleIdx < numX());
240 const auto& colSamplePoints = samples_.at(xSampleIdx);
241
242 assert(colSamplePoints.size() >= 2);
243 assert(extrapolate || (yMin(xSampleIdx) <= y && y <= yMax(xSampleIdx)));
244
245 if (y <= std::get<1>(colSamplePoints[1]))
246 return 0;
247 else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
248 return colSamplePoints.size() - 2;
249 else {
250 assert(colSamplePoints.size() >= 3);
251
252 // bisection
253 unsigned lowerIdx = 1;
254 unsigned upperIdx = colSamplePoints.size() - 2;
255 while (lowerIdx + 1 < upperIdx) {
256 unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
257 if (y < std::get<1>(colSamplePoints[pivotIdx]))
258 upperIdx = pivotIdx;
259 else
260 lowerIdx = pivotIdx;
261 }
262
263 return lowerIdx;
264 }
265 }
266
273 template <class Evaluation>
274 Evaluation yToBeta(const Evaluation& y, unsigned xSampleIdx, unsigned ySegmentIdx) const
275 {
276 assert(xSampleIdx < numX());
277 assert(ySegmentIdx < numY(xSampleIdx) - 1);
278
279 const auto& colSamplePoints = samples_.at(xSampleIdx);
280
281 Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
282 Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
283
284 return (y - y1)/(y2 - y1);
285 }
286
290 template <class Evaluation>
291 bool applies(const Evaluation& x, const Evaluation& y) const
292 {
293 if (x < xMin() || xMax() < x)
294 return false;
295
296 unsigned i = xSegmentIndex(x, /*extrapolate=*/false);
297 Scalar alpha = xToAlpha(decay<Scalar>(x), i);
298
299 const auto& col1SamplePoints = samples_.at(i);
300 const auto& col2SamplePoints = samples_.at(i + 1);
301
302 Scalar minY =
303 alpha*std::get<1>(col1SamplePoints.front()) +
304 (1 - alpha)*std::get<1>(col2SamplePoints.front());
305
306 Scalar maxY =
307 alpha*std::get<1>(col1SamplePoints.back()) +
308 (1 - alpha)*std::get<1>(col2SamplePoints.back());
309
310 return minY <= y && y <= maxY;
311 }
318 template <class Evaluation>
319 Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
320 {
321 Evaluation alpha, beta1, beta2;
322 unsigned i, j1, j2;
323 findPoints(i, j1, j2, alpha, beta1, beta2, x, y, extrapolate);
324 return eval(i, j1, j2, alpha, beta1, beta2);
325 }
326
327 template <class Evaluation>
328 void findPoints(unsigned& i,
329 unsigned& j1,
330 unsigned& j2,
331 Evaluation& alpha,
332 Evaluation& beta1,
333 Evaluation& beta2,
334 const Evaluation& x,
335 const Evaluation& y,
336 bool extrapolate) const
337 {
338#ifndef NDEBUG
339 if (!extrapolate && !applies(x, y)) {
340 if constexpr (std::is_floating_point_v<Evaluation>) {
341 throw NumericalProblem("Attempt to get undefined table value (" +
342 std::to_string(x) + ", " +
343 std::to_string(y) + ")");
344 } else {
345 throw NumericalProblem("Attempt to get undefined table value (" +
346 std::to_string(x.value()) + ", " +
347 std::to_string(y.value()) + ")");
348 }
349 };
350#endif
351
352 // bi-linear interpolation: first, calculate the x and y indices in the lookup
353 // table ...
354 i = xSegmentIndex(x, extrapolate);
355 alpha = xToAlpha(x, i);
356 // The 'shift' is used to shift the points used to interpolate within
357 // the (i) and (i+1) sets of sample points, so that when approaching
358 // the boundary of the domain given by the samples, one gets the same
359 // value as one would get by interpolating along the boundary curve
360 // itself.
361 Evaluation shift = 0.0;
362 if (interpolationGuide_ == InterpolationPolicy::Vertical) {
363 // Shift is zero, no need to reset it.
364 } else {
365 // find upper and lower y value
366 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
367 // The domain is above the boundary curve, up to y = infinity.
368 // The shift is therefore the same for all values of y.
369 shift = yPos_[i+1] - yPos_[i];
370 } else {
371 assert(interpolationGuide_ == InterpolationPolicy::RightExtreme);
372 // The domain is below the boundary curve, down to y = 0.
373 // The shift is therefore no longer the the same for all
374 // values of y, since at y = 0 the shift must be zero.
375 // The shift is computed by linear interpolation between
376 // the maximal value at the domain boundary curve, and zero.
377 shift = yPos_[i+1] - yPos_[i];
378 auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
379 if (yEnd > 0.) {
380 shift = shift * y / yEnd;
381 } else {
382 shift = 0.;
383 }
384 }
385 }
386 auto yLower = y - alpha*shift;
387 auto yUpper = y + (1-alpha)*shift;
388
389 j1 = ySegmentIndex(yLower, i, extrapolate);
390 j2 = ySegmentIndex(yUpper, i + 1, extrapolate);
391 beta1 = yToBeta(yLower, i, j1);
392 beta2 = yToBeta(yUpper, i + 1, j2);
393 }
394
395 template <class Evaluation>
396 Evaluation eval(const unsigned& i, const unsigned& j1, const unsigned& j2, const Evaluation& alpha,const Evaluation& beta1,const Evaluation& beta2) const
397 {
398 // evaluate the two function values for the same y value ...
399 const Evaluation& s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
400 const Evaluation& s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
401
402 Valgrind::CheckDefined(s1);
403 Valgrind::CheckDefined(s2);
404
405 // ... and combine them using the x position
406 const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
407 Valgrind::CheckDefined(result);
408
409 return result;
410 }
411
417 size_t appendXPos(Scalar nextX)
418 {
419 if (xPos_.empty() || xPos_.back() < nextX) {
420 xPos_.push_back(nextX);
421 yPos_.push_back(std::numeric_limits<Scalar>::lowest() / 2);
422 samples_.push_back({});
423 return xPos_.size() - 1;
424 }
425 else if (xPos_.front() > nextX) {
426 // this is slow, but so what?
427 xPos_.insert(xPos_.begin(), nextX);
428 yPos_.insert(yPos_.begin(), std::numeric_limits<Scalar>::lowest() / 2);
429 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
430 return 0;
431 }
432 throw std::invalid_argument("Sampling points should be specified either monotonically "
433 "ascending or descending.");
434 }
435
441 size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
442 {
443 assert(i < numX());
444 Scalar x = iToX(i);
445 if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
446 samples_[i].push_back(SamplePoint(x, y, value));
447 if (interpolationGuide_ == InterpolationPolicy::RightExtreme) {
448 yPos_[i] = y;
449 }
450 return samples_[i].size() - 1;
451 }
452 else if (std::get<1>(samples_[i].front()) > y) {
453 // slow, but we still don't care...
454 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
455 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
456 yPos_[i] = y;
457 }
458 return 0;
459 }
460
461 throw std::invalid_argument("Sampling points must be specified in either monotonically "
462 "ascending or descending order.");
463 }
464
471 void print(std::ostream& os) const;
472
473 bool operator==(const UniformXTabulated2DFunction<Scalar>& data) const {
474 return this->xPos() == data.xPos() &&
475 this->yPos() == data.yPos() &&
476 this->samples() == data.samples() &&
477 this->interpolationGuide() == data.interpolationGuide();
478 }
479
480private:
481 // the vector which contains the values of the sample points
482 // f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
483 // instead!
484 std::vector<std::vector<SamplePoint> > samples_;
485
486 // the position of each vertical line on the x-axis
487 std::vector<Scalar> xPos_;
488 // the position on the y-axis of the guide point
489 std::vector<Scalar> yPos_;
490 InterpolationPolicy interpolationGuide_;
491};
492} // namespace Opm
493
494#endif
Provides the OPM specific exception classes.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Some templates to wrap the valgrind client request macros.
Definition Exceptions.hpp:40
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition UniformXTabulated2DFunction.hpp:55
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:98
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
Append a sample point.
Definition UniformXTabulated2DFunction.hpp:441
Evaluation xToAlpha(const Evaluation &x, unsigned segmentIdx) const
Return the relative position of an x value in an intervall.
Definition UniformXTabulated2DFunction.hpp:225
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
Scalar xAt(size_t i) const
Returns the value of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:104
Scalar jToY(unsigned i, unsigned j) const
Return the position on the y-axis of the j-th interval.
Definition UniformXTabulated2DFunction.hpp:176
InterpolationPolicy
Indicates how interpolation will be performed.
Definition UniformXTabulated2DFunction.hpp:69
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition UniformXTabulated2DFunction.hpp:116
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
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:92
void print(std::ostream &os) const
Print the table for debugging purposes.
Definition UniformXTabulated2DFunction.cpp:32
Scalar yMin(unsigned i) const
Returns the minimum of the Y coordinate of the sampling points for a given column.
Definition UniformXTabulated2DFunction.hpp:128
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
size_t numY(unsigned i) const
Returns the number of sampling points in Y direction a given column.
Definition UniformXTabulated2DFunction.hpp:140
size_t numX() const
Returns the number of sampling points in X direction.
Definition UniformXTabulated2DFunction.hpp:122
Scalar yMax(unsigned i) const
Returns the maximum of the Y coordinate of the sampling points for a given column.
Definition UniformXTabulated2DFunction.hpp:134
Scalar iToX(unsigned i) const
Return the position on the x-axis of the i-th interval.
Definition UniformXTabulated2DFunction.hpp:146
Scalar yAt(size_t i, size_t j) const
Returns the value of the Y coordinate of a sampling point.
Definition UniformXTabulated2DFunction.hpp:110
size_t appendXPos(Scalar nextX)
Set the x-position of a vertical line.
Definition UniformXTabulated2DFunction.hpp:417
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
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true iff a coordinate lies in the tabulated range.
Definition UniformXTabulated2DFunction.hpp:291
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30