refactor the 2D tabulation classes
- the StaticTabulated2DFunction class and the base class (Tabulated2DFunction) are gone - the DynamicTabulated2DFunction class has been renamed to UniformTabulated2DFunction - a new class called UniformXTabulated2DFunction has been introduced. Like UniformTabulated2DFunction, it assumes uniform intervalls of the sampling points in X direction, but in contrast to UniformTabulated2DFunction, the Y locations of the sampling points can be set freely (as long as they are specified in increasing order for each x value) - add a unit test for the two tabulation classes
This commit is contained in:
Andreas Lauser
parent
a7126a7a13
commit
8f52b80448
@ -1,109 +0,0 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::StaticTabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_STATIC_TABULATED_2D_FUNCTION_HPP
|
||||
#define OPM_STATIC_TABULATED_2D_FUNCTION_HPP
|
||||
|
||||
#include <opm/core/utility/Exceptions.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/material/Tabulated2dFunction.hpp>
|
||||
|
||||
#include <assert.h>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \copydoc Opm::Tabulated2DFunction
|
||||
*
|
||||
* This class can be used when the sampling points are calculated at
|
||||
* compile time.
|
||||
*/
|
||||
template <class Traits>
|
||||
class StaticTabulated2DFunction
|
||||
: public Tabulated2DFunction<typename Traits::Scalar, StaticTabulated2DFunction<Traits> >
|
||||
{
|
||||
typedef typename Traits::Scalar Scalar;
|
||||
|
||||
public:
|
||||
StaticTabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the X coordinate of the sampling points.
|
||||
*/
|
||||
Scalar xMin() const
|
||||
{ return Traits::xMin; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the X coordinate of the sampling points.
|
||||
*/
|
||||
Scalar xMax() const
|
||||
{ return Traits::xMax; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the Y coordinate of the sampling points.
|
||||
*/
|
||||
Scalar yMin() const
|
||||
{ return Traits::yMin; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the Y coordinate of the sampling points.
|
||||
*/
|
||||
Scalar yMax() const
|
||||
{ return Traits::yMax; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in X direction.
|
||||
*/
|
||||
int numX() const
|
||||
{ return Traits::numX; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in Y direction.
|
||||
*/
|
||||
int numY() const
|
||||
{ return Traits::numY; }
|
||||
|
||||
/*!
|
||||
* \brief Get the value of the sample point which is at the
|
||||
* intersection of the \f$i\f$-th interval of the x-Axis
|
||||
* and the \f$j\f$-th of the y-Axis.
|
||||
*/
|
||||
Scalar getSamplePoint(int i, int j) const
|
||||
{
|
||||
#if !defined NDEBUG
|
||||
if (i < 0 || i >= Traits::numX ||
|
||||
j < 0 || j >= Traits::numY) {
|
||||
OPM_THROW(NumericalProblem,
|
||||
"Attempt to access element ("
|
||||
<< i << ", " << j
|
||||
<< ") on a " << Traits::name << " table of size ("
|
||||
<< Traits::numX << ", " << Traits::numY
|
||||
<< ")\n");
|
||||
};
|
||||
#endif
|
||||
return Traits::vals[i][j];
|
||||
}
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@ -1,144 +0,0 @@
|
||||
/*
|
||||
Copyright (C) 2012-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/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::Tabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_TABULATED_2D_FUNCTION_HPP
|
||||
#define OPM_TABULATED_2D_FUNCTION_HPP
|
||||
|
||||
#include <opm/core/utility/Exceptions.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <assert.h>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \brief A generic class that represents tabulated 2 dimensional functions
|
||||
*
|
||||
* This class can be used to tabulate a two dimensional function
|
||||
* \f$f(x, y)\f$ over the range \f$[x_{min}, x_{max}] \times [y_{min},
|
||||
* y_{max}]\f$. For this, the ranges of the \f$x\f$ and \f$y\f$ axes are
|
||||
* divided into \f$m\f$ and \f$n\f$ sub-intervals and the values of
|
||||
* \f$f(x_i, y_j)\f$ need to be provided. Here, \f$x_i\f$ and
|
||||
* \f$y_j\f$ are the largest positions of the \f$i\f$-th and
|
||||
* \f$j\f$-th intervall. Between these sampling points this tabulation
|
||||
* class uses linear interpolation.
|
||||
*
|
||||
* If the class is queried for a value outside of the tabulated range,
|
||||
* a \c Opm::NumericalProblem exception is thrown.
|
||||
*/
|
||||
template <class Scalar, class Implementation>
|
||||
class Tabulated2DFunction
|
||||
{
|
||||
public:
|
||||
Tabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the x-axis of the i-th interval.
|
||||
*/
|
||||
Scalar iToX(int i) const
|
||||
{
|
||||
assert(0 <= i && i < asImp_().numX());
|
||||
|
||||
return asImp_().xMin() + i*(asImp_().xMax() - asImp_().xMin())/(asImp_().numX() - 1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the y-axis of the j-th interval.
|
||||
*/
|
||||
Scalar jToY(int j) const
|
||||
{
|
||||
assert(0 <= j && j < asImp_().numY());
|
||||
|
||||
return asImp_().yMin() + j*(asImp_().yMax() - asImp_().yMin())/(asImp_().numY() - 1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the x-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the x value between the i-th and the (i+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar xToI(Scalar x) const
|
||||
{ return (x - asImp_().xMin())/(asImp_().xMax() - asImp_().xMin())*(asImp_().numX() - 1); }
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the y-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the y value between the j-th and the (j+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar yToJ(Scalar y) const
|
||||
{ return (y - asImp_().yMin())/(asImp_().yMax() - asImp_().yMin())*(asImp_().numY() - 1); }
|
||||
|
||||
/*!
|
||||
* \brief Returns true iff a coordinate lies in the tabulated range
|
||||
*/
|
||||
bool applies(Scalar x, Scalar y) const
|
||||
{ return asImp_().xMin() <= x && x <= asImp_().xMax() && asImp_().yMin() <= y && y <= asImp_().yMax(); }
|
||||
|
||||
/*!
|
||||
* \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.
|
||||
*/
|
||||
Scalar eval(Scalar x, Scalar y) const
|
||||
{
|
||||
#ifndef NDEBUG
|
||||
if (!applies(x,y))
|
||||
{
|
||||
OPM_THROW(NumericalProblem,
|
||||
"Attempt to get tabulated value for ("
|
||||
<< x << ", " << y
|
||||
<< ") on a table of extend "
|
||||
<< asImp_().xMin() << " to " << asImp_().xMax() << " times "
|
||||
<< asImp_().yMin() << " to " << asImp_().yMax());
|
||||
};
|
||||
#endif
|
||||
|
||||
Scalar alpha = xToI(x);
|
||||
Scalar beta = yToJ(y);
|
||||
|
||||
int i = std::max(0, std::min(asImp_().numX(), static_cast<int>(alpha)));
|
||||
int j = std::max(0, std::min(asImp_().numY(), static_cast<int>(beta)));
|
||||
|
||||
alpha -= i;
|
||||
beta -= j;
|
||||
|
||||
// bi-linear interpolation
|
||||
Scalar s1 = asImp_().getSamplePoint(i, j)*(1.0 - alpha) + asImp_().getSamplePoint(i + 1, j)*alpha;
|
||||
Scalar s2 = asImp_().getSamplePoint(i, j + 1)*(1.0 - alpha) + asImp_().getSamplePoint(i + 1, j + 1)*alpha;
|
||||
return s1*(1.0 - beta) + s2*beta;
|
||||
}
|
||||
|
||||
private:
|
||||
const Implementation &asImp_() const
|
||||
{ return *static_cast<const Implementation*>(this); }
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@ -19,14 +19,13 @@
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::DynamicTabulated2DFunction
|
||||
* \copydoc Opm::UniformTabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_DYNAMIC_TABULATED_2D_FUNCTION_HPP
|
||||
#define OPM_DYNAMIC_TABULATED_2D_FUNCTION_HPP
|
||||
#ifndef OPM_UNIFORM_TABULATED_2D_FUNCTION_HPP
|
||||
#define OPM_UNIFORM_TABULATED_2D_FUNCTION_HPP
|
||||
|
||||
#include <opm/core/utility/Exceptions.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/material/Tabulated2dFunction.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
@ -35,24 +34,24 @@
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \copydoc Opm::Tabulated2DFunction
|
||||
* \brief Implements a scalar function that depends on two variables and which is sampled
|
||||
* on an uniform X-Y grid.
|
||||
*
|
||||
* This class can be used when the sampling points are calculated at
|
||||
* run time.
|
||||
*/
|
||||
template <class Scalar>
|
||||
class DynamicTabulated2DFunction
|
||||
: public Tabulated2DFunction<Scalar, DynamicTabulated2DFunction<Scalar> >
|
||||
class UniformTabulated2DFunction
|
||||
{
|
||||
public:
|
||||
DynamicTabulated2DFunction()
|
||||
UniformTabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \brief Constructor where the tabulation parameters are already
|
||||
* provided.
|
||||
*/
|
||||
DynamicTabulated2DFunction(Scalar xMin, Scalar xMax, int m,
|
||||
UniformTabulated2DFunction(Scalar xMin, Scalar xMax, int m,
|
||||
Scalar yMin, Scalar yMax, int n)
|
||||
{
|
||||
resize(xMin, xMax, m, yMin, yMax, n);
|
||||
@ -112,6 +111,89 @@ public:
|
||||
int numY() const
|
||||
{ return n_; }
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the x-axis of the i-th interval.
|
||||
*/
|
||||
Scalar iToX(int i) const
|
||||
{
|
||||
assert(0 <= i && i < numX());
|
||||
|
||||
return xMin() + i*(xMax() - xMin())/(numX() - 1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the y-axis of the j-th interval.
|
||||
*/
|
||||
Scalar jToY(int j) const
|
||||
{
|
||||
assert(0 <= j && j < numY());
|
||||
|
||||
return yMin() + j*(yMax() - yMin())/(numY() - 1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the x-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the x value between the i-th and the (i+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar xToI(Scalar x) const
|
||||
{ return (x - xMin())/(xMax() - xMin())*(numX() - 1); }
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the y-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the y value between the j-th and the (j+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar yToJ(Scalar 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(); }
|
||||
|
||||
/*!
|
||||
* \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.
|
||||
*/
|
||||
Scalar eval(Scalar x, Scalar y) const
|
||||
{
|
||||
#ifndef NDEBUG
|
||||
if (!applies(x,y))
|
||||
{
|
||||
OPM_THROW(NumericalProblem,
|
||||
"Attempt to get tabulated value for ("
|
||||
<< x << ", " << y
|
||||
<< ") on a table of extend "
|
||||
<< xMin() << " to " << xMax() << " times "
|
||||
<< yMin() << " to " << yMax());
|
||||
};
|
||||
#endif
|
||||
|
||||
Scalar alpha = xToI(x);
|
||||
Scalar 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)));
|
||||
|
||||
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;
|
||||
return s1*(1.0 - beta) + s2*beta;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Get the value of the sample point which is at the
|
||||
* intersection of the \f$i\f$-th interval of the x-Axis
|
||||
322
opm/material/UniformXTabulated2DFunction.hpp
Normal file
322
opm/material/UniformXTabulated2DFunction.hpp
Normal file
@ -0,0 +1,322 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::UniformXTabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
|
||||
#define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
|
||||
|
||||
#include <opm/core/utility/Exceptions.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <limits>
|
||||
#include <tuple>
|
||||
|
||||
#include <assert.h>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \brief Implements a scalar function that depends on two variables and which is sampled
|
||||
* uniformly in the X direction, but non-uniformly on the Y axis-
|
||||
*
|
||||
* "Uniform on the X-axis" means that all Y sampling points must be located along a line
|
||||
* for this value. This class can be used when the sampling points are calculated at run
|
||||
* time.
|
||||
*/
|
||||
template <class Scalar>
|
||||
class UniformXTabulated2DFunction
|
||||
{
|
||||
typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
|
||||
|
||||
public:
|
||||
UniformXTabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \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 number of sampling points in X direction.
|
||||
*/
|
||||
int numX() const
|
||||
{ return xPos_.size(); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the Y coordinate of the sampling points for a given column.
|
||||
*/
|
||||
Scalar yMin(int i) const
|
||||
{ return std::get<1>(samples_.at(i).front()); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the Y coordinate of the sampling points for a given column.
|
||||
*/
|
||||
Scalar yMax(int i) const
|
||||
{ return std::get<1>(samples_.at(i).back()); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in Y direction a given column.
|
||||
*/
|
||||
int numY(int i) const
|
||||
{ return samples_.at(i).size(); }
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the x-axis of the i-th interval.
|
||||
*/
|
||||
Scalar iToX(int i) const
|
||||
{
|
||||
assert(0 <= i && i < numX());
|
||||
|
||||
return xPos_.at(i);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the y-axis of the j-th interval.
|
||||
*/
|
||||
Scalar jToY(int i, int 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));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the x-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the x value between the i-th and the (i+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar xToI(Scalar x, bool extrapolate = false) const
|
||||
{
|
||||
assert(extrapolate || (xMin() <= x && x <= xMax()));
|
||||
|
||||
// interval halving
|
||||
int lowerIdx = 0;
|
||||
int upperIdx = xPos_.size() - 2;
|
||||
int pivotIdx = (lowerIdx + upperIdx) / 2;
|
||||
while (lowerIdx + 1 < upperIdx) {
|
||||
if (x < xPos_[pivotIdx])
|
||||
upperIdx = pivotIdx;
|
||||
else
|
||||
lowerIdx = pivotIdx;
|
||||
|
||||
pivotIdx = (lowerIdx + upperIdx) / 2;
|
||||
}
|
||||
|
||||
Scalar x1 = xPos_[lowerIdx];
|
||||
Scalar x2 = xPos_[lowerIdx + 1];
|
||||
return lowerIdx + (x - x1)/(x2 - x1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the y-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the y value between the j-th and the (j+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar yToJ(int i, Scalar y, bool extrapolate = false) const
|
||||
{
|
||||
assert(0 <= i && i < numX());
|
||||
const auto &colSamplePoints = samples_.at(i);
|
||||
|
||||
assert(extrapolate || (yMin(i) <= y && y <= yMax(i)));
|
||||
|
||||
// interval halving
|
||||
int lowerIdx = 0;
|
||||
int upperIdx = colSamplePoints.size() - 2;
|
||||
int pivotIdx = (lowerIdx + upperIdx) / 2;
|
||||
while (lowerIdx + 1 < upperIdx) {
|
||||
if (y < std::get<1>(colSamplePoints[pivotIdx]))
|
||||
upperIdx = pivotIdx;
|
||||
else
|
||||
lowerIdx = pivotIdx;
|
||||
|
||||
pivotIdx = (lowerIdx + upperIdx) / 2;
|
||||
}
|
||||
|
||||
Scalar y1 = std::get<1>(colSamplePoints[lowerIdx]);
|
||||
Scalar y2 = std::get<1>(colSamplePoints[lowerIdx + 1]);
|
||||
return lowerIdx + (y - y1)/(y2 - y1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns true iff a coordinate lies in the tabulated range
|
||||
*/
|
||||
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(int(i));
|
||||
const auto &col2SamplePoints = samples_.at(int(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;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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.
|
||||
*/
|
||||
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);
|
||||
int i = std::max(0, std::min(numX() - 2, static_cast<int>(alpha)));
|
||||
alpha -= i;
|
||||
|
||||
Scalar beta1 = yToJ(i, y, extrapolate);
|
||||
Scalar beta2 = yToJ(i + 1, y, extrapolate);
|
||||
|
||||
int j1 = std::max(0, std::min(numY(i) - 2, static_cast<int>(beta1)));
|
||||
int j2 = std::max(0, std::min(numY(i + 1) - 2, static_cast<int>(beta2)));
|
||||
|
||||
beta1 -= j1;
|
||||
beta2 -= j2;
|
||||
|
||||
// bi-linear interpolation
|
||||
Scalar s1 = getSamplePoint(i, j1)*(1.0 - beta1) + getSamplePoint(i, j1 + 1)*beta1;
|
||||
Scalar s2 = getSamplePoint(i + 1, j2)*(1.0 - beta2) + getSamplePoint(i + 1, j2 + 1)*beta2;
|
||||
return s1*(1.0 - alpha) + s2*alpha;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Get the value of the sample point which is at the
|
||||
* intersection of the \f$i\f$-th interval of the x-Axis
|
||||
* and the \f$j\f$-th of the y-Axis.
|
||||
*/
|
||||
Scalar getSamplePoint(int i, int j) const
|
||||
{
|
||||
assert(0 <= i && i < numX());
|
||||
|
||||
const auto &colSamples = samples_[i];
|
||||
|
||||
assert(0 <= j && size_t(j) < colSamples.size());
|
||||
|
||||
return std::get<2>(colSamples[j]);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Set the x-position of a vertical line.
|
||||
*/
|
||||
void appendXPos(Scalar nextX)
|
||||
{
|
||||
#ifndef NDEBUG
|
||||
if (xPos_.size())
|
||||
assert(xPos_.back() < nextX);
|
||||
#endif
|
||||
|
||||
xPos_.push_back(nextX);
|
||||
samples_.resize(xPos_.size());
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Append a sample point.
|
||||
*/
|
||||
void appendSamplePoint(int i, Scalar y, Scalar value)
|
||||
{
|
||||
assert(0 <= i && i < numX());
|
||||
|
||||
Scalar x = iToX(i);
|
||||
samples_[i].push_back(SamplePoint(x, y, value));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Print the table for debugging purposes.
|
||||
*
|
||||
* It will produce the data in CSV format on stdout, so that it can be visualized
|
||||
* using e.g. gnuplot.
|
||||
*/
|
||||
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
|
||||
@ -10,23 +10,6 @@
|
||||
*
|
||||
* ./extractproperties 280.0 400.0 200 1e5 100e6 500
|
||||
*/
|
||||
|
||||
|
||||
class HiResDummy {
|
||||
public:
|
||||
|
||||
HiResDummy() {}
|
||||
|
||||
bool applies(double x, double y) const
|
||||
{ return false; }
|
||||
|
||||
bool hiresWeight(double x, double y) const
|
||||
{ return 0.0; }
|
||||
|
||||
bool at(double x, double y) const
|
||||
{ return 0.0; }
|
||||
};
|
||||
|
||||
struct TabulatedDensityTraits {
|
||||
typedef double Scalar;
|
||||
static const char *name;
|
||||
@ -36,7 +19,7 @@ struct TabulatedDensityTraits {
|
||||
static const int numY = 500;
|
||||
static const Scalar yMin;
|
||||
static const Scalar yMax;
|
||||
static const HiResDummy hires;
|
||||
|
||||
static const Scalar vals[numX][numY];
|
||||
};
|
||||
|
||||
@ -45,7 +28,6 @@ const double TabulatedDensityTraits::xMax = 4.000000000000000e+02;
|
||||
const double TabulatedDensityTraits::yMin = 1.000000000000000e+05;
|
||||
const double TabulatedDensityTraits::yMax = 1.000000000000000e+08;
|
||||
const char *TabulatedDensityTraits::name = "density";
|
||||
const HiResDummy TabulatedDensityTraits::hires;
|
||||
|
||||
const double TabulatedDensityTraits::vals[200][500] =
|
||||
{
|
||||
@ -20451,8 +20433,6 @@ const double TabulatedDensityTraits::vals[200][500] =
|
||||
}
|
||||
};
|
||||
|
||||
typedef Opm::StaticTabulated2DFunction< TabulatedDensityTraits > TabulatedDensity;
|
||||
|
||||
struct TabulatedEnthalpyTraits {
|
||||
typedef double Scalar;
|
||||
static const char *name;
|
||||
@ -20462,7 +20442,6 @@ struct TabulatedEnthalpyTraits {
|
||||
static const int numY = 500;
|
||||
static const Scalar yMin;
|
||||
static const Scalar yMax;
|
||||
static const HiResDummy hires;
|
||||
static const Scalar vals[numX][numY];
|
||||
};
|
||||
|
||||
@ -20471,7 +20450,6 @@ const double TabulatedEnthalpyTraits::xMax = 4.000000000000000e+02;
|
||||
const double TabulatedEnthalpyTraits::yMin = 1.000000000000000e+05;
|
||||
const double TabulatedEnthalpyTraits::yMax = 1.000000000000000e+08;
|
||||
const char *TabulatedEnthalpyTraits::name = "enthalpy";
|
||||
const HiResDummy TabulatedEnthalpyTraits::hires;
|
||||
|
||||
const double TabulatedEnthalpyTraits::vals[200][500] =
|
||||
{
|
||||
@ -40877,16 +40855,45 @@ const double TabulatedEnthalpyTraits::vals[200][500] =
|
||||
}
|
||||
};
|
||||
|
||||
typedef Opm::StaticTabulated2DFunction< TabulatedEnthalpyTraits > TabulatedEnthalpy;
|
||||
|
||||
typedef Opm::UniformTabulated2DFunction< double > TabulatedFunction;
|
||||
|
||||
// this class collects all the tabulated quantities in one convenient place
|
||||
struct CO2Tables {
|
||||
static const TabulatedEnthalpy tabulatedEnthalpy;
|
||||
static const TabulatedDensity tabulatedDensity;
|
||||
static TabulatedFunction tabulatedEnthalpy;
|
||||
static TabulatedFunction tabulatedDensity;
|
||||
static const double brineSalinity;
|
||||
};
|
||||
|
||||
const TabulatedEnthalpy CO2Tables::tabulatedEnthalpy;
|
||||
const TabulatedDensity CO2Tables::tabulatedDensity;
|
||||
TabulatedFunction CO2Tables::tabulatedEnthalpy;
|
||||
TabulatedFunction CO2Tables::tabulatedDensity;
|
||||
const double CO2Tables::brineSalinity = 1.000000000000000e-01;
|
||||
|
||||
// initialize the static tables once. this is a bit hacky in so far as it uses some
|
||||
// advanced C++ features (static initializer functions)
|
||||
int initCO2Tables_()
|
||||
{
|
||||
CO2Tables::tabulatedEnthalpy.resize(TabulatedEnthalpyTraits::xMin,
|
||||
TabulatedEnthalpyTraits::xMax,
|
||||
TabulatedEnthalpyTraits::numX,
|
||||
TabulatedEnthalpyTraits::yMin,
|
||||
TabulatedEnthalpyTraits::yMax,
|
||||
TabulatedEnthalpyTraits::numY);
|
||||
|
||||
for (int i = 0; i < TabulatedEnthalpyTraits::numX; ++i)
|
||||
for (int j = 0; j < TabulatedEnthalpyTraits::numY; ++j)
|
||||
CO2Tables::tabulatedEnthalpy.setSamplePoint(i, j, TabulatedEnthalpyTraits::vals[i][j]);
|
||||
|
||||
CO2Tables::tabulatedDensity.resize(TabulatedDensityTraits::xMin,
|
||||
TabulatedDensityTraits::xMax,
|
||||
TabulatedDensityTraits::numX,
|
||||
TabulatedDensityTraits::yMin,
|
||||
TabulatedDensityTraits::yMax,
|
||||
TabulatedDensityTraits::numY);
|
||||
|
||||
for (int i = 0; i < TabulatedDensityTraits::numX; ++i)
|
||||
for (int j = 0; j < TabulatedDensityTraits::numY; ++j)
|
||||
CO2Tables::tabulatedDensity.setSamplePoint(i, j, TabulatedDensityTraits::vals[i][j]);
|
||||
|
||||
return 0;
|
||||
}
|
||||
static int co2TablesInitDummy__ = initCO2Tables_();
|
||||
|
||||
@ -26,7 +26,7 @@
|
||||
|
||||
#include <opm/material/fluidstates/TemperatureOverlayFluidState.hpp>
|
||||
#include <opm/material/IdealGas.hpp>
|
||||
#include <opm/material/DynamicTabulated2dFunction.hpp>
|
||||
#include <opm/material/UniformTabulated2DFunction.hpp>
|
||||
|
||||
#include <opm/core/utility/Unused.hpp>
|
||||
#include <opm/core/utility/PolynomialUtils.hpp>
|
||||
@ -518,9 +518,9 @@ protected:
|
||||
{ return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
|
||||
|
||||
/*
|
||||
static DynamicTabulated2DFunction<Scalar> criticalTemperature_;
|
||||
static DynamicTabulated2DFunction<Scalar> criticalPressure_;
|
||||
static DynamicTabulated2DFunction<Scalar> criticalMolarVolume_;
|
||||
static UniformTabulated2DFunction<Scalar> criticalTemperature_;
|
||||
static UniformTabulated2DFunction<Scalar> criticalPressure_;
|
||||
static UniformTabulated2DFunction<Scalar> criticalMolarVolume_;
|
||||
*/
|
||||
};
|
||||
|
||||
@ -529,13 +529,13 @@ const Scalar PengRobinson<Scalar>::R = Opm::Constants<Scalar>::R;
|
||||
|
||||
/*
|
||||
template <class Scalar>
|
||||
DynamicTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalTemperature_;
|
||||
UniformTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalTemperature_;
|
||||
|
||||
template <class Scalar>
|
||||
DynamicTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalPressure_;
|
||||
UniformTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalPressure_;
|
||||
|
||||
template <class Scalar>
|
||||
DynamicTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalMolarVolume_;
|
||||
UniformTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalMolarVolume_;
|
||||
*/
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
371
tests/test_2dtables.cpp
Normal file
371
tests/test_2dtables.cpp
Normal file
@ -0,0 +1,371 @@
|
||||
/*
|
||||
Copyright (C) 2014 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/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \brief This is the unit test for the 2D tabulation classes.
|
||||
*
|
||||
* I.e., for the UniformTabulated2DFunction and UniformXTabulated2DFunction classes.
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/material/UniformXTabulated2DFunction.hpp>
|
||||
#include <opm/material/UniformTabulated2DFunction.hpp>
|
||||
|
||||
#include <memory>
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
|
||||
typedef double Scalar;
|
||||
|
||||
Scalar testFn1(Scalar x, Scalar y)
|
||||
{ return x; }
|
||||
|
||||
Scalar testFn2(Scalar x, Scalar y)
|
||||
{ return y; }
|
||||
|
||||
Scalar testFn3(Scalar x, Scalar y)
|
||||
{ return x*y; }
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformTabulated2DFunction<Scalar> >
|
||||
createUniformTabulatedFunction(Fn &f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
Scalar m = 50;
|
||||
|
||||
Scalar yMin = -1/2.0;
|
||||
Scalar yMax = 1/3.0;
|
||||
Scalar n = 40;
|
||||
|
||||
auto tab = std::make_shared<Opm::UniformTabulated2DFunction<Scalar>>(
|
||||
xMin, xMax, m,
|
||||
yMin, yMax, n);
|
||||
for (int i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
for (int j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n - 1) * (yMax - yMin);
|
||||
tab->setSamplePoint(i, j, f(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformXTabulated2DFunction<Scalar> >
|
||||
createUniformXTabulatedFunction(Fn &f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
Scalar m = 50;
|
||||
|
||||
Scalar yMin = -1/2.0;
|
||||
Scalar yMax = 1/3.0;
|
||||
Scalar n = 40;
|
||||
|
||||
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
|
||||
for (int i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
tab->appendXPos(x);
|
||||
for (int j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
tab->appendSamplePoint(i, y, f(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
|
||||
template <class Fn>
|
||||
std::shared_ptr<Opm::UniformXTabulated2DFunction<Scalar> >
|
||||
createUniformXTabulatedFunction2(Fn &f)
|
||||
{
|
||||
Scalar xMin = -2.0;
|
||||
Scalar xMax = 3.0;
|
||||
Scalar m = 50;
|
||||
|
||||
|
||||
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>();
|
||||
for (int i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
|
||||
tab->appendXPos(x);
|
||||
|
||||
Scalar n = i + 10;
|
||||
Scalar yMin = - (x + 1);
|
||||
Scalar yMax = (x + 1);
|
||||
|
||||
for (int j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin);
|
||||
tab->appendSamplePoint(i, y, f(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
return tab;
|
||||
}
|
||||
|
||||
template <class Fn, class Table>
|
||||
bool compareTableWithAnalyticFn(const Table &table,
|
||||
Scalar xMin,
|
||||
Scalar xMax,
|
||||
int numX,
|
||||
|
||||
Scalar yMin,
|
||||
Scalar yMax,
|
||||
int numY,
|
||||
|
||||
Fn &f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure that the tabulated function evaluates to the same thing as the analytic
|
||||
// one (modulo tolerance)
|
||||
for (int i = 1; i <= numX; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
|
||||
|
||||
for (int j = 0; j < numY; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
|
||||
if (std::abs(table->eval(x, y) - f(x, y)) > tolerance) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << table->eval(x,y) << " != " << f(x,y) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class UniformTablePtr, class UniformXTablePtr, class Fn>
|
||||
bool compareTables(const UniformTablePtr uTable,
|
||||
const UniformXTablePtr uXTable,
|
||||
Fn &f,
|
||||
Scalar tolerance = 1e-8)
|
||||
{
|
||||
// make sure the uniform and the non-uniform tables exhibit the same dimensions
|
||||
if (std::abs(uTable->xMin() - uXTable->xMin()) > 1e-8) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMin() != uXTable->xMin(): " << uTable->xMin() << " != " << uXTable->xMin() << "\n";
|
||||
return false;
|
||||
}
|
||||
if (std::abs(uTable->xMax() - uXTable->xMax()) > 1e-8) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMax() != uXTable->xMax(): " << uTable->xMax() << " != " << uXTable->xMax() << "\n";
|
||||
return false;
|
||||
}
|
||||
if (uTable->numX() != uXTable->numX()) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numX() != uXTable->numX(): " << uTable->numX() << " != " << uXTable->numX() << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
for (int i = 0; i < uTable->numX(); ++i) {
|
||||
if (std::abs(uTable->yMin() - uXTable->yMin(i)) > 1e-8) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMin() != uXTable->yMin("<<i<<"): " << uTable->yMin() << " != " << uXTable->yMin(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
if (std::abs(uTable->yMax() - uXTable->yMax(i)) > 1e-8) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMax() != uXTable->yMax("<<i<<"): " << uTable->yMax() << " != " << uXTable->yMax(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
if (uTable->numY() != uXTable->numY(i)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numY() != uXTable->numY("<<i<<"): " << uTable->numY() << " != " << uXTable->numY(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
// make sure that the x and y values are identical
|
||||
for (int i = 0; i < uTable->numX(); ++i) {
|
||||
if (std::abs(uTable->iToX(i) - uXTable->iToX(i)) > 1e-8) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->iToX("<<i<<") != uXTable->iToX("<<i<<"): " << uTable->iToX(i) << " != " << uXTable->iToX(i) << "\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
for (int j = 0; j < uTable->numY(); ++j) {
|
||||
if (std::abs(uTable->jToY(j) - uXTable->jToY(i, j)) > 1e-8) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->jToY("<<j<<") != uXTable->jToY("<<i<<","<<j<<"): " << uTable->jToY(i) << " != " << uXTable->jToY(i, j) << "\n";
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// check that the appicable range is correct. Note that due to rounding errors it is
|
||||
// undefined whether the table applies to the boundary of the tabulated domain or not
|
||||
Scalar xMin = uTable->xMin();
|
||||
Scalar yMin = uTable->yMin();
|
||||
Scalar xMax = uTable->xMax();
|
||||
Scalar yMax = uTable->yMax();
|
||||
|
||||
Scalar x = xMin - 1e-8;
|
||||
Scalar y = yMin - 1e-8;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMin - 1e-8;
|
||||
y = yMin + 1e-8;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMin + 1e-8;
|
||||
y = yMin - 1e-8;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMin + 1e-8;
|
||||
y = yMin + 1e-8;
|
||||
if (!uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (!uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax + 1e-8;
|
||||
y = yMax + 1e-8;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax - 1e-8;
|
||||
y = yMax + 1e-8;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax + 1e-8;
|
||||
y = yMax - 1e-8;
|
||||
if (uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
x = xMax - 1e-8;
|
||||
y = yMax - 1e-8;
|
||||
if (!uTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
if (!uXTable->applies(x, y)) {
|
||||
std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<<x<<","<<y<<")\n";
|
||||
return false;
|
||||
}
|
||||
|
||||
// make sure that the function values at the sampling points are identical and that
|
||||
// they correspond to the analytic function
|
||||
int m2 = uTable->numX()*5;
|
||||
int n2 = uTable->numY()*5;
|
||||
if (!compareTableWithAnalyticFn(uTable,
|
||||
xMin, xMax, m2,
|
||||
yMin, yMax, n2,
|
||||
f,
|
||||
tolerance))
|
||||
return false;
|
||||
if (!compareTableWithAnalyticFn(uXTable,
|
||||
xMin, xMax, m2,
|
||||
yMin, yMax, n2,
|
||||
f,
|
||||
tolerance))
|
||||
return false;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
auto uniformTab = createUniformTabulatedFunction(testFn1);
|
||||
auto uniformXTab = createUniformXTabulatedFunction(testFn1);
|
||||
if (!compareTables(uniformTab, uniformXTab, testFn1, /*tolerance=*/1e-12))
|
||||
return 1;
|
||||
|
||||
uniformTab = createUniformTabulatedFunction(testFn2);
|
||||
uniformXTab = createUniformXTabulatedFunction(testFn2);
|
||||
if (!compareTables(uniformTab, uniformXTab, testFn2, /*tolerance=*/1e-12))
|
||||
return 1;
|
||||
|
||||
uniformTab = createUniformTabulatedFunction(testFn3);
|
||||
uniformXTab = createUniformXTabulatedFunction(testFn3);
|
||||
if (!compareTables(uniformTab, uniformXTab, testFn3, /*tolerance=*/1e-2))
|
||||
return 1;
|
||||
|
||||
uniformXTab = createUniformXTabulatedFunction2(testFn3);
|
||||
if (!compareTableWithAnalyticFn(uniformXTab,
|
||||
-10, 10, 100,
|
||||
-10, 10, 100,
|
||||
testFn3,
|
||||
/*tolerance=*/1e-2))
|
||||
return 1;
|
||||
|
||||
// CSV output for debugging
|
||||
#if 0
|
||||
int m = 100;
|
||||
int n = 100;
|
||||
Scalar xMin = -3.0;
|
||||
Scalar xMax = 4.0;
|
||||
|
||||
Scalar yMin = -1;
|
||||
Scalar yMax = 1;
|
||||
for (int i = 0; i < m; ++i) {
|
||||
Scalar x = xMin + Scalar(i)/m * (xMax - xMin);
|
||||
|
||||
for (int j = 0; j < n; ++j) {
|
||||
Scalar y = yMin + Scalar(j)/n * (yMax - yMin);
|
||||
|
||||
std::cout << x << " "
|
||||
<< y << " "
|
||||
<< uniformXTab->eval(x,y,true) << "\n";
|
||||
}
|
||||
std::cout << "\n";
|
||||
}
|
||||
#endif
|
||||
|
||||
return 0;
|
||||
}
|
||||
@ -47,7 +47,7 @@
|
||||
|
||||
// include the tables for CO2 which are delivered with opm-material by
|
||||
// default
|
||||
#include <opm/material/StaticTabulated2dFunction.hpp>
|
||||
#include <opm/material/UniformTabulated2DFunction.hpp>
|
||||
|
||||
namespace Opm {
|
||||
namespace FluidSystemsTest {
|
||||
|
||||
Loading…
Reference in New Issue
Block a user