OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
991e2f570fc516ddc1a0581161a654858799daa5
opm-common/tests/test_2dtables.cpp
Markus Blatt 868531e414 Extrapolate for values not in domain (UniformTabulated2DFunction/liquidDensity). 
Extrapolate liquidDensity for unreasonable
pressure/temperature (Brine) and in UniformTabulated2DFunction.
In the current code the interpolation actually already works if the
values are outside of the tabulated region.

With this change there is now an additional function parameter
If it is true we will interpolate for every value instead of throwing
and aborting (was the case always before).
2021-09-28 16:34:42 +02:00

486 lines
17 KiB
C++
Raw Blame History

// -*- 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
*
* \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/common/UniformXTabulated2DFunction.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>
#include <opm/material/common/IntervalTabulated2DFunction.hpp>
#include <dune/common/parallel/mpihelper.hh>
#include <memory>
#include <cmath>
#include <iostream>
template <class ScalarT>
struct Test
{
typedef ScalarT Scalar;
static Scalar testFn1(Scalar x, Scalar /* y */)
{ return x; }
static Scalar testFn2(Scalar /* x */, Scalar y)
{ return y; }
static 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;
unsigned m = 50;
Scalar yMin = -1/2.0;
Scalar yMax = 1/3.0;
unsigned n = 40;
auto tab = std::make_shared<Opm::UniformTabulated2DFunction<Scalar>>(
xMin, xMax, m,
yMin, yMax, n);
for (unsigned i = 0; i < m; ++i) {
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
for (unsigned 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;
unsigned m = 50;
Scalar yMin = -1/2.0;
Scalar yMax = 1/3.0;
unsigned n = 40;
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>(Opm::UniformXTabulated2DFunction<Scalar>::InterpolationPolicy::Vertical);
for (unsigned i = 0; i < m; ++i) {
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
tab->appendXPos(x);
for (unsigned 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;
Scalar yMin = - 4.0;
Scalar yMax = 5.0;
auto tab = std::make_shared<Opm::UniformXTabulated2DFunction<Scalar>>(Opm::UniformXTabulated2DFunction<Scalar>::InterpolationPolicy::Vertical);
for (unsigned i = 0; i < m; ++i) {
Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
tab->appendXPos(x);
Scalar n = i + 10;
for (unsigned 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::IntervalTabulated2DFunction<Scalar> >
createIntervalTabulated2DFunction(Fn& f)
{
const Scalar xMin = -2.0;
const Scalar xMax = 3.0;
const unsigned m = 50;
const Scalar yMin = -1/2.0;
const Scalar yMax = 1/3.0;
const unsigned n = 40;
std::vector<Scalar> xSamples(m);
std::vector<Scalar> ySamples(n);
std::vector<std::vector<Scalar>> data(m, std::vector<Scalar>(n));
for (unsigned i = 0; i < m; ++i) {
xSamples[i] = xMin + Scalar(i)/(m - 1) * (xMax - xMin);
for (unsigned j = 0; j < n; ++j) {
ySamples[j] = yMin + Scalar(j)/(n -1) * (yMax - yMin);
data[i][j] = f(xSamples[i], ySamples[j]);
}
}
return std::make_shared<Opm::IntervalTabulated2DFunction<Scalar>>(xSamples, ySamples, data, true, true);
}
template <class Fn, class Table>
bool compareTableWithAnalyticFn(const Table& table,
Scalar xMin,
Scalar xMax,
unsigned numX,
Scalar yMin,
Scalar yMax,
unsigned 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 (unsigned i = 1; i <= numX; ++i) {
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
for (unsigned j = 0; j < numY; ++j) {
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
if (std::abs(table->eval(x, y, false) - f(x, y)) > tolerance) {
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << table->eval(x,y, false) << " != " << f(x,y) << "\n";
return false;
}
}
}
return true;
}
template <class Fn, class Table>
bool compareTableWithAnalyticFn2(const Table& table,
const Scalar xMin,
const Scalar xMax,
unsigned numX,
const Scalar yMin,
const Scalar yMax,
unsigned 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 (unsigned i = 1; i <= numX; ++i) {
Scalar x = xMin + Scalar(i)/numX*(xMax - xMin);
for (unsigned j = 0; j < numY; ++j) {
Scalar y = yMin + Scalar(j)/numY*(yMax - yMin);
Scalar result = table->eval(x, y);
if (std::abs(result - f(x, y)) > tolerance) {
std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<<x<<","<<y<<") != f("<<x<<","<<y<<"): " << result << " != " << 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()) > tolerance) {
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMin() != uXTable->xMin(): " << uTable->xMin() << " != " << uXTable->xMin() << "\n";
return false;
}
if (std::abs(uTable->xMax() - uXTable->xMax()) > tolerance) {
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 (unsigned i = 0; i < uTable->numX(); ++i) {
if (std::abs(uTable->yMin() - uXTable->yMin(i)) > tolerance) {
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)) > tolerance) {
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 (unsigned i = 0; i < uTable->numX(); ++i) {
if (std::abs(uTable->iToX(i) - uXTable->iToX(i)) > tolerance) {
std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->iToX("<<i<<") != uXTable->iToX("<<i<<"): " << uTable->iToX(i) << " != " << uXTable->iToX(i) << "\n";
return false;
}
for (unsigned j = 0; j < uTable->numY(); ++j) {
if (std::abs(uTable->jToY(j) - uXTable->jToY(i, j)) > tolerance) {
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 - tolerance;
Scalar y = yMin - tolerance;
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 - tolerance;
y = yMin + tolerance;
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 + tolerance;
y = yMin - tolerance;
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 + tolerance;
y = yMin + tolerance;
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 + tolerance;
y = yMax + tolerance;
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 - tolerance;
y = yMax + tolerance;
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 + tolerance;
y = yMax - tolerance;
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 - tolerance;
y = yMax - tolerance;
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
unsigned m2 = uTable->numX()*5;
unsigned 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;
}
};
template <class TestType>
inline int testAll(const typename TestType::Scalar tolerance = 1e-6)
{
TestType test;
auto uniformTab = test.createUniformTabulatedFunction(TestType::testFn1);
auto uniformXTab = test.createUniformXTabulatedFunction(TestType::testFn1);
if (!test.compareTables(uniformTab, uniformXTab, TestType::testFn1, tolerance))
return 1;
uniformTab = test.createUniformTabulatedFunction(TestType::testFn2);
uniformXTab = test.createUniformXTabulatedFunction(TestType::testFn2);
if (!test.compareTables(uniformTab, uniformXTab, TestType::testFn2, tolerance))
return 1;
uniformTab = test.createUniformTabulatedFunction(TestType::testFn3);
uniformXTab = test.createUniformXTabulatedFunction(TestType::testFn3);
if (!test.compareTables(uniformTab, uniformXTab, TestType::testFn3, /*tolerance=*/1e-2))
return 1;
uniformXTab = test.createUniformXTabulatedFunction2(TestType::testFn3);
if (!test.compareTableWithAnalyticFn(uniformXTab,
-2.0, 3.0, 100,
-4.0, 5.0, 100,
TestType::testFn3,
/*tolerance=*/1e-2))
return 1;
{
using ScalarType = typename TestType::Scalar;
auto xytab = test.createIntervalTabulated2DFunction(TestType::testFn1);
const ScalarType xMin = -4.0;
const ScalarType xMax = 8.0;
const unsigned m = 250;
const ScalarType yMin = -2. / 2.0;
const ScalarType yMax = 3. / 3.0;
const unsigned n = 170;
// extrapolation and interpolation involved, the tolerance needs to be bigger
const ScalarType tmpTolerance = 1000. * tolerance;
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn1, tmpTolerance))
return 1;
xytab = test.createIntervalTabulated2DFunction(TestType::testFn2);
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn2, tmpTolerance))
return 1;
xytab = test.createIntervalTabulated2DFunction(TestType::testFn3);
if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn3, tmpTolerance))
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;
}
int main(int argc, char **argv)
{
Dune::MPIHelper::instance(argc, argv);
if (testAll<Test<double> >(1e-12))
return 1;
if (testAll<Test<float> >(1e-6))
return 1;
return 0;
}
Reference in New Issue View Git Blame Copy Permalink