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:
/*
  Copyright (C) 2013 by Andreas Lauser

  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/>.
*/
#ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
#define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP

#include <opm/material/common/Valgrind.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/material/common/Unused.hpp>
#include <opm/material/common/MathToolbox.hpp>

#include <iostream>
#include <vector>
#include <limits>
#include <tuple>

#include <assert.h>

namespace Opm {

template <class Scalar>
class UniformXTabulated2DFunction
{
    typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;

public:
    UniformXTabulated2DFunction()
    { }

    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(size_t i) const
    { return std::get<1>(samples_.at(i).front()); }

    Scalar yMax(size_t i) const
    { return std::get<1>(samples_.at(i).back()); }

    size_t numY(size_t i) const
    { return samples_.at(i).size(); }

    Scalar iToX(size_t i) const
    {
        assert(0 <= i && i < numX());

        return xPos_.at(i);
    }

    Scalar jToY(size_t i, size_t j) const
    {
        assert(0 <= i && i < numX());
        assert(0 <= j && size_t(j) < samples_[i].size());

        return std::get<1>(samples_.at(i).at(j));
    }

    Scalar xToI(Scalar x, OPM_OPTIM_UNUSED bool extrapolate = false) const
    {
        assert(extrapolate || (xMin() <= x && x <= xMax()));

        // we need at least two sampling points!
        assert(xPos_.size() >= 2);

        size_t segmentIdx;
        if (x <= xPos_[1])
            segmentIdx = 0;
        else if (x >= xPos_[xPos_.size() - 2])
            segmentIdx = xPos_.size() - 2;
        else {
            // bisection
            segmentIdx = 1;
            size_t upperIdx = xPos_.size() - 2;
            while (segmentIdx + 1 < upperIdx) {
                size_t pivotIdx = (segmentIdx + upperIdx) / 2;
                if (x < xPos_[pivotIdx])
                    upperIdx = pivotIdx;
                else
                    segmentIdx = pivotIdx;
            }

            assert(xPos_[segmentIdx] <= x);
            assert(x <= xPos_[segmentIdx + 1]);
        }

        Scalar x1 = xPos_[segmentIdx];
        Scalar x2 = xPos_[segmentIdx + 1];
        return segmentIdx + (x - x1)/(x2 - x1);
    }

    Scalar yToJ(size_t i, Scalar y, OPM_OPTIM_UNUSED bool extrapolate = false) const
    {
        assert(0 <= i && i < numX());
        const auto &colSamplePoints = samples_.at(i);

        assert(extrapolate || (yMin(i) <= y && y <= yMax(i)));

        Scalar y1;
        Scalar y2;

        // interval halving
        size_t lowerIdx = 0;
        size_t upperIdx = colSamplePoints.size() - 1;
        size_t pivotIdx = (lowerIdx + upperIdx) / 2;
        while (lowerIdx + 1 < upperIdx) {
            if (y < std::get<1>(colSamplePoints[pivotIdx]))
                upperIdx = pivotIdx;
            else
                lowerIdx = pivotIdx;
            pivotIdx = (lowerIdx + upperIdx) / 2;
        }

        y1 = std::get<1>(colSamplePoints[lowerIdx]);
        y2 = std::get<1>(colSamplePoints[lowerIdx + 1]);

        assert(y1 <= y || (extrapolate && lowerIdx == 0));
        assert(y <= y2 || (extrapolate && lowerIdx == colSamplePoints.size() - 2));

        return lowerIdx + (y - y1)/(y2 - y1);
    }

    bool applies(Scalar x, Scalar y) const
    {
        if (x < xMin() || xMax() < x)
            return false;

        Scalar i = xToI(x, /*extrapolate=*/false);
        const auto &col1SamplePoints = samples_.at(unsigned(i));
        const auto &col2SamplePoints = samples_.at(unsigned(i));
        Scalar alpha = i - int(i);

        Scalar yMin =
                alpha*std::get<1>(col1SamplePoints.front()) +
                (1 - alpha)*std::get<1>(col2SamplePoints.front());

        Scalar yMax =
                alpha*std::get<1>(col1SamplePoints.back()) +
                (1 - alpha)*std::get<1>(col2SamplePoints.back());

        return yMin <= y && y <= yMax;
    }

    Scalar eval(Scalar x, Scalar y, bool extrapolate = true) const
    {
#ifndef NDEBUG
        if (!extrapolate && !applies(x,y))
        {
            OPM_THROW(NumericalProblem,
                       "Attempt to get tabulated value for ("
                       << x << ", " << y
                       << ") on table");
        };
#endif

        Scalar alpha = xToI(x, extrapolate);
        size_t i = std::max<size_t>(0, std::min(numX() - 2, static_cast<size_t>(alpha)));
        alpha -= i;

        Scalar beta1 = yToJ(i, y, extrapolate);
        Scalar beta2 = yToJ(i + 1, y, extrapolate);

        size_t j1 = std::max<size_t>(0, std::min(numY(i) - 2, static_cast<size_t>(beta1)));
        size_t j2 = std::max<size_t>(0, std::min(numY(i + 1) - 2, static_cast<size_t>(beta2)));

        beta1 -= j1;
        beta2 -= j2;

        // bi-linear interpolation
        Scalar s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
        Scalar s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
        return s1*(1.0 - alpha) + s2*alpha;
    }

    template <class Evaluation>
    Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
    {
#ifndef NDEBUG
        if (!extrapolate && !applies(x.value, y.value)) {
            OPM_THROW(NumericalProblem,
                      "Attempt to get undefined table value (" << x << ", " << y << ")");
        };
#endif

        // bi-linear interpolation: first, calculate the x and y indices in the lookup
        // table ...
        Evaluation alpha = Evaluation::createConstant(xToI(x.value, extrapolate));
        unsigned i = static_cast<unsigned>(
            std::max<int>(0, std::min<int>(static_cast<int>(numX() - 2),
                                           static_cast<int>(alpha.value))));
        alpha -= i;

        Evaluation beta1;
        Evaluation beta2;

        beta1.value = yToJ(i, y.value, extrapolate);
        beta2.value = yToJ(i + 1, y.value, extrapolate);

        unsigned j1 = static_cast<unsigned>(
            std::max(0, std::min<int>(static_cast<int>(numY(i) - 2),
                                      static_cast<int>(beta1.value))));
        unsigned j2 = static_cast<unsigned>(
            std::max(0, std::min<int>(static_cast<int>(numY(i + 1) - 2),
                                      static_cast<int>(beta2.value))));

        beta1.value -= j1;
        beta2.value -= j2;

        // set the correct derivatives of alpha and the betas
        for (unsigned varIdx = 0; varIdx < x.derivatives.size(); ++varIdx) {
            alpha.derivatives[varIdx] = x.derivatives[varIdx]/(xAt(i + 1) - xAt(i));

            beta1.derivatives[varIdx] = y.derivatives[varIdx]/(yAt(i, j1 + 1) - yAt(i, j1));
            beta2.derivatives[varIdx] = y.derivatives[varIdx]/(yAt(i + 1, j2 + 1) - yAt(i + 1, j2));
        }

        Valgrind::CheckDefined(alpha);
        Valgrind::CheckDefined(beta1);
        Valgrind::CheckDefined(beta2);

        // ... then evaluate the two function values for the same y value ...
        Evaluation s1, s2;
        s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
        s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;

        Valgrind::CheckDefined(s1);
        Valgrind::CheckDefined(s2);

        // ... and finally combine them using x the position
        Evaluation result;
        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);
            samples_.resize(xPos_.size());
            return xPos_.size() - 1;
        }
        else if (xPos_.front() > nextX) {
            // this is slow, but so what?
            xPos_.insert(xPos_.begin(), nextX);
            samples_.insert(samples_.begin(), std::vector<SamplePoint>());
            return 0;
        }
        OPM_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(0 <= i && 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));
            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));
            return 0;
        }

        OPM_THROW(std::invalid_argument,
                  "Sampling points should be specified either monotonically "
                  "ascending or descending.");
    }

    void print(std::ostream &os = std::cout) const
    {
        Scalar x0 = xMin();
        Scalar x1 = xMax();
        int m = numX();

        Scalar y0 = 1e100;
        Scalar y1 = -1e100;
        int 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 (int j = 0; j <= n; ++j) {
                Scalar y = y0 + (y1 - y0)*j/n;
                os << x << " " << y << " " << eval(x, y) << "\n";
            }
            os << "\n";
        }
    }

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_;
};
} // namespace Opm

#endif