From 35bbae78ab0ffd190d2285fe8124925ebe4556cd Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 21 May 2015 15:33:05 +0200 Subject: [PATCH] make the linear 1D and 2D tabulation classes local-AD aware --- opm/material/common/Tabulated1DFunction.hpp | 38 +++++++++++ .../common/UniformTabulated2DFunction.hpp | 34 ++++++---- .../common/UniformXTabulated2DFunction.hpp | 63 +++++++++++++++++++ 3 files changed, 124 insertions(+), 11 deletions(-) diff --git a/opm/material/common/Tabulated1DFunction.hpp b/opm/material/common/Tabulated1DFunction.hpp index d0c2bcc28..6204f06e1 100644 --- a/opm/material/common/Tabulated1DFunction.hpp +++ b/opm/material/common/Tabulated1DFunction.hpp @@ -253,6 +253,44 @@ public: return y0 + (y1 - y0)*(x - x0)/(x1 - x0); } + /*! + * \brief Evaluate the spline at a given position. + * + * \param x The value on the abscissa where the function ought to be evaluated + * \param extrapolate If this parameter is set to true, the function will be extended + * beyond its range by straight lines, if false calling + * extrapolate for \f$ x \not [x_{min}, x_{max}]\f$ will cause a + * failed assertation. + */ + template + Evaluation eval(const Evaluation& x, bool extrapolate=false) const + { + int segIdx; + if (extrapolate && x.value < xValues_.front()) + segIdx = 0; + else if (extrapolate && x.value > xValues_.back()) + segIdx = numSamples() - 2; + else { + assert(xValues_.front() <= x.value && x.value <= xValues_.back()); + segIdx = findSegmentIndex_(x.value); + } + + Scalar x0 = xValues_[segIdx]; + Scalar x1 = xValues_[segIdx + 1]; + + Scalar y0 = yValues_[segIdx]; + Scalar y1 = yValues_[segIdx + 1]; + + Scalar m = (y1 - y0)/(x1 - x0); + + Evaluation result; + result.value = y0 + m*(x.value - x0); + for (unsigned varIdx = 0; varIdx < result.derivatives.size(); ++varIdx) + result.derivatives[varIdx] = m*x.derivatives[varIdx]; + + return result; + } + /*! * \brief Evaluate the spline's derivative at a given position. * diff --git a/opm/material/common/UniformTabulated2DFunction.hpp b/opm/material/common/UniformTabulated2DFunction.hpp index 1944677e8..1dfbde17c 100644 --- a/opm/material/common/UniformTabulated2DFunction.hpp +++ b/opm/material/common/UniformTabulated2DFunction.hpp @@ -26,6 +26,8 @@ #include #include +#include + #include @@ -139,7 +141,8 @@ public: * position of the x value between the i-th and the (i+1)-th * sample point. */ - Scalar xToI(Scalar x) const + template + Evaluation xToI(const Evaluation& x) const { return (x - xMin())/(xMax() - xMin())*(numX() - 1); } /*! @@ -150,14 +153,20 @@ public: * position of the y value between the j-th and the (j+1)-th * sample point. */ - Scalar yToJ(Scalar y) const + template + Evaluation yToJ(const Evaluation& y) const { return (y - yMin())/(yMax() - yMin())*(numY() - 1); } /*! * \brief Returns true iff a coordinate lies in the tabulated range */ - bool applies(Scalar x, Scalar y) const - { return xMin() <= x && x <= xMax() && yMin() <= y && y <= yMax(); } + template + bool applies(const Evaluation& x, const Evaluation& y) const + { + return + xMin() <= x && x <= xMax() && + yMin() <= y && y <= yMax(); + } /*! * \brief Evaluate the function at a given (x,y) position. @@ -165,8 +174,11 @@ public: * If this method is called for a value outside of the tabulated * range, a \c Opm::NumericalIssue exception is thrown. */ - Scalar eval(Scalar x, Scalar y) const + template + Evaluation eval(const Evaluation& x, const Evaluation& y) const { + typedef MathToolbox Toolbox; + #ifndef NDEBUG if (!applies(x,y)) { @@ -179,18 +191,18 @@ public: }; #endif - Scalar alpha = xToI(x); - Scalar beta = yToJ(y); + Evaluation alpha = xToI(x); + Evaluation beta = yToJ(y); - int i = std::max(0, std::min(numX() - 2, static_cast(alpha))); - int j = std::max(0, std::min(numY() - 2, static_cast(beta))); + int i = std::max(0, std::min(numX() - 2, Toolbox::value(alpha))); + int j = std::max(0, std::min(numY() - 2, Toolbox::value(beta))); alpha -= i; beta -= j; // bi-linear interpolation - Scalar s1 = getSamplePoint(i, j)*(1.0 - alpha) + getSamplePoint(i + 1, j)*alpha; - Scalar s2 = getSamplePoint(i, j + 1)*(1.0 - alpha) + getSamplePoint(i + 1, j + 1)*alpha; + const Evaluation& s1 = getSamplePoint(i, j)*(1.0 - alpha) + getSamplePoint(i + 1, j)*alpha; + const Evaluation& s2 = getSamplePoint(i, j + 1)*(1.0 - alpha) + getSamplePoint(i + 1, j + 1)*alpha; return s1*(1.0 - beta) + s2*beta; } diff --git a/opm/material/common/UniformXTabulated2DFunction.hpp b/opm/material/common/UniformXTabulated2DFunction.hpp index 92133b6b8..88966f5aa 100644 --- a/opm/material/common/UniformXTabulated2DFunction.hpp +++ b/opm/material/common/UniformXTabulated2DFunction.hpp @@ -24,6 +24,7 @@ #ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP #define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP +#include #include #include @@ -269,6 +270,68 @@ public: return s1*(1.0 - alpha) + s2*alpha; } + /*! + * \brief Evaluate the function at a given (x,y) position. + * + * If this method is called for a value outside of the tabulated + * range, a \c Opm::NumericalProblem exception is thrown. + */ + template + Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const + { +#ifndef NDEBUG + if (!extrapolate && !applies(x.value, y.value)) { + OPM_THROW(NumericalIssue, + "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)); + int i = std::max(0, std::min(numX() - 2, static_cast(alpha.value))); + alpha -= i; + + Evaluation beta1; + Evaluation beta2; + + beta1.value = yToJ(i, y.value, extrapolate); + beta2.value = yToJ(i + 1, y.value, extrapolate); + + int j1 = std::max(0, std::min(numY(i) - 2, static_cast(beta1.value))); + int j2 = std::max(0, std::min(numY(i + 1) - 2, static_cast(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; + } + /*! * \brief Set the x-position of a vertical line. *