diff --git a/opm/material/StaticTabulated2dFunction.hpp b/opm/material/StaticTabulated2dFunction.hpp deleted file mode 100644 index 6a6671f03..000000000 --- a/opm/material/StaticTabulated2dFunction.hpp +++ /dev/null @@ -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 . -*/ -/*! - * \file - * - * \copydoc Opm::StaticTabulated2DFunction - */ -#ifndef OPM_STATIC_TABULATED_2D_FUNCTION_HPP -#define OPM_STATIC_TABULATED_2D_FUNCTION_HPP - -#include -#include -#include - -#include - -namespace Opm { -/*! - * \copydoc Opm::Tabulated2DFunction - * - * This class can be used when the sampling points are calculated at - * compile time. - */ -template -class StaticTabulated2DFunction - : public Tabulated2DFunction > -{ - 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 diff --git a/opm/material/Tabulated2dFunction.hpp b/opm/material/Tabulated2dFunction.hpp deleted file mode 100644 index d8cd5754f..000000000 --- a/opm/material/Tabulated2dFunction.hpp +++ /dev/null @@ -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 . -*/ -/*! - * \file - * - * \copydoc Opm::Tabulated2DFunction - */ -#ifndef OPM_TABULATED_2D_FUNCTION_HPP -#define OPM_TABULATED_2D_FUNCTION_HPP - -#include -#include - -#include - -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 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(alpha))); - int j = std::max(0, std::min(asImp_().numY(), static_cast(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(this); } -}; -} // namespace Opm - -#endif diff --git a/opm/material/DynamicTabulated2dFunction.hpp b/opm/material/UniformTabulated2DFunction.hpp similarity index 54% rename from opm/material/DynamicTabulated2dFunction.hpp rename to opm/material/UniformTabulated2DFunction.hpp index 8f1a66d5a..eab75ab20 100644 --- a/opm/material/DynamicTabulated2dFunction.hpp +++ b/opm/material/UniformTabulated2DFunction.hpp @@ -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 #include -#include #include @@ -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 DynamicTabulated2DFunction - : public Tabulated2DFunction > +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(alpha))); + int j = std::max(0, std::min(numY() - 2, static_cast(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 diff --git a/opm/material/UniformXTabulated2DFunction.hpp b/opm/material/UniformXTabulated2DFunction.hpp new file mode 100644 index 000000000..5f3739822 --- /dev/null +++ b/opm/material/UniformXTabulated2DFunction.hpp @@ -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 . +*/ +/*! + * \file + * + * \copydoc Opm::UniformXTabulated2DFunction + */ +#ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP +#define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP + +#include +#include + +#include +#include +#include +#include + +#include + +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 UniformXTabulated2DFunction +{ + typedef std::tuple 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(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(beta1))); + int j2 = std::max(0, std::min(numY(i + 1) - 2, static_cast(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 > samples_; + + // the position of each vertical line on the x-axis + std::vector xPos_; +}; +} // namespace Opm + +#endif diff --git a/opm/material/components/co2tables.inc b/opm/material/components/co2tables.inc index 4e6dfee8d..f980bf501 100644 --- a/opm/material/components/co2tables.inc +++ b/opm/material/components/co2tables.inc @@ -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_(); diff --git a/opm/material/eos/PengRobinson.hpp b/opm/material/eos/PengRobinson.hpp index 4de6c380c..7e69fd66e 100644 --- a/opm/material/eos/PengRobinson.hpp +++ b/opm/material/eos/PengRobinson.hpp @@ -26,7 +26,7 @@ #include #include -#include +#include #include #include @@ -518,9 +518,9 @@ protected: { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); } /* - static DynamicTabulated2DFunction criticalTemperature_; - static DynamicTabulated2DFunction criticalPressure_; - static DynamicTabulated2DFunction criticalMolarVolume_; + static UniformTabulated2DFunction criticalTemperature_; + static UniformTabulated2DFunction criticalPressure_; + static UniformTabulated2DFunction criticalMolarVolume_; */ }; @@ -529,13 +529,13 @@ const Scalar PengRobinson::R = Opm::Constants::R; /* template -DynamicTabulated2DFunction PengRobinson::criticalTemperature_; +UniformTabulated2DFunction PengRobinson::criticalTemperature_; template -DynamicTabulated2DFunction PengRobinson::criticalPressure_; +UniformTabulated2DFunction PengRobinson::criticalPressure_; template -DynamicTabulated2DFunction PengRobinson::criticalMolarVolume_; +UniformTabulated2DFunction PengRobinson::criticalMolarVolume_; */ } // namespace Opm diff --git a/tests/test_2dtables.cpp b/tests/test_2dtables.cpp new file mode 100644 index 000000000..65f8207a1 --- /dev/null +++ b/tests/test_2dtables.cpp @@ -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 . +*/ +/*! + * \file + * + * \brief This is the unit test for the 2D tabulation classes. + * + * I.e., for the UniformTabulated2DFunction and UniformXTabulated2DFunction classes. + */ +#include "config.h" + +#include +#include + +#include +#include +#include + +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 +std::shared_ptr > +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>( + 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 +std::shared_ptr > +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>(); + 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 +std::shared_ptr > +createUniformXTabulatedFunction2(Fn &f) +{ + Scalar xMin = -2.0; + Scalar xMax = 3.0; + Scalar m = 50; + + + auto tab = std::make_shared>(); + 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 +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("<eval(x,y) << " != " << f(x,y) << "\n"; + return false; + } + } + } + + return true; +} + +template +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("<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("<yMax() << " != " << uXTable->yMax(i) << "\n"; + return false; + } + + if (uTable->numY() != uXTable->numY(i)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numY() != uXTable->numY("<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("<iToX("<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("<jToY("<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("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<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; +} diff --git a/tests/test_fluidsystems.cpp b/tests/test_fluidsystems.cpp index 6b2195595..ec44d4d31 100644 --- a/tests/test_fluidsystems.cpp +++ b/tests/test_fluidsystems.cpp @@ -47,7 +47,7 @@ // include the tables for CO2 which are delivered with opm-material by // default -#include +#include namespace Opm { namespace FluidSystemsTest {