OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make the linear 1D and 2D tabulation classes local-AD aware

Browse Source
This commit is contained in:
Andreas Lauser
2015-05-21 15:33:05 +02:00
parent 2739b031b7
commit 35bbae78ab
3 changed files with 124 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
opm/material/common/Tabulated1DFunction.hpp
View File

@@ -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 <class Evaluation>
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.
*

34
opm/material/common/UniformTabulated2DFunction.hpp
View File

@@ -26,6 +26,8 @@
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/ErrorMacros.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <vector>
@@ -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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
Evaluation eval(const Evaluation& x, const Evaluation& y) const
{
typedef MathToolbox<Evaluation> 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<int>(alpha)));
int j = std::max(0, std::min(numY() - 2, static_cast<int>(beta)));
int i = std::max(0, std::min<int>(numX() - 2, Toolbox::value(alpha)));
int j = std::max(0, std::min<int>(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;
}

63
opm/material/common/UniformXTabulated2DFunction.hpp
View File

@@ -24,6 +24,7 @@
#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/ErrorMacros.hpp>
@@ -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 <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(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<int>(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<int>(beta1.value)));
int j2 = std::max(0, std::min(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;
}
/*!
* \brief Set the x-position of a vertical line.
*