Merge pull request #317 from GitPaean/adding_2D_cartisian_grid

Adding regular 2d cartesian interpolation for the injectivity simulation.

opm/material/common/Tabulated1DFunction.hpp

@@ -470,15 +470,17 @@ class Tabulated1DFunction
     }
 
     template <class Evaluation>
-    Evaluation evalDerivative_(const Evaluation& x OPM_UNUSED, size_t segIdx) const
+    Evaluation evalDerivative_(const Evaluation& x, size_t segIdx) const
     {
         Scalar x0 = xValues_[segIdx];
         Scalar x1 = xValues_[segIdx + 1];
 
         Scalar y0 = yValues_[segIdx];
         Scalar y1 = yValues_[segIdx + 1];
 
-        return (y1 - y0)/(x1 - x0);
+        Evaluation ret = blank(x);
+        ret = (y1 - y0)/(x1 - x0);
+        return ret;
     }
 
     // returns the monotonicity of a segment

opm/material/common/XYTabulated2DFunction.hpp

@@ -0,0 +1,271 @@
+// -*- 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.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::XYTabulated2DFunction.hpp
+ */
+#ifndef OPM_XY_TABULATED_2D_FUNCTION_HPP
+#define OPM_XY_TABULATED_2D_FUNCTION_HPP
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Exceptions.hpp>
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/MathToolbox.hpp>
+
+#include <vector>
+#include <limits>
+#include <sstream>
+#include <cassert>
+#include <algorithm>
+
+namespace Opm {
+
+/*!
+ * \brief Implements a scalar function that depends on two variables, which are sampled
+ *        in X and Y directions with sampling points, respectively.
+ *        The table might be able to be extrapolated in certian directions.
+ */
+
+template <class Scalar>
+class XYTabulated2DFunction
+{
+
+public:
+    XYTabulated2DFunction()
+    { }
+
+    XYTabulated2DFunction(const std::vector<Scalar>& x_pos,
+                          const std::vector<Scalar>& y_pos,
+                          const std::vector<std::vector<Scalar> >& data,
+                          const bool x_extrapolate = false,
+                          const bool y_extrapolate = false)
+       : xPos_(x_pos)
+       , yPos_(y_pos)
+       , samples_(data)
+       , xExtrapolate_(x_extrapolate)
+       , yExtrapolate_(y_extrapolate)
+    {
+        // make sure the size is correct
+        if (numX() != samples_.size()) {
+            throw std::runtime_error("numX() is not equal to the number of rows of the sample points");
+        }
+
+        for (unsigned ix = 0; ix < numX(); ++ix) {
+            if (samples_[ix].size() != numY()) {
+                std::ostringstream oss;
+                oss << "the " << ix << "th row of the sample points has different number of data from numY() ";
+                throw std::runtime_error(oss.str());
+            }
+        }
+    }
+
+    /*!
+     * \brief Evaluate the function at a given (x,y) position.
+     *
+     * If this method is called for a value outside of the tabulated
+     * range, and extrapolation is not allowed in the corresponding direction,
+     * a \c Opm::NumericalIssue exception is thrown.
+     */
+    template <typename DataType>
+    void eval(const DataType& x, const DataType& y, DataType& result) const
+    {
+        if ( (!xExtrapolate_ && !appliesX(x)) ||
+             (!yExtrapolate_ && !appliesY(y)) ) {
+            std::ostringstream oss;
+            oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
+            throw NumericalIssue(oss.str());
+        };
+
+        // bi-linear interpolation: first, calculate the x and y indices in the lookup
+        // table ...
+        const unsigned i = xSegmentIndex(x);
+        const unsigned j = ySegmentIndex(y);
+
+        // bi-linear interpolation / extrapolation
+        const DataType alpha = xToAlpha(x, i);
+        const DataType beta = yToBeta(y, j);
+
+        const DataType s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta;
+        const DataType s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta;
+
+        Valgrind::CheckDefined(s1);
+        Valgrind::CheckDefined(s2);
+
+        // ... and combine them using the x position
+        result = s1 * (1.0 - alpha) + s2 * alpha;
+        Valgrind::CheckDefined(result);
+    }
+
+private:
+    // the sampling points in the x-drection
+    std::vector<Scalar> xPos_;
+    // the sampling points in the y-drection
+    std::vector<Scalar> yPos_;
+    // data at the sampling points
+    std::vector<std::vector<Scalar> > samples_;
+
+    bool xExtrapolate_ = false;
+    bool yExtrapolate_ = false;
+
+    /*!
+     * \brief Returns the minimum of the X coordinate of the sampling points.
+     */
+    Scalar xMin() const
+    { return xPos_.front(); }
+
+    /*!
+     * \brief Returns the maximum of the X coordinate of the sampling points.
+     */
+    Scalar xMax() const
+    { return xPos_.back(); }
+
+    /*!
+     * \brief Returns the minimum of the Y coordinate of the sampling points.
+     */
+    Scalar yMin() const
+    { return yPos_.front(); }
+
+    /*!
+     * \brief Returns the maximum of the Y coordinate of the sampling points.
+     */
+    Scalar yMax() const
+    { return yPos_.back(); }
+
+    /*!
+     * \brief Returns the value of a sampling point.
+     */
+    Scalar valueAt(size_t i, size_t j) const
+    { return samples_[i][j]; }
+
+    /*!
+     * \brief Returns true if a coordinate lies in the tabulated range on the x direction
+     */
+    template <class Evaluation>
+    bool appliesX(const Evaluation& x) const
+    {
+        if (x < xMin() || xMax() < x) {
+            return false;
+        } else {
+            return true;
+        }
+    }
+
+    /*!
+     * \brief Returns true if a coordinate lies in the tabulated range on the y direction
+     */
+    template <class Evaluation>
+    bool appliesY(const Evaluation& y) const
+    {
+        if (y < yMin() || yMax() < y) {
+            return false;
+        } else {
+            return true;
+        }
+    }
+
+    /*!
+     * \brief Returns the number of sampling points in X direction.
+     */
+    size_t numX() const
+    { return xPos_.size(); }
+
+    /*!
+     * \brief Returns the number of sampling points in Y direction.
+     */
+    size_t numY() const
+    { return yPos_.size(); }
+
+    /*!
+     * \brief Return the interval index of a given position on the x-axis.
+     */
+    template <class Evaluation>
+    unsigned xSegmentIndex(const Evaluation& x) const
+    {
+        assert(xExtrapolate_ || appliesX(x) );
+
+        return segmentIndex(x, xPos_);
+    }
+
+    /*!
+     * \brief Return the interval index of a given position on the y-axis.
+     */
+    template <class Evaluation>
+    unsigned ySegmentIndex(const Evaluation& y) const
+    {
+        assert(yExtrapolate_ || appliesY(y) );
+
+        return segmentIndex(y, yPos_);
+    }
+
+
+    template <class Evaluation>
+    static unsigned segmentIndex(const Evaluation& y, const std::vector<Scalar>& pos)
+    {
+
+        const unsigned num_pos = pos.size();
+        assert(num_pos >= 2);
+
+        if (y <= pos.front() || num_pos == 2) {
+            return 0;
+        } else if (y >= pos.back()) {
+            return num_pos - 2;
+        } else {
+            assert(num_pos >= 3);
+
+            return --std::lower_bound(pos.begin(), pos.end(), y) - pos.begin();
+        }
+    }
+
+    /*!
+     * \brief Return the relative position of an x value in an intervall
+     *
+     * The returned value can be larger than 1 or smaller than zero if it is outside of
+     * the range of the segment. In particular this happens for the extrapolation case.
+     */
+    template <class Evaluation>
+    Evaluation xToAlpha(const Evaluation& x, unsigned xSegmentIdx) const
+    {
+        Scalar x1 = xPos_[xSegmentIdx];
+        Scalar x2 = xPos_[xSegmentIdx + 1];
+        return (x - x1)/(x2 - x1);
+    }
+
+    /*!
+     * \brief Return the relative position of an y value in an interval
+     *
+     * The returned value can be larger than 1 or smaller than zero if it is outside of
+     * the range of the segment. In particular this happens for the extrapolation case.
+     */
+    template <class Evaluation>
+    Evaluation yToBeta(const Evaluation& y, unsigned ySegmentIdx) const
+    {
+        Scalar y1 = yPos_[ySegmentIdx];
+        Scalar y2 = yPos_[ySegmentIdx + 1];
+        return (y - y1)/(y2 - y1);
+    }
+
+};
+} // namespace Opm
+
+#endif

opm/material/densead/DynamicEvaluation.hpp

@@ -361,7 +361,7 @@ class Evaluation<ValueT, DynamicSize>
     // negation (unary minus) operator
     Evaluation operator-() const
     {
-        Evaluation result;
+        Evaluation result(*this);
 
         // set value and derivatives to negative
         for (int i = 0; i < length_(); ++i)

opm/material/densead/Evaluation.hpp

@@ -538,9 +538,7 @@ Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation
 template <class RhsValueType, class ValueType, int numVars>
 Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
 {
-    Evaluation<ValueType, numVars> result(a);
-    result -= b;
-    return result;
+    return -(b - a);
 }
 
 template <class RhsValueType, class ValueType, int numVars>