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/material/common/Valgrind.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/MathToolbox.hpp>
 
#include <iostream>
#include <vector>
#include <limits>
#include <tuple>
#include <sstream>
#include <cassert>
#include <cmath>
 
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
    {
#ifndef NDEBUG
        if (!extrapolate && !applies(x, y)) {
            std::ostringstream oss;
            oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
            throw NumericalIssue(oss.str());
        };
#endif
 
        // bi-linear interpolation: first, calculate the x and y indices in the lookup
        // table ...
        unsigned i = xSegmentIndex(x, extrapolate);
        const Evaluation& 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;
 
        unsigned j1 = ySegmentIndex(yLower, i, extrapolate);
        unsigned j2 = ySegmentIndex(yUpper, i + 1, extrapolate);
        const Evaluation& beta1 = yToBeta(yLower, i, j1);
        const Evaluation& beta2 = yToBeta(yUpper, i + 1, j2);
 
        // 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(-1e100);
            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(), -1e100);
            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() || std::get<1>(samples_[i].back()) < y) {
            samples_[i].push_back(SamplePoint(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].insert(samples_[i].begin(), SamplePoint(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 = std::cout) const
    {
        Scalar x0 = xMin();
        Scalar x1 = xMax();
        int m = numX();
 
        Scalar y0 = 1e30;
        Scalar y1 = -1e30;
        size_t n = 0;
        for (int i = 0; i < m; ++ i) {
            y0 = std::min(y0, yMin(i));
            y1 = std::max(y1, yMax(i));
            n = std::max(n, numY(i));
        }
 
        m *= 3;
        n *= 3;
        for (int i = 0; i <= m; ++i) {
            Scalar x = x0 + (x1 - x0)*i/m;
            for (size_t j = 0; j <= n; ++j) {
                Scalar y = y0 + (y1 - y0)*j/n;
                os << x << " " << y << " " << eval(x, y) << "\n";
            }
            os << "\n";
        }
    }
 
    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