diff --git a/opm/common/dynamictabulated2dfunction.hh b/opm/common/dynamictabulated2dfunction.hh new file mode 100644 index 000000000..98a8ee166 --- /dev/null +++ b/opm/common/dynamictabulated2dfunction.hh @@ -0,0 +1,162 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \copydoc Opm::DynamicTabulated2DFunction + */ +#ifndef OPM_DYNAMIC_TABULATED_2D_FUNCTION_HH +#define OPM_DYNAMIC_TABULATED_2D_FUNCTION_HH + +#include +#include + +#include + +#include + +namespace Opm { + +/*! + * \copydoc Opm::Tabulated2DFunction + * + * This class can be used when the sampling points are calculated at + * run time. + */ +template +class DynamicTabulated2DFunction + : public Tabulated2DFunction > +{ +public: + DynamicTabulated2DFunction() + { } + + /*! + * \brief Constructor where the tabulation parameters are already + * provided. + */ + DynamicTabulated2DFunction(Scalar xMin, Scalar xMax, int m, + Scalar yMin, Scalar yMax, int n) + { + resize(xMin, xMax, m, yMin, yMax, n); + } + + /*! + * \brief Resize the tabulation to a new range. + */ + void resize(Scalar xMin, Scalar xMax, int m, + Scalar yMin, Scalar yMax, int n) + { + samples_.resize(m*n); + + m_ = m; + n_ = n; + + xMin_ = xMin; + xMax_ = xMax; + + yMin_ = yMin; + yMax_ = yMax; + } + + /*! + * \brief Returns the minimum of the X coordinate of the sampling points. + */ + Scalar xMin() const + { return xMin_; } + + /*! + * \brief Returns the maximum of the X coordinate of the sampling points. + */ + Scalar xMax() const + { return xMax_; } + + /*! + * \brief Returns the minimum of the Y coordinate of the sampling points. + */ + Scalar yMin() const + { return yMin_; } + + /*! + * \brief Returns the maximum of the Y coordinate of the sampling points. + */ + Scalar yMax() const + { return yMax_; } + + /*! + * \brief Returns the number of sampling points in X direction. + */ + int numX() const + { return m_; } + + /*! + * \brief Returns the number of sampling points in Y direction. + */ + int numY() const + { return n_; } + + /*! + * \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 < m_); + assert(0 <= j && j < n_); + + return samples_[j*m_ + i]; + } + + /*! + * \brief Set 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. + */ + void setSamplePoint(int i, int j, Scalar value) + { + assert(0 <= i && i < m_); + assert(0 <= j && j < n_); + + samples_[j*m_ + i] = value; + } + +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 number of sample points in x direction + int m_; + + // the number of sample points in y direction + int n_; + + // the range of the tabulation on the x axis + Scalar xMin_; + Scalar xMax_; + + // the range of the tabulation on the y axis + Scalar yMin_; + Scalar yMax_; +}; +} + +#endif diff --git a/opm/common/exceptions.hh b/opm/common/exceptions.hh new file mode 100644 index 000000000..3004455fb --- /dev/null +++ b/opm/common/exceptions.hh @@ -0,0 +1,77 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2008-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief Some exceptions thrown by opm-material. + */ +#ifndef OPM_EXCEPTIONS_HH +#define OPM_EXCEPTIONS_HH + +#include + +#include + +namespace Opm { +/*! + * \ingroup Exception + * \brief Exception thrown if a fixable numerical problem occurs. + * + * (e.g. time step too big, etc.) + */ +class NumericalProblem : public Dune::Exception +{ +public: + // copy constructor + NumericalProblem(const NumericalProblem &v) + : Dune::Exception(v) + {} + + // default constructor + NumericalProblem() + {} + + // constructor with error message + NumericalProblem(const std::string &s) + { this->message(s); } +}; + +/*! + * \ingroup Exception + * \brief Exception thrown if a run-time parameter is not specified correctly. + */ +class ParameterException : public Dune::Exception +{ +public: + // copy constructor + ParameterException(const ParameterException &v) + : Dune::Exception(v) + {} + + // default constructor + ParameterException() + {} + + // constructor with error message + ParameterException(const std::string &s) + { this->message(s); } +}; + +} // namespace Opm + +#endif diff --git a/opm/common/fixedlengthspline_.hh b/opm/common/fixedlengthspline_.hh new file mode 100644 index 000000000..8e090580c --- /dev/null +++ b/opm/common/fixedlengthspline_.hh @@ -0,0 +1,412 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc Opm::FixedLengthSpline_ + */ +#ifndef OPM_FIXED_LENGTH_SPLINE_HH +#define OPM_FIXED_LENGTH_SPLINE_HH + +#include +#include +#include + +#include "splinecommon_.hh" + +namespace Opm +{ +/*! + * \brief The common code for all 3rd order polynomial splines with + * more than two sampling points. + */ +template +class FixedLengthSpline_ + : public SplineCommon_ > +{ + friend class SplineCommon_ >; + + typedef ScalarT Scalar; + typedef Dune::FieldVector Vector; + typedef Dune::BlockVector > BlockVector; + typedef Dune::BTDMatrix > BTDMatrix; + +protected: + FixedLengthSpline_() + : m_(nSamples) + {} + +public: + /*! + * \brief Returns the number of sampling points. + */ + size_t numSamples() const + { return nSamples; } + + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + // Full splines // + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using C-style arrays. + * + * This method uses separate array-like objects for the values of + * the X and Y coordinates. In this context 'array-like' means + * that an access to the members is provided via the [] + * operator. (e.g. C arrays, std::vector, std::array, etc.) Each + * array must be of size 'numberSamples' at least. Also, the number of + * sampling points must be larger than 1. + */ + template + void setXYArrays(size_t numberSamples, + const ScalarArrayX &x, + const ScalarArrayY &y, + Scalar m0, Scalar m1) + { + assert(numberSamples == numSamples()); + + for (size_t i = 0; i < numberSamples; ++i) { + xPos_[i] = x[i]; + yPos_[i] = y[i]; + } + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using STL-compatible containers. + * + * This method uses separate STL-compatible containers for the + * values of the sampling points' X and Y + * coordinates. "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y, + Scalar m0, Scalar m1) + { + assert(x.size() == y.size()); + assert(x.size() == numSamples()); + + std::copy(x.begin(), x.end(), xPos_.begin()); + std::copy(y.begin(), y.end(), yPos_.begin()); + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a C-style array. + * + * This method uses a single array of sampling points, which are + * seen as an array-like object which provides access to the X and + * Y coordinates. In this context 'array-like' means that an + * access to the members is provided via the [] operator. (e.g. C + * arrays, std::vector, std::array, etc.) The array containing + * the sampling points must be of size 'numberSamples' at least. Also, + * the number of sampling points must be larger than 1. + */ + template + void setArrayOfPoints(size_t numberSamples, + const PointArray &points, + Scalar m0, + Scalar m1) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(numberSamples == numSamples()); + + for (size_t i = 0; i < numberSamples; ++i) { + xPos_[i] = points[i][0]; + yPos_[i] = points[i][1]; + } + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a STL-compatible container of + * array-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be array-like objects storing the + * X and Y coordinates. "STL-compatible" means that the container + * provides access to iterators using the begin(), end() methods + * and also provides a size() method. Also, the number of entries + * in the X and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfPoints(const XYContainer &points, + Scalar m0, + Scalar m1) + { + assert(points.size() == numSamples()); + + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = (*it)[0]; + yPos_[i] = (*it)[1]; + } + + // make a full spline + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a STL-compatible container of + * tuple-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be tuple-like objects storing the + * X and Y coordinates. "tuple-like" means that the objects + * provide access to the x values via std::get<0>(obj) and to the + * y value via std::get<1>(obj) (e.g. std::tuple or + * std::pair). "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfTuples(const XYContainer &points, + Scalar m0, + Scalar m1) + { + assert(points.size() == numSamples()); + + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = std::get<0>(*it); + yPos_[i] = std::get<1>(*it); + } + + // make a full spline + makeFullSpline_(m0, m1); + } + + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + // Natural splines // + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + /*! + * \brief Set the sampling points natural spline using C-style arrays. + * + * This method uses separate array-like objects for the values of + * the X and Y coordinates. In this context 'array-like' means + * that an access to the members is provided via the [] + * operator. (e.g. C arrays, std::vector, std::array, etc.) Each + * array must be of size 'numberSamples' at least. Also, the number of + * sampling points must be larger than 1. + */ + template + void setXYArrays(size_t numberSamples, + const ScalarArrayX &x, + const ScalarArrayY &y) + { + assert(numberSamples == numSamples()); + + for (size_t i = 0; i < numberSamples; ++i) { + xPos_[i] = x[i]; + yPos_[i] = y[i]; + } + + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using + * STL-compatible containers. + * + * This method uses separate STL-compatible containers for the + * values of the sampling points' X and Y + * coordinates. "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y) + { + assert(x.size() == y.size()); + assert(x.size() > 1); + + setNumSamples_(x.size()); + std::copy(x.begin(), x.end(), xPos_.begin()); + std::copy(y.begin(), y.end(), yPos_.begin()); + makeNaturalSpline_(); + } + + + /*! + * \brief Set the sampling points of a natural spline using a + * C-style array. + * + * This method uses a single array of sampling points, which are + * seen as an array-like object which provides access to the X and + * Y coordinates. In this context 'array-like' means that an + * access to the members is provided via the [] operator. (e.g. C + * arrays, std::vector, std::array, etc.) The array containing + * the sampling points must be of size 'numberSamples' at least. Also, + * the number of sampling points must be larger than 1. + */ + template + void setArrayOfPoints(int numberSamples, + const PointArray &points) + { + assert(numberSamples == numSamples()); + + for (int i = 0; i < numberSamples; ++i) { + xPos_[i] = points[i][0]; + yPos_[i] = points[i][1]; + } + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * STL-compatible container of array-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be array-like objects storing the + * X and Y coordinates. "STL-compatible" means that the container + * provides access to iterators using the begin(), end() methods + * and also provides a size() method. Also, the number of entries + * in the X and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfPoints(const XYContainer &points) + { + assert(points.size() == numSamples()); + + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++ i, ++it) { + xPos_[i] = (*it)[0]; + yPos_[i] = (*it)[1]; + } + + // make a natural spline + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * STL-compatible container of tuple-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be tuple-like objects storing the + * X and Y coordinates. "tuple-like" means that the objects + * provide access to the x values via std::get<0>(obj) and to the + * y value via std::get<1>(obj) (e.g. std::tuple or + * std::pair). "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfTuples(const XYContainer &points) + { + assert(points.size() == numSamples()); + + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = std::get<0>(*it); + yPos_[i] = std::get<1>(*it); + } + + // make a natural spline + makeNaturalSpline_(); + } + +protected: + /*! + * \brief Create a full spline from the already set sampling points. + * + * Also creates temporary matrix and right hand side vector. + */ + void makeFullSpline_(Scalar m0, Scalar m1) + { + BTDMatrix M(numSamples()); + BlockVector d(numSamples()); + + // create linear system of equations + this->makeFullSystem_(M, d, m0, m1); + + // solve for the moments + M.solve(m_, d); + } + + /*! + * \brief Create a natural spline from the already set sampling points. + * + * Also creates temporary matrix and right hand side vector. + */ + void makeNaturalSpline_() + { + BTDMatrix M(numSamples()); + BlockVector d(numSamples()); + + // create linear system of equations + this->makeNaturalSystem_(M, d); + + // solve for the moments + M.solve(m_, d); + } + + /*! + * \brief Returns the x coordinate of the i-th sampling point. + */ + Scalar x_(int i) const + { return xPos_[i]; } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar y_(int i) const + { return yPos_[i]; } + + /*! + * \brief Returns the moment (i.e. second derivative) of the + * spline at the i-th sampling point. + */ + Scalar moment_(int i) const + { return m_[i]; } + + Vector xPos_; + Vector yPos_; + BlockVector m_; +}; + +} + +#endif diff --git a/opm/common/math.hh b/opm/common/math.hh new file mode 100644 index 000000000..55b2376d2 --- /dev/null +++ b/opm/common/math.hh @@ -0,0 +1,392 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief Define some often used mathematical functions + */ +#ifndef OPM_MATH_HH +#define OPM_MATH_HH + +#include +#include + +#include +#include + +namespace Opm +{ +/*! + * \ingroup Math + * \brief Calculate the harmonic mean of two scalar values. + * + * \param x The first input value + * \param y The second input value + */ +template +Scalar harmonicMean(Scalar x, Scalar y) +{ + if (x*y <= 0) + return 0; + return (2*x*y)/(x + y); +} + +/*! + * \ingroup Math + * \brief Calculate the geometric mean of two scalar values. + * + * \param x The first input value + * \param y The second input value + */ +template +Scalar geometricMean(Scalar x, Scalar y) +{ + if (x*y <= 0) + return 0; + return std::sqrt(x*y)*((x < 0)?-1:1); +} + +/*! + * \ingroup Math + * \brief Calculate the harmonic mean of a fixed-size matrix. + * + * This is done by calculating the harmonic mean for each entry + * individually. The result is stored the first argument. + * + * \param K The matrix for storing the result + * \param Ki The first input matrix + * \param Kj The second input matrix + */ +template +void harmonicMeanMatrix(Dune::FieldMatrix &K, + const Dune::FieldMatrix &Ki, + const Dune::FieldMatrix &Kj) +{ + for (int rowIdx=0; rowIdx < m; rowIdx++){ + for (int colIdx=0; colIdx< n; colIdx++){ + if (Ki[rowIdx][colIdx] != Kj[rowIdx][colIdx]) { + K[rowIdx][colIdx] = + harmonicMean(Ki[rowIdx][colIdx], + Kj[rowIdx][colIdx]); + } + else + K[rowIdx][colIdx] = Ki[rowIdx][colIdx]; + } + } +} + +/*! + * \ingroup Math + * \brief Invert a linear polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x + b \f] + * + * This method Returns the number of solutions which are in the real + * numbers, i.e. 1 except if the slope of the line is 0. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the linear term + * \param b The coefficient for the constant term + */ +template +int invertLinearPolynomial(SolContainer &sol, + Scalar a, + Scalar b) +{ + if (a == 0.0) + return 0; + + sol[0] = -b/a; + return 1; +} + +/*! + * \ingroup Math + * \brief Invert a quadratic polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x^2 + + b\;x + c \f] + * + * This method teturns the number of solutions which are in the real + * numbers. The "sol" argument contains the real roots of the parabola + * in order with the smallest root first. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the quadratic term + * \param b The coefficient for the linear term + * \param c The coefficient for the constant term + */ +template +int invertQuadraticPolynomial(SolContainer &sol, + Scalar a, + Scalar b, + Scalar c) +{ + // check for a line + if (a == 0.0) + return invertLinearPolynomial(sol, b, c); + + // discriminant + Scalar Delta = b*b - 4*a*c; + if (Delta < 0) + return 0; // no real roots + + Delta = std::sqrt(Delta); + sol[0] = (- b + Delta)/(2*a); + sol[1] = (- b - Delta)/(2*a); + + // sort the result + if (sol[0] > sol[1]) + std::swap(sol[0], sol[1]); + return 2; // two real roots +} + +//! \cond SKIP_THIS +template +void invertCubicPolynomialPostProcess_(SolContainer &sol, + int numSol, + Scalar a, + Scalar b, + Scalar c, + Scalar d) +{ + // do one Newton iteration on the analytic solution if the + // precision is increased + for (int i = 0; i < numSol; ++i) { + Scalar x = sol[i]; + Scalar fOld = d + x*(c + x*(b + x*a)); + + Scalar fPrime = c + x*(2*b + x*3*a); + if (fPrime == 0.0) + continue; + x -= fOld/fPrime; + + Scalar fNew = d + x*(c + x*(b + x*a)); + if (std::abs(fNew) < std::abs(fOld)) + sol[i] = x; + } +} +//! \endcond + +/*! + * \ingroup Math + * \brief Invert a cubic polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f] + * + * This method teturns the number of solutions which are in the real + * numbers. The "sol" argument contains the real roots of the cubic + * polynomial in order with the smallest root first. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the cubic term + * \param b The coefficient for the quadratic term + * \param c The coefficient for the linear term + * \param d The coefficient for the constant term + */ +template +int invertCubicPolynomial(SolContainer *sol, + Scalar a, + Scalar b, + Scalar c, + Scalar d) +{ + // reduces to a quadratic polynomial + if (a == 0) + return invertQuadraticPolynomial(sol, b, c, d); + + // normalize the polynomial + b /= a; + c /= a; + d /= a; + a = 1; + + // get rid of the quadratic term by subsituting x = t - b/3 + Scalar p = c - b*b/3; + Scalar q = d + (2*b*b*b - 9*b*c)/27; + + if (p != 0.0 && q != 0.0) { + // At this point + // + // t^3 + p*t + q = 0 + // + // with p != 0 and q != 0 holds. Introducing the variables u and v + // with the properties + // + // u + v = t and 3*u*v + p = 0 + // + // leads to + // + // u^3 + v^3 + q = 0 . + // + // multiplying both sides with u^3 and taking advantage of the + // fact that u*v = -p/3 leads to + // + // u^6 + q*u^3 - p^3/27 = 0 + // + // Now, substituting u^3 = w yields + // + // w^2 + q*w - p^3/27 = 0 + // + // This is a quadratic equation with the solutions + // + // w = -q/2 + sqrt(q^2/4 + p^3/27) and + // w = -q/2 - sqrt(q^2/4 + p^3/27) + // + // Since w is equivalent to u^3 it is sufficient to only look at + // one of the two cases. Then, there are still 2 cases: positive + // and negative discriminant. + Scalar wDisc = q*q/4 + p*p*p/27; + if (wDisc >= 0) { // the positive discriminant case: + // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27) + Scalar u = - q/2 + std::sqrt(wDisc); + if (u < 0) u = - std::pow(-u, 1.0/3); + else u = std::pow(u, 1.0/3); + + // at this point, u != 0 since p^3 = 0 is necessary in order + // for u = 0 to hold, so + sol[0] = u - p/(3*u) - b/3; + // does not produce a division by zero. the remaining two + // roots of u are rotated by +- 2/3*pi in the complex plane + // and thus not considered here + invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d); + return 1; + } + else { // the negative discriminant case: + Scalar uCubedRe = - q/2; + Scalar uCubedIm = std::sqrt(-wDisc); + // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27) + Scalar uAbs = std::pow(std::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3); + Scalar phi = std::atan2(uCubedIm, uCubedRe)/3; + + // calculate the length and the angle of the primitive root + + // with the definitions from above it follows that + // + // x = u - p/(3*u) - b/3 + // + // where x and u are complex numbers. Rewritten in polar form + // this is equivalent to + // + // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 . + // + // Factoring out the e^ terms and subtracting the additional + // terms, yields + // + // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3 + // + // with + // + // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) . + // + // The crucial observation is the fact that y is the conjugate + // of - x + b/3. This means that after taking advantage of the + // relation + // + // e^(i*phi) + e^(-i*phi) = 2*cos(phi) + // + // the equation + // + // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3 + // + // holds. Since |u|, p, b and cos(phi) are real numbers, it + // follows that Im(x) = - Im(x) and thus Im(x) = 0. This + // implies + // + // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 . + // + // Considering the fact that u is a cubic root, we have three + // values for phi which differ by 2/3*pi. This allows to + // calculate the three real roots of the polynomial: + for (int i = 0; i < 3; ++i) { + sol[i] = std::cos(phi)*(uAbs - p/(3*uAbs)) - b/3; + phi += 2*M_PI/3; + } + + // post process the obtained solution to increase numerical + // precision + invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d); + + // sort the result + std::sort(sol, sol + 3); + + return 3; + } + } + // Handle some (pretty unlikely) special cases required to avoid + // divisions by zero in the code above... + else if (p == 0.0 && q == 0.0) { + // t^3 = 0, i.e. triple root at t = 0 + sol[0] = sol[1] = sol[2] = 0.0 - b/3; + return 3; + } + else if (p == 0.0 && q != 0.0) { + // t^3 + q = 0, + // + // i. e. single real root at t=curt(q) + Scalar t; + if (-q > 0) t = std::pow(-q, 1./3); + else t = - std::pow(q, 1./3); + sol[0] = t - b/3; + + return 1; + } + + assert(p != 0.0 && q == 0.0); + + // t^3 + p*t = 0 = t*(t^2 + p), + // + // i. e. roots at t = 0, t^2 + p = 0 + if (p > 0) { + sol[0] = 0.0 - b/3; + return 1; // only a single real root at t=0 + } + + // two additional real roots at t = sqrt(-p) and t = -sqrt(-p) + sol[0] = -std::sqrt(-p) - b/3;; + sol[1] = 0.0 - b/3; + sol[2] = std::sqrt(-p) - b/3; + + return 3; +} + +/*! + * \brief Cross product of two vectors in three-dimensional Euclidean space + * + * \param vec1 The first vector + * \param vec2 The second vector + */ +template +Dune::FieldVector crossProduct(const Dune::FieldVector &vec1, + const Dune::FieldVector &vec2) +{ + Dune::FieldVector result; + + result[0] = vec1[1]*vec2[2] - vec1[2]*vec2[1]; + result[1] = vec1[2]*vec2[0] - vec1[0]*vec2[2]; + result[2] = vec1[0]*vec2[1] - vec1[1]*vec2[0]; + + return result; +} + +} + + +#endif diff --git a/opm/common/spline.hh b/opm/common/spline.hh new file mode 100644 index 000000000..60f6d25ed --- /dev/null +++ b/opm/common/spline.hh @@ -0,0 +1,532 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc Opm::Spline + */ +#ifndef OPM_SPLINE_HH +#define OPM_SPLINE_HH + +#include "fixedlengthspline_.hh" +#include "variablelengthspline_.hh" +#include "splinecommon_.hh" + +#include +#include + +namespace Opm +{ +/*! + * \ingroup Spline + * \brief A 3rd order polynomial spline. + + * This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling + * points \f$x_1, \dots, x_n\f$, the following conditions hold + *\f{align*}{ + s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \}\\ + s'(x_1) & = m_1 \\ + s'(x_n) & = m_n + \f} +* +* for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. Alternatively, natural +* splines are supported which are defined by +*\f{align*}{ + s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\ + s''(x_1) & = 0 \\ + s''(x_n) & = 0 +\f} + */ +template +class Spline : public FixedLengthSpline_ +{ +public: + /*! + * \brief Default constructor for a spline. + * + * To specfiy the acutal curve, use one of the set() methods. + */ + Spline() + { } + + /*! + * \brief Convenience constructor for a full spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + */ + template + Spline(const ScalarArray &x, + const ScalarArray &y) + { this->setXYArrays(numSamples, x, y); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + */ + template + Spline(const PointArray &points) + { this->setArrayOfPoints(numSamples, points); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(const ScalarArray &x, + const ScalarArray &y, + Scalar m0, + Scalar m1) + { this->setXYArrays(numSamples, x, y, m0, m1); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(const PointArray &points, + Scalar m0, + Scalar m1) + { this->setArrayOfPoints(numSamples, points, m0, m1); } +}; + +/*! + * \brief Specialization of a spline with the number of sampling + * points only known at run time. + * + * This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling + * points \f$x_1, \dots, x_n\f$, the following conditions hold + *\f{align*}{ + s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \}\\ + s'(x_1) & = m_1 \\ + s'(x_n) & = m_n + \f} + * + * for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. Alternatively, natural + * splines are supported which are defined by + *\f{align*}{ + s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\ + s''(x_1) & = 0 \\ + s''(x_n) & = 0 + \f} +*/ +template +class Spline : public VariableLengthSpline_ +{ +public: + /*! + * \brief Default constructor for a spline. + * + * To specfiy the acutal curve, use one of the set() methods. + */ + Spline() + { } + + /*! + * \brief Convenience constructor for a natural spline + * + * \param nSamples The number of sampling points (must be > 2) + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + */ + template + Spline(int nSamples, + const ScalarArrayX &x, + const ScalarArrayY &y) + { this->setXYArrays(nSamples, x, y); } + + /*! + * \brief Convenience constructor for a natural spline + * + * \param nSamples The number of sampling points (must be > 2) + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + */ + template + Spline(int nSamples, + const PointArray &points) + { this->setArrayOfPoints(nSamples, points); } + + /*! + * \brief Convenience constructor for a natural spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method) + * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method) + */ + template + Spline(const ScalarContainer &x, + const ScalarContainer &y) + { this->setXYContainers(x, y); } + + /*! + * \brief Convenience constructor for a natural spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method) + */ + template + Spline(const PointContainer &points) + { this->setContainerOfPoints(points); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param nSamples The number of sampling points (must be >= 2) + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(int nSamples, + const ScalarArray &x, + const ScalarArray &y, + Scalar m0, + Scalar m1) + { this->setXYArrays(nSamples, x, y, m0, m1); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param nSamples The number of sampling points (must be >= 2) + * \param points An array containing the \f$x\f$ and \f$x\f$ values of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(int nSamples, + const PointArray &points, + Scalar m0, + Scalar m1) + { this->setArrayOfPoints(nSamples, points, m0, m1); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method) + * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method) + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(const ScalarContainerX &x, + const ScalarContainerY &y, + Scalar m0, + Scalar m1) + { this->setXYContainers(x, y, m0, m1); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method) + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(const PointContainer &points, + Scalar m0, + Scalar m1) + { this->setContainerOfPoints(points, m0, m1); } +}; + +/*! + * \brief Do not allow splines with zero sampling points + */ +template +class Spline +// Splines with zero sampling points do not make sense! +{ private: Spline() { } }; + +/*! + * \brief Do not allow splines with one sampling point + */ +template +class Spline +// Splines with one sampling point do not make sense! +{ private: Spline() { } }; + +/*! + * \brief Spline for two sampling points. + * + * For this type of spline there is no natural spline. + */ +template +class Spline : public SplineCommon_ > +{ + friend class SplineCommon_ >; + typedef Dune::FieldVector Vector; + typedef Dune::FieldMatrix Matrix; + +public: + Spline() + {} + + /*! + * \brief Convenience constructor for a full spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(const ScalarArrayX &x, + const ScalarArrayY &y, + Scalar m0, Scalar m1) + { setXYArrays(2, x, y, m0, m1); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + template + Spline(const PointArray &points, + Scalar m0, + Scalar m1) + { this->setArrayOfPoints(2, points, m0, m1); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param x0 The \f$x\f$ value of the first sampling point + * \param x1 The \f$x\f$ value of the second sampling point + * \param y0 The \f$y\f$ value of the first sampling point + * \param y1 The \f$y\f$ value of the second sampling point + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + */ + Spline(Scalar x0, Scalar x1, + Scalar y0, Scalar y1, + Scalar m0, Scalar m1) + { + set(x0, x1, + y0, y1, + m0, m1); + } + + /*! + * \brief Returns the number of sampling points. + */ + int numSamples() const + { return 2; } + + /*! + * \brief Set the sampling points and the boundary slopes of the + * spline. + * + * \param x0 The \f$x\f$ value of the first sampling point + * \param x1 The \f$x\f$ value of the second sampling point + * \param y0 The \f$y\f$ value of the first sampling point + * \param y1 The \f$y\f$ value of the second sampling point + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + void set(Scalar x0, Scalar x1, + Scalar y0, Scalar y1, + Scalar m0, Scalar m1) + { + Matrix M(numSamples()); + Vector d; + assignXY_(x0, x1, y0, y1); + this->makeFullSystem_(M, d, m0, m1); + + // solve for the moments + M.solve(m_, d); + } + + /*! + * \brief Set the sampling points and the boundary slopes of the + * spline. + * + * \param nSamples The number of sampling points (must be >= 2) + * \param x An array containing the \f$x\f$ values of the sampling points + * \param y An array containing the \f$y\f$ values of the sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + template + void setXYArrays(int nSamples, + const ScalarContainer &x, + const ScalarContainer &y, + Scalar m0, Scalar m1) + { + assert(nSamples == 2); + set(x[0], x[1], y[0], y[1], m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of the + * spline. + * + * \param x An array containing the \f$x\f$ values of the sampling points + * \param y An array containing the \f$y\f$ values of the sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y, + Scalar m0, Scalar m1) + { + assert(x.size() == y.size()); + assert(x.size() == 2); + + Matrix M(numSamples()); + Vector d; + + typename ScalarContainerX::const_iterator xIt0 = x.begin(); + typename ScalarContainerX::const_iterator xIt1 = xIt0; + ++xIt1; + typename ScalarContainerY::const_iterator yIt0 = y.begin(); + typename ScalarContainerY::const_iterator yIt1 = yIt0; + ++yIt1; + set(*xIt0, *xIt1, *yIt0, *yIt1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of the + * spline. + * + * \param nSamples The number of sampling points (must be >= 2) + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + template + void setArrayOfPoints(int nSamples, + const PointArray &points, + Scalar m0, + Scalar m1) + { + assert(nSamples == 2); + + set(points[0][0], + points[1][0], + points[0][1], + points[1][1], + m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes from an + * STL-like container of points. + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + template + void setContainerOfPoints(const PointContainer &points, + Scalar m0, + Scalar m1) + { + assert(points.size() == 2); + + Matrix M; + Vector d; + typename PointContainer::const_iterator it0 = points.begin(); + typename PointContainer::const_iterator it1 = it0; + ++it1; + + set((*it0)[0], + (*it0)[1], + (*it1)[0], + (*it1)[1], + m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes from an + * STL-like container of tuples. + * + * \param tuples An array of \f$(x,y)\f$ tuples of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + template + void setContainerOfTuples(const TupleContainer &tuples, + Scalar m0, + Scalar m1) + { + assert(tuples.size() == 2); + + typename TupleContainer::const_iterator it0 = tuples.begin(); + typename TupleContainer::const_iterator it1 = it0; + ++it1; + + set(std::get<0>(*it0), + std::get<1>(*it0), + std::get<0>(*it1), + std::get<1>(*it1), + m0, m1); + } + +protected: + void assignXY_(Scalar x0, Scalar x1, + Scalar y0, Scalar y1) + { + if (x0 > x1) { + xPos_[0] = x1; + xPos_[1] = x0; + yPos_[0] = y1; + yPos_[1] = y0; + } + else { + xPos_[0] = x0; + xPos_[1] = x1; + yPos_[0] = y0; + yPos_[1] = y1; + } + } + + /*! + * \brief Returns the x coordinate of the i-th sampling point. + */ + Scalar x_(int i) const + { return xPos_[i]; } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar y_(int i) const + { return yPos_[i]; } + + /*! + * \brief Returns the moment (i.e. second derivative) of the + * spline at the i-th sampling point. + */ + Scalar moment_(int i) const + { return m_[i]; } + + Vector xPos_; + Vector yPos_; + Vector m_; +}; + +} + +#endif diff --git a/opm/common/splinecommon_.hh b/opm/common/splinecommon_.hh new file mode 100644 index 000000000..11f79e15b --- /dev/null +++ b/opm/common/splinecommon_.hh @@ -0,0 +1,685 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc Opm::SplineCommon_ + */ +#ifndef OPM_SPLINE_COMMON__HH +#define OPM_SPLINE_COMMON__HH + +#include "valgrind.hh" +#include "math.hh" + +#include + +#include +#include +#include +#include + +namespace Opm +{ +/*! + * \brief The common code for all 3rd order polynomial splines. + */ +template +class SplineCommon_ +{ + typedef ScalarT Scalar; + typedef ImplementationT Implementation; + + Implementation &asImp_() + { return *static_cast(this); } + const Implementation &asImp_() const + { return *static_cast(this); } + +public: + /*! + * \brief Return true iff the given x is in range [x1, xn]. + */ + bool applies(Scalar x) const + { + return x_(0) <= x && x <= x_(numSamples_() - 1); + } + + /*! + * \brief Return the x value of the leftmost sampling point. + */ + Scalar xMin() const + { return x_(0); } + + /*! + * \brief Return the x value of the rightmost sampling point. + */ + Scalar xMax() const + { return x_(numSamples_() - 1); } + + /*! + * \brief Prints k tuples of the format (x, y, dx/dy, isMonotonic) + * to stdout. + * + * If the spline does not apply for parts of [x0, x1] it is + * extrapolated using a straight line. The result can be inspected + * using the following commands: + * + ----------- snip ----------- + ./yourProgramm > spline.csv + gnuplot + + gnuplot> plot "spline.csv" using 1:2 w l ti "Curve", \ + "spline.csv" using 1:3 w l ti "Derivative", \ + "spline.csv" using 1:4 w p ti "Monotonic" + ----------- snap ----------- + */ + void printCSV(Scalar xi0, Scalar xi1, int k) const + { + Scalar x0 = std::min(xi0, xi1); + Scalar x1 = std::max(xi0, xi1); + const int n = numSamples_() - 1; + for (int i = 0; i <= k; ++i) { + double x = i*(x1 - x0)/k + x0; + double x_p1 = x + (x1 - x0)/k; + double y; + double dy_dx; + double mono = 1; + if (!applies(x)) { + if (x < 0) { + dy_dx = evalDerivative(x_(0)); + y = (x - x_(0))*dy_dx + y_(0); + mono = (dy_dx>0)?1:-1; + } + else if (x > x_(n)) { + dy_dx = evalDerivative(x_(n)); + y = (x - x_(n))*dy_dx + y_(n); + mono = (dy_dx>0)?1:-1; + } + else { + DUNE_THROW(Dune::InvalidStateException, + "The sampling points given to a spline must be sorted by their x value!"); + } + } + else { + y = eval(x); + dy_dx = evalDerivative(x); + mono = monotonic(std::max(x_(0), x), std::min(x_(n), x_p1)); + } + + std::cout << x << " " << y << " " << dy_dx << " " << mono << "\n"; + } + } + + /*! + * \brief Evaluate the spline at a given position. + * + * \param x The value on the abscissa where the spline ought to be + * evaluated + * \param extrapolate If this parameter is set to true, the spline + * will be extended beyond its range by + * straight lines, if false calling extrapolate + * for \f$ x \not [x_{min}, x_{max}]\f$ will + * cause a failed assertation. + */ + Scalar eval(Scalar x, bool extrapolate=false) const + { + assert(extrapolate || applies(x)); + + // handle extrapolation + if (extrapolate) { + if (x < xMin()) { + Scalar m = evalDerivative(xMin(), /*segmentIdx=*/0); + Scalar y0 = y_(0); + return y0 + m*(x - xMin()); + } + else if (x > xMax()) { + Scalar m = evalDerivative(xMax(), /*segmentIdx=*/numSamples_()-2); + Scalar y0 = y_(numSamples_() - 1); + return y0 + m*(x - xMax()); + } + } + + return eval_(x, segmentIdx_(x)); + } + + /*! + * \brief Evaluate the spline's derivative at a given position. + * + * \param x The value on the abscissa where the spline's + * derivative ought to be evaluated + * + * \param extrapolate If this parameter is set to true, the spline + * will be extended beyond its range by + * straight lines, if false calling extrapolate + * for \f$ x \not [x_{min}, x_{max}]\f$ will + * cause a failed assertation. + + */ + Scalar evalDerivative(Scalar x, bool extrapolate=false) const + { + assert(extrapolate || applies(x)); + if (extrapolate) { + if (x < xMin()) + evalDerivative_(xMin(), 0); + else if (x > xMax()) + evalDerivative_(xMax(), numSamples_() - 2); + } + + return evalDerivative_(x, segmentIdx_(x)); + } + + /*! + * \brief Find the intersections of the spline with a cubic + * polynomial in the whole intervall, throws + * Dune::MathError exception if there is more or less than + * one solution. + */ + Scalar intersect(Scalar a, Scalar b, Scalar c, Scalar d) const + { + return intersectIntervall(xMin(), xMax(), a, b, c, d); + } + + /*! + * \brief Find the intersections of the spline with a cubic + * polynomial in a sub-intervall of the spline, throws + * Dune::MathError exception if there is more or less than + * one solution. + */ + Scalar intersectInterval(Scalar x0, Scalar x1, + Scalar a, Scalar b, Scalar c, Scalar d) const + { + assert(applies(x0) && applies(x1)); + + Scalar tmpSol[3]; + int nSol = 0; + int iFirst = segmentIdx_(x0); + int iLast = segmentIdx_(x1); + for (int i = iFirst; i <= iLast; ++i) + { + nSol += intersectSegment_(tmpSol, i, a, b, c, d, x0, x1); + + if (nSol > 1) { + DUNE_THROW(Dune::MathError, + "Spline has more than one intersection"); //< x1) + std::swap(x0, x1); + + assert(x0 < x1); + + // corner case where the whole spline is a constant + if (moment_(0) == 0 && + moment_(1) == 0 && + y_(0) == y_(1)) + { + // actually the is monotonically increasing as well as + // monotonously decreasing + return 2; + } + + + int i = segmentIdx_(x0); + if (x_(i) <= x0 && x1 <= x_(i+1)) { + // the interval to check is completely included in a + // single segment + return monotonic_(i, x0, x1); + } + + // make sure that the segments which are completly in the + // interval [x0, x1] all exhibit the same monotonicity. + int iEnd = segmentIdx_(x1); + int r = monotonic_(i, x0, x_(1)); + for (; i < iEnd - 1; ++i) + if (r != monotonic_(i, x_(i), x_(i + 1))) + return 0; + + // check for the last segment + if (x_(iEnd) < x1 && r != monotonic_(iEnd, x_(iEnd), x1)) + { return 0; } + + return r; + } + + /*! + * \brief Same as monotonic(x0, x1), but with the entire range of the + * spline as interval. + */ + int monotonic() const + { return monotonic(xMin(), xMax()); } + +protected: + // this is an internal class, so everything is protected! + SplineCommon_() + { Valgrind::SetUndefined(asImp_()); } + + /*! + * \brief Set the sampling point vectors. + * + * This takes care that the order of the x-values is ascending, + * although the input must be ordered already! + */ + template + void assignSamplingPoints_(DestVector &destX, + DestVector &destY, + const SourceVector &srcX, + const SourceVector &srcY, + int numSamples) + { + assert(numSamples >= 2); + + // copy sample points, make sure that the first x value is + // smaller than the last one + for (int i = 0; i < numSamples; ++i) { + int idx = i; + if (srcX[0] > srcX[numSamples - 1]) + idx = numSamples - i - 1; + destX[i] = srcX[idx]; + destY[i] = srcY[idx]; + } + } + + template + void assignFromArrayList_(DestVector &destX, + DestVector &destY, + const ListIterator &srcBegin, + const ListIterator &srcEnd, + int numSamples) + { + assert(numSamples >= 2); + + // find out wether the x values are in reverse order + ListIterator it = srcBegin; + ++it; + bool reverse = false; + if ((*srcBegin)[0] > (*it)[0]) + reverse = true; + --it; + + // loop over all sampling points + for (int i = 0; it != srcEnd; ++i, ++it) { + int idx = i; + if (reverse) + idx = numSamples - i - 1; + destX[i] = (*it)[0]; + destY[i] = (*it)[1]; + } + } + + /*! + * \brief Set the sampling points. + * + * Here we assume that the elements of the source vector have an + * [] operator where v[0] is the x value and v[1] is the y value + * if the sampling point. + */ + template + void assignFromTupleList_(DestVector &destX, + DestVector &destY, + ListIterator srcBegin, + ListIterator srcEnd, + int numSamples) + { + assert(numSamples >= 2); + + // copy sample points, make sure that the first x value is + // smaller than the last one + + // find out wether the x values are in reverse order + ListIterator it = srcBegin; + ++it; + bool reverse = false; + if (std::get<0>(*srcBegin) > std::get<0>(*it)) + reverse = true; + --it; + + // loop over all sampling points + for (int i = 0; it != srcEnd; ++i, ++it) { + int idx = i; + if (reverse) + idx = numSamples - i - 1; + destX[i] = std::get<0>(*it); + destY[i] = std::get<1>(*it); + } + } + + + /*! + * \brief Make the linear system of equations Mx = d which results + * in the moments of the periodic spline. + * + * When solving Mx = d, it should be noted that x[0] is trash and + * needs to be set to x[n-1] + */ + template + void makePeriodicSystem_(Matrix &M, Vector &d) + { + // a periodic spline only makes sense if the first and the + // last sampling point have the same y value + assert(y_(0) == y_(numSamples_() - 1)); + + makeNaturalSystem(M, d); + + int n = numSamples_() - 1; + + Scalar lambda_n = h_(1)/(h_(n) + h_(1)); + Scalar lambda_1 = M[1][2]; + Scalar mu_1 = 1 - lambda_1; + Scalar mu_n = 1 - lambda_n; + + Scalar d_n = + 6 / (h_(n) + h_(1)) + * + ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n)); + + // insert lambda_n and mu_1 + M[1][n] = mu_1; + M[n][1] = lambda_n; + M[n][n-1] = mu_n; + d[n] = d_n; + } + + /*! + * \brief Make the linear system of equations Mx = d which results + * in the moments of the full spline. + */ + template + void makeFullSystem_(Matrix &M, Vector &d, Scalar m0, Scalar m1) + { + makeNaturalSystem_(M, d); + + int n = numSamples_() - 1; + // first row + M[0][1] = 1; + d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0); + + // last row + M[n][n - 1] = 1; + + // right hand side + d[n] = + 6/h_(n) + * + (m1 - (y_(n) - y_(n - 1))/h_(n)); + } + + /*! + * \brief Make the linear system of equations Mx = d which results + * in the moments of the natural spline. + */ + template + void makeNaturalSystem_(Matrix &M, Vector &d) + { + M = 0; + d = 0; + + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 111 + const int n = numSamples_() - 1; + + // second to next to last rows + for (int i = 1; i < n; ++i) { + const Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i+1)); + const Scalar mu_i = 1 - lambda_i; + const Scalar d_i = + 6 / (h_(i) + h_(i+1)) + * + ( (y_(i+1) - y_(i))/h_(i+1) - (y_(i) - y_(i-1))/h_(i)); + + M[i][i-1] = mu_i; + M[i][i] = 2; + M[i][i + 1] = lambda_i; + d[i] = d_i; + }; + + // first row + M[0][0] = 2; + + // last row + M[n][n] = 2; + // right hand side + d[0] = 0.0; + d[n] = 0.0; + } + + // evaluate the spline at a given the position and given the + // segment index + Scalar eval_(Scalar x, int i) const + { + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 109 + Scalar h_i1 = h_(i + 1); + Scalar x_i = x - x_(i); + Scalar x_i1 = x_(i+1) - x; + + Scalar A_i = + (y_(i+1) - y_(i))/h_i1 + - + h_i1/6*(moment_(i+1) - moment_(i)); + Scalar B_i = y_(i) - moment_(i)* (h_i1*h_i1) / 6; + + return + moment_(i)* x_i1*x_i1*x_i1 / (6 * h_i1) + + + moment_(i + 1)* x_i*x_i*x_i / (6 * h_i1) + + + A_i*x_i + + + B_i; + } + + // evaluate the derivative of a spline given the actual position + // and the segment index + Scalar evalDerivative_(Scalar x, int i) const + { + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 109 + Scalar h_i1 = h_(i + 1); + Scalar x_i = x - x_(i); + Scalar x_i1 = x_(i+1) - x; + + Scalar A_i = + (y_(i+1) - y_(i))/h_i1 + - + h_i1/6*(moment_(i+1) - moment_(i)); + + return + -moment_(i) * x_i1*x_i1 / (2 * h_i1) + + + moment_(i + 1) * x_i*x_i / (2 * h_i1) + + + A_i; + } + + + // returns the monotonicality of an interval of a spline segment + // + // The return value have the following meaning: + // + // 1: spline is monotonously increasing in the specified interval + // 0: spline is not monotonic in the specified interval + // -1: spline is monotonously decreasing in the specified interval + int monotonic_(int i, Scalar x0, Scalar x1) const + { + // shift the interval so that it is consistent with the + // definitions by Stoer + x0 = x0 - x_(i); + x1 = x1 - x_(i); + + Scalar a = a_(i); + Scalar b = b_(i); + Scalar c = c_(i); + + Scalar disc = 4*b*b - 12*a*c; + if (disc < 0) { + // discriminant is smaller than 0, i.e. the segment does + // not exhibit any extrema. + return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1; + } + disc = std::sqrt(disc); + Scalar xE1 = (-2*b + disc)/(6*a); + Scalar xE2 = (-2*b - disc)/(6*a); + + if (disc == 0) { + // saddle point -> no extrema + if (xE1 == x0) + // make sure that we're not picking the saddle point + // to determine whether we're monotonically increasing + // or decreasing + x0 = x1; + return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1; + }; + if ((x0 < xE1 && xE1 < x1) || + (x0 < xE2 && xE2 < x1)) + { + // there is an extremum in the range (x0, x1) + return 0; + } + // no extremum in range (x0, x1) + x0 = (x0 + x1)/2; // pick point in the middle of the interval + // to avoid extrema on the boundaries + return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1; + } + + /*! + * \brief Find all the intersections of a segment of the spline + * with a cubic polynomial within a specified interval. + */ + int intersectSegment_(Scalar *sol, + int segIdx, + Scalar a, Scalar b, Scalar c, Scalar d, + Scalar x0 = -1e100, Scalar x1 = 1e100) const + { + int n = Opm::invertCubicPolynomial(sol, + a_(segIdx) - a, + b_(segIdx) - b, + c_(segIdx) - c, + d_(segIdx) - d); + x0 = std::max(x_(segIdx), x0); + x1 = std::max(x_(segIdx+1), x1); + + // filter the intersections outside of the specified intervall + int k = 0; + for (int j = 0; j < n; ++j) { + sol[j] += x_(segIdx); // add the offset of the intervall. For details see Stoer + if (x0 <= sol[j] && sol[j] <= x1) { + sol[k] = sol[j]; + ++k; + } + } + return k; + } + + // find the segment index for a given x coordinate + int segmentIdx_(Scalar x) const + { + // bisection + int iLow = 0; + int iHigh = numSamples_() - 1; + + while (iLow + 1 < iHigh) { + int i = (iLow + iHigh) / 2; + if (x_(i) > x) + iHigh = i; + else + iLow = i; + }; + return iLow; + } + + /*! + * \brief Returns x[i] - x[i - 1] + */ + Scalar h_(int i) const + { + assert(x_(i) > x_(i-1)); // the sampling points must be given + // in ascending order + return x_(i) - x_(i - 1); + } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar x_(int i) const + { return asImp_().x_(i); } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar y_(int i) const + { return asImp_().y_(i); } + + /*! + * \brief Returns the moment (i.e. second derivative) of the + * spline at the i-th sampling point. + */ + Scalar moment_(int i) const + { return asImp_().moment_(i); } + + // returns the coefficient in front of the x^0 term. In Stoer this + // is delta. + Scalar a_(int i) const + { return (moment_(i+1) - moment_(i))/(6*h_(i+1)); } + + // returns the coefficient in front of the x^2 term In Stoer this + // is gamma. + Scalar b_(int i) const + { return moment_(i)/2; } + + // returns the coefficient in front of the x^1 term. In Stoer this + // is beta. + Scalar c_(int i) const + { return (y_(i+1) - y_(i))/h_(i + 1) - h_(i+1)/6*(2*moment_(i) + moment_(i+1)); } + + // returns the coefficient in front of the x^0 term. In Stoer this + // is alpha. + Scalar d_(int i) const + { return y_(i); } + + /*! + * \brief Returns the number of sampling points. + */ + int numSamples_() const + { return asImp_().numSamples(); } +}; + +} + +#endif diff --git a/opm/common/statictabulated2dfunction.hh b/opm/common/statictabulated2dfunction.hh new file mode 100644 index 000000000..ba57abf12 --- /dev/null +++ b/opm/common/statictabulated2dfunction.hh @@ -0,0 +1,108 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \copydoc Opm::StaticTabulated2DFunction + */ +#ifndef OPM_STATIC_TABULATED_2D_FUNCTION_HH +#define OPM_STATIC_TABULATED_2D_FUNCTION_HH + +#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) { + DUNE_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]; + } +}; +} + +#endif diff --git a/opm/common/tabulated2dfunction.hh b/opm/common/tabulated2dfunction.hh new file mode 100644 index 000000000..48f7d064a --- /dev/null +++ b/opm/common/tabulated2dfunction.hh @@ -0,0 +1,143 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2011-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \copydoc Opm::Tabulated2DFunction + */ +#ifndef OPM_TABULATED_2D_FUNCTION_HH +#define OPM_TABULATED_2D_FUNCTION_HH + +#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)) + { + DUNE_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); } +}; +} + +#endif diff --git a/opm/common/valgrind.hh b/opm/common/valgrind.hh new file mode 100644 index 000000000..990d55334 --- /dev/null +++ b/opm/common/valgrind.hh @@ -0,0 +1,285 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief Some templates to wrap the valgrind macros + */ +#ifndef OPM_VALGRIND_HH +#define OPM_VALGRIND_HH + +#if ! HAVE_VALGRIND && ! defined(DOXYGEN) +namespace Valgrind +{ +bool boolBlubb(bool value) { return value; } +void voidBlubb() { } + +#define SetUndefined(t) voidBlubb() +#define SetDefined(t) voidBlubb() +#define CheckDefined(t) boolBlubb(true) +#define SetNoAccess(t) voidBlubb() +#define IsRunning() boolBlubb(false) +} + +#else + +#include + +namespace Valgrind +{ +/*! + * \ingroup Valgrind + * \brief Returns whether the program is running under Valgrind or not. + */ +inline bool IsRunning() +{ +#if !defined NDEBUG && HAVE_VALGRIND + return RUNNING_ON_VALGRIND; +#else + return false; +#endif +} + + +/*! + * \ingroup Valgrind + * \brief Make valgrind complain if any of the memory occupied by an object + * is undefined. + * + * Please note that this does not check whether the destinations of an + * object's pointers or references are defined. Also, for performance + * reasons the compiler might insert "padding bytes" between within + * the objects which leads to false positives. + * + * Example: + * + * \code + * int i; + * Opm::Valgrind::CheckDefined(i); // Valgrind complains! + * \endcode + * + * \tparam T The type of the object which ought to be checked + * + * \param value the object which valgrind should check + * + * \return true iff there are no undefined bytes in the memory + * occupied by the object. + */ +template +inline bool CheckDefined(const T &value) +{ +#if !defined NDEBUG && HAVE_VALGRIND + unsigned int tmp = VALGRIND_CHECK_MEM_IS_DEFINED(&value, sizeof(T)); + return tmp == 0; +#else + return true; +#endif +} +/*! + * \ingroup Valgrind + * + * * \brief Make valgrind complain if any of the the memory occupied + * by a C-style array objects is undefined. + * + * Please note that this does not check whether the destinations of an + * object's pointers or references are defined. Also, for performance + * reasons the compiler might insert "padding bytes" between within + * the objects which leads to false positives. + * + * Example: + * + * \code + * int i[2]; + * Opm::Valgrind::CheckDefined(i, 2); // Valgrind complains! + * \endcode + * + * \tparam T The type of the object which ought to be checked + * + * \param value Pointer to the first object of the array. + * \param size The size of the array in number of objects + * + * \return true iff there are no undefined bytes in the memory + * occupied by the array. + */ +template +inline bool CheckDefined(const T *value, int size) +{ +#if !defined NDEBUG && HAVE_VALGRIND + unsigned int tmp = VALGRIND_CHECK_MEM_IS_DEFINED(value, size*sizeof(T)); + return tmp == 0; +#else + return true; +#endif +} + +/*! + * \ingroup Valgrind + * \brief Make the memory on which an object resides undefined in + * valgrind runs. + * + * Example: + * + * \code + * int i = 0; + * Opm::Valgrind::SetUndefined(i); + * Opm::Valgrind::CheckDefined(i); // Valgrind complains! + * \endcode + * + * \tparam T The type of the object which ought to be set to undefined + * + * \param value The object which's memory valgrind should be told is undefined + */ +template +inline void SetUndefined(const T &value) +{ +#if !defined NDEBUG && HAVE_VALGRIND + auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(&value, sizeof(T)); +#endif +} + +/*! + * \ingroup Valgrind + * \brief Make the memory on which an array of object resides + * undefined in valgrind runs. + * + * Example: + * + * \code + * int i[3] = {0, 1, 3}; + * Opm::Valgrind::SetUndefined(&i[1], 2); + * Opm::Valgrind::CheckDefined(i, 3); // Valgrind complains! + * \endcode + * + * \tparam T The type of the object which ought to be set to undefined + * + * \param value Pointer to the first object of the array. + * \param size The size of the array in number of objects + */ +template +inline void SetUndefined(const T *value, int size) +{ +#if !defined NDEBUG && HAVE_VALGRIND + auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(value, size*sizeof(T)); +#endif +} + +/*! + * \ingroup Valgrind + * \brief Make the memory on which an object resides defined. + * + * Example: + * + * \code + * int i; + * Opm::Valgrind::SetDefined(i); + * Opm::Valgrind::CheckDefined(i); // Valgrind does not complain! + * \endcode + * + * \tparam T The type of the object which valgrind should consider as defined + * + * \param value The object which's memory valgrind should consider as defined + */ +template +inline void SetDefined(const T &value) +{ +#if !defined NDEBUG && HAVE_VALGRIND + auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(&value, sizeof(T)); +#endif +} + +/*! + * \ingroup Valgrind + * \brief Make the memory on which a C-style array of objects resides + * defined. + * + * Example: + * + * \code + * int i[3]; + * Opm::Valgrind::SetDefined(i, 3); + * Opm::Valgrind::CheckDefined(i, 3); // Valgrind does not complain! + * \endcode + * + * \tparam T The type of the object which valgrind should consider as defined + * + * \param value Pointer to the first object of the array. + * \param n The size of the array in number of objects + */ +template +inline void SetDefined(const T *value, int n) +{ +#if !defined NDEBUG && HAVE_VALGRIND + auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(value, n*sizeof(T)); +#endif +} + +/*! + * \ingroup Valgrind + * \brief Make valgrind complain if an object's memory is accessed. + * + * Example: + * + * \code + * int i = 1; + * Opm::Valgrind::SetNoAccess(i); + * int j = i; // Valgrind complains! + * \endcode + * + * \tparam T The type of the object which valgrind should complain if accessed + * + * \param value The object which's memory valgrind should complain if accessed + */ +template +inline void SetNoAccess(const T &value) +{ +#if !defined NDEBUG && HAVE_VALGRIND + auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(&value, sizeof(T)); +#endif +} + +/*! + * \ingroup Valgrind + * \brief Make valgrind complain if the memory of a C-style array of + * objects is accessed. + * + * Example: + * + * \code + * int i[3] = {0, 1, 2}; + * Opm::Valgrind::SetNoAccess(i, 2); + * int j = i[1]; // Valgrind complains! + * \endcode + * + * \param value Pointer to the first object of the array. + * \param size The size of the array in number of objects + * + * \param value The object which's memory valgrind should complain if accessed + */ +template +inline void SetNoAccess(const T *value, int size) +{ +#if !defined NDEBUG && HAVE_VALGRIND + auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(value, size*sizeof(T)); +#endif +} + +} + +#endif + +#endif diff --git a/opm/common/variablelengthspline_.hh b/opm/common/variablelengthspline_.hh new file mode 100644 index 000000000..6e743e253 --- /dev/null +++ b/opm/common/variablelengthspline_.hh @@ -0,0 +1,479 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc Opm::VariableLengthSpline_ + */ +#ifndef OPM_VARIABLE_LENGTH_SPLINE_HH +#define OPM_VARIABLE_LENGTH_SPLINE_HH + +#include +#include +#include +#include + +#include "splinecommon_.hh" + +namespace Opm +{ +/*! + * \brief The common code for all 3rd order polynomial splines with + * where the number of sampling points only known at run-time. + */ +template +class VariableLengthSpline_ + : public SplineCommon_ > +{ + friend class SplineCommon_ >; + + typedef ScalarT Scalar; + typedef Dune::BlockVector > Vector; + typedef Dune::BTDMatrix > BTDMatrix; + +public: + /*! + * \brief Returns the number of sampling points. + */ + int numSamples() const + { return xPos_.size(); } + + + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + // Full splines // + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using C-style arrays. + * + * This method uses separate array-like objects for the values of + * the X and Y coordinates. In this context 'array-like' means + * that an access to the members is provided via the [] + * operator. (e.g. C arrays, std::vector, std::array, etc.) Each + * array must be of size 'nSamples' at least. Also, the number of + * sampling points must be larger than 1. + */ + template + void setXYArrays(int nSamples, + const ScalarArrayX &x, + const ScalarArrayY &y, + Scalar m0, Scalar m1) + { + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = x[i]; + yPos_[i] = y[i]; + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using STL-compatible containers. + * + * This method uses separate STL-compatible containers for the + * values of the sampling points' X and Y + * coordinates. "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y, + Scalar m0, Scalar m1) + { + assert(x.size() == y.size()); + assert(x.size() > 1); + + setNumSamples_(x.size()); + + std::copy(x.begin(), x.end(), xPos_.begin()); + std::copy(y.begin(), y.end(), yPos_.begin()); + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a C-style array. + * + * This method uses a single array of sampling points, which are + * seen as an array-like object which provides access to the X and + * Y coordinates. In this context 'array-like' means that an + * access to the members is provided via the [] operator. (e.g. C + * arrays, std::vector, std::array, etc.) The array containing + * the sampling points must be of size 'nSamples' at least. Also, + * the number of sampling points must be larger than 1. + */ + template + void setArrayOfPoints(int nSamples, + const PointArray &points, + Scalar m0, + Scalar m1) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = points[i][0]; + yPos_[i] = points[i][1]; + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a STL-compatible container of + * array-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be array-like objects storing the + * X and Y coordinates. "STL-compatible" means that the container + * provides access to iterators using the begin(), end() methods + * and also provides a size() method. Also, the number of entries + * in the X and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfPoints(const XYContainer &points, + Scalar m0, + Scalar m1) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(points.size() > 1); + + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = (*it)[0]; + yPos_[i] = (*it)[1]; + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + // make a full spline + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a STL-compatible container of + * tuple-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be tuple-like objects storing the + * X and Y coordinates. "tuple-like" means that the objects + * provide access to the x values via std::get<0>(obj) and to the + * y value via std::get<1>(obj) (e.g. std::tuple or + * std::pair). "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfTuples(const XYContainer &points, + Scalar m0, + Scalar m1) + { + // resize internal arrays + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = std::get<0>(*it); + yPos_[i] = std::get<1>(*it); + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + // make a full spline + makeFullSpline_(m0, m1); + } + + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + // Natural splines // + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + /*! + * \brief Set the sampling points natural spline using C-style arrays. + * + * This method uses separate array-like objects for the values of + * the X and Y coordinates. In this context 'array-like' means + * that an access to the members is provided via the [] + * operator. (e.g. C arrays, std::vector, std::array, etc.) Each + * array must be of size 'nSamples' at least. Also, the number of + * sampling points must be larger than 1. + */ + template + void setXYArrays(int nSamples, + const ScalarArrayX &x, + const ScalarArrayY &y) + { + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = x[i]; + yPos_[i] = y[i]; + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using + * STL-compatible containers. + * + * This method uses separate STL-compatible containers for the + * values of the sampling points' X and Y + * coordinates. "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y) + { + assert(x.size() == y.size()); + assert(x.size() > 1); + + setNumSamples_(x.size()); + std::copy(x.begin(), x.end(), xPos_.begin()); + std::copy(y.begin(), y.end(), yPos_.begin()); + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * C-style array. + * + * This method uses a single array of sampling points, which are + * seen as an array-like object which provides access to the X and + * Y coordinates. In this context 'array-like' means that an + * access to the members is provided via the [] operator. (e.g. C + * arrays, std::vector, std::array, etc.) The array containing + * the sampling points must be of size 'nSamples' at least. Also, + * the number of sampling points must be larger than 1. + */ + template + void setArrayOfPoints(int nSamples, + const PointArray &points) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = points[i][0]; + yPos_[i] = points[i][1]; + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * STL-compatible container of array-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be array-like objects storing the + * X and Y coordinates. "STL-compatible" means that the container + * provides access to iterators using the begin(), end() methods + * and also provides a size() method. Also, the number of entries + * in the X and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfPoints(const XYContainer &points) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(points.size() > 1); + + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++ i, ++it) { + xPos_[i] = (*it)[0]; + yPos_[i] = (*it)[1]; + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + // make a natural spline + makeNaturalSpline_(); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * STL-compatible container of tuple-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be tuple-like objects storing the + * X and Y coordinates. "tuple-like" means that the objects + * provide access to the x values via std::get<0>(obj) and to the + * y value via std::get<1>(obj) (e.g. std::tuple or + * std::pair). "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfTuples(const XYContainer &points) + { + // resize internal arrays + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = std::get<0>(*it); + yPos_[i] = std::get<1>(*it); + } + + if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + // make a natural spline + makeNaturalSpline_(); + } + +protected: + /*! + * \brief Reverse order of the elements in the arrays which + * contain the sampling points. + */ + void reverseSamplingPoints_() + { + // reverse the arrays + int n = numSamples(); + for (int i = 0; i <= (n - 1)/2; ++i) { + std::swap(xPos_[i], xPos_[n - i - 1]); + std::swap(yPos_[i], yPos_[n - i - 1]); + } + } + + /*! + * \brief Resizes the internal vectors to store the sample points. + */ + void setNumSamples_(int nSamples) + { + xPos_.resize(nSamples); + yPos_.resize(nSamples); + m_.resize(nSamples); + } + + /*! + * \brief Create a natural spline from the already set sampling points. + * + * Also creates temporary matrix and right hand side vector. + */ + void makeFullSpline_(Scalar m0, Scalar m1) + { + BTDMatrix M(numSamples()); + Vector d(numSamples()); + + // create linear system of equations + this->makeFullSystem_(M, d, m0, m1); + + // solve for the moments + M.solve(m_, d); + } + + /*! + * \brief Create a natural spline from the already set sampling points. + * + * Also creates temporary matrix and right hand side vector. + */ + void makeNaturalSpline_() + { + BTDMatrix M(numSamples()); + Vector d(numSamples()); + + // create linear system of equations + this->makeNaturalSystem_(M, d); + + // solve for the moments + M.solve(m_, d); + } + + /*! + * \brief Returns the x coordinate of the i-th sampling point. + */ + Scalar x_(int i) const + { return xPos_[i]; } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar y_(int i) const + { return yPos_[i]; } + + /*! + * \brief Returns the moment (i.e. second derivative) of the + * spline at the i-th sampling point. + */ + Scalar moment_(int i) const + { return m_[i]; } + + Vector xPos_; + Vector yPos_; + Vector m_; +}; +} + +#endif diff --git a/opm/material/binarycoefficients/air_mesitylene.hh b/opm/material/binarycoefficients/air_mesitylene.hh index 1f7023592..db9c166a9 100644 --- a/opm/material/binarycoefficients/air_mesitylene.hh +++ b/opm/material/binarycoefficients/air_mesitylene.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BinaryCoeff::Air_Mesitylene */ -#ifndef EWOMS_BINARY_COEFF_AIR_MESITYLENE_HH -#define EWOMS_BINARY_COEFF_AIR_MESITYLENE_HH +#ifndef OPM_BINARY_COEFF_AIR_MESITYLENE_HH +#define OPM_BINARY_COEFF_AIR_MESITYLENE_HH #include #include diff --git a/opm/material/binarycoefficients/air_xylene.hh b/opm/material/binarycoefficients/air_xylene.hh index 19f693a6a..82a19c7e4 100644 --- a/opm/material/binarycoefficients/air_xylene.hh +++ b/opm/material/binarycoefficients/air_xylene.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BinaryCoeff::Air_Xylene */ -#ifndef EWOMS_BINARY_COEFF_AIR_XYLENE_HH -#define EWOMS_BINARY_COEFF_AIR_XYLENE_HH +#ifndef OPM_BINARY_COEFF_AIR_XYLENE_HH +#define OPM_BINARY_COEFF_AIR_XYLENE_HH #include #include diff --git a/opm/material/binarycoefficients/brine_co2.hh b/opm/material/binarycoefficients/brine_co2.hh index ff430bf4d..8460806e1 100644 --- a/opm/material/binarycoefficients/brine_co2.hh +++ b/opm/material/binarycoefficients/brine_co2.hh @@ -21,8 +21,8 @@ * * \copydoc Opm::BinaryCoeff::Brine_CO2 */ -#ifndef EWOMS_BINARY_COEFF_BRINE_CO2_HH -#define EWOMS_BINARY_COEFF_BRINE_CO2_HH +#ifndef OPM_BINARY_COEFF_BRINE_CO2_HH +#define OPM_BINARY_COEFF_BRINE_CO2_HH #include #include diff --git a/opm/material/binarycoefficients/fullermethod.hh b/opm/material/binarycoefficients/fullermethod.hh index 24a1c3122..958fb1b6f 100644 --- a/opm/material/binarycoefficients/fullermethod.hh +++ b/opm/material/binarycoefficients/fullermethod.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::BinaryCoeff::fullerMethod */ -#ifndef EWOMS_FULLERMETHOD_HH -#define EWOMS_FULLERMETHOD_HH +#ifndef OPM_FULLERMETHOD_HH +#define OPM_FULLERMETHOD_HH #include diff --git a/opm/material/binarycoefficients/h2o_air.hh b/opm/material/binarycoefficients/h2o_air.hh index 6be1796f9..e7d81126e 100644 --- a/opm/material/binarycoefficients/h2o_air.hh +++ b/opm/material/binarycoefficients/h2o_air.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BinaryCoeff::H2O_Air */ -#ifndef EWOMS_BINARY_COEFF_H2O_AIR_HH -#define EWOMS_BINARY_COEFF_H2O_AIR_HH +#ifndef OPM_BINARY_COEFF_H2O_AIR_HH +#define OPM_BINARY_COEFF_H2O_AIR_HH #include diff --git a/opm/material/binarycoefficients/h2o_co2.hh b/opm/material/binarycoefficients/h2o_co2.hh index b8ce03ab7..fc0bde485 100644 --- a/opm/material/binarycoefficients/h2o_co2.hh +++ b/opm/material/binarycoefficients/h2o_co2.hh @@ -21,8 +21,8 @@ * * \copydoc Opm::BinaryCoeff::H2O_CO2 */ -#ifndef EWOMS_BINARY_COEFF_H2O_CO2_HH -#define EWOMS_BINARY_COEFF_H2O_CO2_HH +#ifndef OPM_BINARY_COEFF_H2O_CO2_HH +#define OPM_BINARY_COEFF_H2O_CO2_HH #include #include diff --git a/opm/material/binarycoefficients/h2o_mesitylene.hh b/opm/material/binarycoefficients/h2o_mesitylene.hh index bd50093f4..51fff1a6e 100644 --- a/opm/material/binarycoefficients/h2o_mesitylene.hh +++ b/opm/material/binarycoefficients/h2o_mesitylene.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BinaryCoeff::H2O_Mesitylene */ -#ifndef EWOMS_BINARY_COEFF_H2O_MESITYLENE_HH -#define EWOMS_BINARY_COEFF_H2O_MESITYLENE_HH +#ifndef OPM_BINARY_COEFF_H2O_MESITYLENE_HH +#define OPM_BINARY_COEFF_H2O_MESITYLENE_HH #include #include @@ -51,9 +51,9 @@ public: // after Sanders Scalar sanderH = 1.7e-1; // [M/atm] //conversion to our Henry definition - Scalar ewomsH = sanderH / 101.325; // has now [(mol/m^3)/Pa] - ewomsH *= 18.02e-6; // multiplied by molar volume of reference phase = water - return 1.0/ewomsH; // [Pa] + Scalar opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa] + opmH *= 18.02e-6; // multiplied by molar volume of reference phase = water + return 1.0/opmH; // [Pa] } /*! diff --git a/opm/material/binarycoefficients/h2o_n2.hh b/opm/material/binarycoefficients/h2o_n2.hh index 74a2ccfab..635d2bc77 100644 --- a/opm/material/binarycoefficients/h2o_n2.hh +++ b/opm/material/binarycoefficients/h2o_n2.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BinaryCoeff::H2O_N2 */ -#ifndef EWOMS_BINARY_COEFF_H2O_N2_HH -#define EWOMS_BINARY_COEFF_H2O_N2_HH +#ifndef OPM_BINARY_COEFF_H2O_N2_HH +#define OPM_BINARY_COEFF_H2O_N2_HH #include "henryiapws.hh" #include "fullermethod.hh" diff --git a/opm/material/binarycoefficients/h2o_xylene.hh b/opm/material/binarycoefficients/h2o_xylene.hh index 781626e62..07430198b 100644 --- a/opm/material/binarycoefficients/h2o_xylene.hh +++ b/opm/material/binarycoefficients/h2o_xylene.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BinaryCoeff::H2O_Xylene */ -#ifndef EWOMS_BINARY_COEFF_H2O_XYLENE_HH -#define EWOMS_BINARY_COEFF_H2O_XYLENE_HH +#ifndef OPM_BINARY_COEFF_H2O_XYLENE_HH +#define OPM_BINARY_COEFF_H2O_XYLENE_HH #include #include @@ -52,9 +52,9 @@ public: // after Sanders Scalar sanderH = 1.5e-1; //[M/atm] //conversion to our Henry definition - Scalar ewomsH = sanderH / 101.325; // has now [(mol/m^3)/Pa] - ewomsH *= 18.02e-6; //multiplied by molar volume of reference phase = water - return 1.0/ewomsH; // [Pa] + Scalar opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa] + opmH *= 18.02e-6; //multiplied by molar volume of reference phase = water + return 1.0/opmH; // [Pa] } /*! diff --git a/opm/material/binarycoefficients/henryiapws.hh b/opm/material/binarycoefficients/henryiapws.hh index b1651266c..bea46ba30 100644 --- a/opm/material/binarycoefficients/henryiapws.hh +++ b/opm/material/binarycoefficients/henryiapws.hh @@ -21,8 +21,8 @@ * * \brief The IAPWS formulation of Henry coefficients in water. */ -#ifndef EWOMS_HENRY_IAPWS_HH -#define EWOMS_HENRY_IAPWS_HH +#ifndef OPM_HENRY_IAPWS_HH +#define OPM_HENRY_IAPWS_HH #include @@ -81,4 +81,4 @@ inline Scalar henryIAPWS(Scalar E, } } -#endif // EWOMS_HENRY_IAPWS_HH +#endif // OPM_HENRY_IAPWS_HH diff --git a/opm/material/components/air.hh b/opm/material/components/air.hh index 27485c0c2..c97e6f44d 100644 --- a/opm/material/components/air.hh +++ b/opm/material/components/air.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::Air */ -#ifndef EWOMS_AIR_HH -#define EWOMS_AIR_HH +#ifndef OPM_AIR_HH +#define OPM_AIR_HH #include #include diff --git a/opm/material/components/brine.hh b/opm/material/components/brine.hh index bc2eb9450..61ee3d39a 100644 --- a/opm/material/components/brine.hh +++ b/opm/material/components/brine.hh @@ -21,8 +21,8 @@ * * \copydoc Opm::Brine */ -#ifndef EWOMS_BRINE_HH -#define EWOMS_BRINE_HH +#ifndef OPM_BRINE_HH +#define OPM_BRINE_HH #include diff --git a/opm/material/components/co2.hh b/opm/material/components/co2.hh index 1f5e63f7e..c6f66bf7c 100644 --- a/opm/material/components/co2.hh +++ b/opm/material/components/co2.hh @@ -21,8 +21,8 @@ * * \copydoc Opm::CO2 */ -#ifndef EWOMS_CO2_HH -#define EWOMS_CO2_HH +#ifndef OPM_CO2_HH +#define OPM_CO2_HH #include #include diff --git a/opm/material/components/component.hh b/opm/material/components/component.hh index c892e2432..cc64d274b 100644 --- a/opm/material/components/component.hh +++ b/opm/material/components/component.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::Component */ -#ifndef EWOMS_COMPONENT_HH -#define EWOMS_COMPONENT_HH +#ifndef OPM_COMPONENT_HH +#define OPM_COMPONENT_HH #include diff --git a/opm/material/components/dnapl.hh b/opm/material/components/dnapl.hh index 513cd1460..781443261 100644 --- a/opm/material/components/dnapl.hh +++ b/opm/material/components/dnapl.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::DNAPL */ -#ifndef EWOMS_DNAPL_HH -#define EWOMS_DNAPL_HH +#ifndef OPM_DNAPL_HH +#define OPM_DNAPL_HH #include #include "component.hh" diff --git a/opm/material/components/h2o.hh b/opm/material/components/h2o.hh index 862d3b40d..56323b557 100644 --- a/opm/material/components/h2o.hh +++ b/opm/material/components/h2o.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::H2O */ -#ifndef EWOMS_H2O_HH -#define EWOMS_H2O_HH +#ifndef OPM_H2O_HH +#define OPM_H2O_HH #include #include diff --git a/opm/material/components/iapws/common.hh b/opm/material/components/iapws/common.hh index 9fd65e55c..10a87171f 100644 --- a/opm/material/components/iapws/common.hh +++ b/opm/material/components/iapws/common.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::IAPWS::Common */ -#ifndef EWOMS_IAPWS_COMMON_HH -#define EWOMS_IAPWS_COMMON_HH +#ifndef OPM_IAPWS_COMMON_HH +#define OPM_IAPWS_COMMON_HH #include diff --git a/opm/material/components/iapws/region1.hh b/opm/material/components/iapws/region1.hh index 6b9c18eb2..4f765b818 100644 --- a/opm/material/components/iapws/region1.hh +++ b/opm/material/components/iapws/region1.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::IAPWS::Region1 */ -#ifndef EWOMS_IAPWS_REGION1_HH -#define EWOMS_IAPWS_REGION1_HH +#ifndef OPM_IAPWS_REGION1_HH +#define OPM_IAPWS_REGION1_HH #include diff --git a/opm/material/components/iapws/region2.hh b/opm/material/components/iapws/region2.hh index 938b30c33..a21c4c704 100644 --- a/opm/material/components/iapws/region2.hh +++ b/opm/material/components/iapws/region2.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::IAPWS::Region2 */ -#ifndef EWOMS_IAPWS_REGION2_HH -#define EWOMS_IAPWS_REGION2_HH +#ifndef OPM_IAPWS_REGION2_HH +#define OPM_IAPWS_REGION2_HH #include diff --git a/opm/material/components/iapws/region4.hh b/opm/material/components/iapws/region4.hh index 0c85f82ee..835065b3b 100644 --- a/opm/material/components/iapws/region4.hh +++ b/opm/material/components/iapws/region4.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::IAPWS::Region4 */ -#ifndef EWOMS_IAPWS_REGION4_HH -#define EWOMS_IAPWS_REGION4_HH +#ifndef OPM_IAPWS_REGION4_HH +#define OPM_IAPWS_REGION4_HH #include diff --git a/opm/material/components/lnapl.hh b/opm/material/components/lnapl.hh index 5b91c41b2..3a3e9a90d 100644 --- a/opm/material/components/lnapl.hh +++ b/opm/material/components/lnapl.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::LNAPL */ -#ifndef EWOMS_LNAPL_HH -#define EWOMS_LNAPL_HH +#ifndef OPM_LNAPL_HH +#define OPM_LNAPL_HH #include "component.hh" diff --git a/opm/material/components/mesitylene.hh b/opm/material/components/mesitylene.hh index 4d5e15e4d..b23e16cbb 100644 --- a/opm/material/components/mesitylene.hh +++ b/opm/material/components/mesitylene.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::Mesitylene */ -#ifndef EWOMS_MESITYLENE_HH -#define EWOMS_MESITYLENE_HH +#ifndef OPM_MESITYLENE_HH +#define OPM_MESITYLENE_HH #include #include diff --git a/opm/material/components/n2.hh b/opm/material/components/n2.hh index 101529a2a..e207a82e2 100644 --- a/opm/material/components/n2.hh +++ b/opm/material/components/n2.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::N2 */ -#ifndef EWOMS_N2_HH -#define EWOMS_N2_HH +#ifndef OPM_N2_HH +#define OPM_N2_HH #include diff --git a/opm/material/components/nullcomponent.hh b/opm/material/components/nullcomponent.hh index 551b7ab5d..519279d36 100644 --- a/opm/material/components/nullcomponent.hh +++ b/opm/material/components/nullcomponent.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::NullComponent */ -#ifndef EWOMS_NULL_COMPONENT_HH -#define EWOMS_NULL_COMPONENT_HH +#ifndef OPM_NULL_COMPONENT_HH +#define OPM_NULL_COMPONENT_HH #include "component.hh" diff --git a/opm/material/components/simpleco2.hh b/opm/material/components/simpleco2.hh index 21a4817fc..db97be569 100644 --- a/opm/material/components/simpleco2.hh +++ b/opm/material/components/simpleco2.hh @@ -23,8 +23,8 @@ * * \copydoc Opm::SimpleCO2 */ -#ifndef EWOMS_SIMPLE_CO2_HH -#define EWOMS_SIMPLE_CO2_HH +#ifndef OPM_SIMPLE_CO2_HH +#define OPM_SIMPLE_CO2_HH #include #include diff --git a/opm/material/components/simpleh2o.hh b/opm/material/components/simpleh2o.hh index bd5c2bb19..8a016c58d 100644 --- a/opm/material/components/simpleh2o.hh +++ b/opm/material/components/simpleh2o.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::SimpleH2O */ -#ifndef EWOMS_SIMPLE_H2O_HH -#define EWOMS_SIMPLE_H2O_HH +#ifndef OPM_SIMPLE_H2O_HH +#define OPM_SIMPLE_H2O_HH #include diff --git a/opm/material/components/tabulatedcomponent.hh b/opm/material/components/tabulatedcomponent.hh index 4d2945355..32beb9186 100644 --- a/opm/material/components/tabulatedcomponent.hh +++ b/opm/material/components/tabulatedcomponent.hh @@ -23,8 +23,8 @@ * * \copydoc Opm::TabulatedComponent */ -#ifndef EWOMS_TABULATED_COMPONENT_HH -#define EWOMS_TABULATED_COMPONENT_HH +#ifndef OPM_TABULATED_COMPONENT_HH +#define OPM_TABULATED_COMPONENT_HH #include #include diff --git a/opm/material/components/unit.hh b/opm/material/components/unit.hh index 115f940da..7f4f1dcdd 100644 --- a/opm/material/components/unit.hh +++ b/opm/material/components/unit.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::Unit */ -#ifndef EWOMS_UNIT_HH -#define EWOMS_UNIT_HH +#ifndef OPM_UNIT_HH +#define OPM_UNIT_HH #include "component.hh" diff --git a/opm/material/components/xylene.hh b/opm/material/components/xylene.hh index 16de79cd1..ef782ec38 100644 --- a/opm/material/components/xylene.hh +++ b/opm/material/components/xylene.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::Xylene */ -#ifndef EWOMS_XYLENE_HH -#define EWOMS_XYLENE_HH +#ifndef OPM_XYLENE_HH +#define OPM_XYLENE_HH #include #include diff --git a/opm/material/constants.hh b/opm/material/constants.hh index 41a640620..c69c0759e 100644 --- a/opm/material/constants.hh +++ b/opm/material/constants.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::Constants */ -#ifndef EWOMS_CONSTANTS_HH -#define EWOMS_CONSTANTS_HH +#ifndef OPM_CONSTANTS_HH +#define OPM_CONSTANTS_HH #include diff --git a/opm/material/constraintsolvers/compositionfromfugacities.hh b/opm/material/constraintsolvers/compositionfromfugacities.hh index d61071151..685169ad7 100644 --- a/opm/material/constraintsolvers/compositionfromfugacities.hh +++ b/opm/material/constraintsolvers/compositionfromfugacities.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::CompositionFromFugacities */ -#ifndef EWOMS_COMPOSITION_FROM_FUGACITIES_HH -#define EWOMS_COMPOSITION_FROM_FUGACITIES_HH +#ifndef OPM_COMPOSITION_FROM_FUGACITIES_HH +#define OPM_COMPOSITION_FROM_FUGACITIES_HH #include #include diff --git a/opm/material/constraintsolvers/computefromreferencephase.hh b/opm/material/constraintsolvers/computefromreferencephase.hh new file mode 100644 index 000000000..8bb312902 --- /dev/null +++ b/opm/material/constraintsolvers/computefromreferencephase.hh @@ -0,0 +1,169 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2011-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc Opm::ComputeFromReferencePhase + */ +#ifndef OPM_COMPUTE_FROM_REFERENCE_PHASE_HH +#define OPM_COMPUTE_FROM_REFERENCE_PHASE_HH + +#include + +#include +#include + +#include + +namespace Opm { + +/*! + * \brief Computes all quantities of a generic fluid state if a + * reference phase has been specified. + * + * This makes it is possible to specify just one phase and let the + * remaining ones be calculated by the constraint solver. This + * constraint solver assumes thermodynamic equilibrium. It assumes the + * following quantities to be set: + * + * - composition (mole+mass fractions) of the *reference* phase + * - temperature of the *reference* phase + * - saturations of *all* phases + * - pressures of *all* phases + * + * after calling the solve() method the following quantities are + * calculated in addition: + * + * - temperature of *all* phases + * - density, molar density, molar volume of *all* phases + * - composition in mole and mass fractions and molaries of *all* phases + * - mean molar masses of *all* phases + * - fugacity coefficients of *all* components in *all* phases + * - if the setViscosity parameter is true, also dynamic viscosities of *all* phases + * - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *all* phases + */ +template +class ComputeFromReferencePhase +{ + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + typedef Opm::CompositionFromFugacities CompositionFromFugacities; + typedef Dune::FieldVector ComponentVector; + +public: + /*! + * \brief Computes all quantities of a generic fluid state if a + * reference phase has been specified. + * + * This makes it is possible to specify just one phase and let the + * remaining ones be calculated by the constraint solver. This + * constraint solver assumes thermodynamic equilibrium. It assumes the + * following quantities to be set: + * + * - composition (mole+mass fractions) of the *reference* phase + * - temperature of the *all* phases + * - saturations of *all* phases + * - pressures of *all* phases + * + * after calling the solve() method the following quantities are + * calculated in addition: + * + * - temperature of *all* phases + * - density, molar density, molar volume of *all* phases + * - composition in mole and mass fractions and molaries of *all* phases + * - mean molar masses of *all* phases + * - fugacity coefficients of *all* components in *all* phases + * - if the setViscosity parameter is true, also dynamic viscosities of *all* phases + * - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *all* phases + * + * \param fluidState Thermodynamic state of the fluids + * \param paramCache Container for cache parameters + * \param refPhaseIdx The phase index of the reference phase + * \param setViscosity Specify whether the dynamic viscosity of + * each phase should also be set. + * \param setEnthalpy Specify whether the specific + * enthalpy/internal energy of each phase + * should also be set. + */ + template + static void solve(FluidState &fluidState, + ParameterCache ¶mCache, + int refPhaseIdx, + bool setViscosity, + bool setEnthalpy) + { + + // compute the density and enthalpy of the + // reference phase + paramCache.updatePhase(fluidState, refPhaseIdx); + fluidState.setDensity(refPhaseIdx, + FluidSystem::density(fluidState, + paramCache, + refPhaseIdx)); + + if (setEnthalpy) + fluidState.setEnthalpy(refPhaseIdx, + FluidSystem::enthalpy(fluidState, + paramCache, + refPhaseIdx)); + + if (setViscosity) + fluidState.setViscosity(refPhaseIdx, + FluidSystem::viscosity(fluidState, + paramCache, + refPhaseIdx)); + + // compute the fugacities of all components in the reference phase + for (int compIdx = 0; compIdx < numComponents; ++compIdx) { + fluidState.setFugacityCoefficient(refPhaseIdx, + compIdx, + FluidSystem::fugacityCoefficient(fluidState, + paramCache, + refPhaseIdx, + compIdx)); + } + + // compute all quantities for the non-reference phases + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (phaseIdx == refPhaseIdx) + continue; // reference phase is already calculated + + ComponentVector fugVec; + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + fugVec[compIdx] = fluidState.fugacity(refPhaseIdx, compIdx); + + CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fugVec); + + if (setViscosity) + fluidState.setViscosity(phaseIdx, + FluidSystem::viscosity(fluidState, + paramCache, + phaseIdx)); + + if (setEnthalpy) + fluidState.setEnthalpy(phaseIdx, + FluidSystem::enthalpy(fluidState, + paramCache, + phaseIdx)); + } + } +}; + +} // end namespace Opm + +#endif diff --git a/opm/material/constraintsolvers/immiscibleflash.hh b/opm/material/constraintsolvers/immiscibleflash.hh index dd4548f92..c5df21a52 100644 --- a/opm/material/constraintsolvers/immiscibleflash.hh +++ b/opm/material/constraintsolvers/immiscibleflash.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::ImmiscibleFlash */ -#ifndef EWOMS_IMMISCIBLE_FLASH_HH -#define EWOMS_IMMISCIBLE_FLASH_HH +#ifndef OPM_IMMISCIBLE_FLASH_HH +#define OPM_IMMISCIBLE_FLASH_HH #include #include diff --git a/opm/material/constraintsolvers/misciblemultiphasecomposition.hh b/opm/material/constraintsolvers/misciblemultiphasecomposition.hh index 47df8f67d..3524cd88e 100644 --- a/opm/material/constraintsolvers/misciblemultiphasecomposition.hh +++ b/opm/material/constraintsolvers/misciblemultiphasecomposition.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::MiscibleMultiPhaseComposition */ -#ifndef EWOMS_MISCIBLE_MULTIPHASE_COMPOSITION_HH -#define EWOMS_MISCIBLE_MULTIPHASE_COMPOSITION_HH +#ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HH +#define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HH #include #include diff --git a/opm/material/constraintsolvers/ncpflash.hh b/opm/material/constraintsolvers/ncpflash.hh index 52380b67b..059fdc590 100644 --- a/opm/material/constraintsolvers/ncpflash.hh +++ b/opm/material/constraintsolvers/ncpflash.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::NcpFlash */ -#ifndef EWOMS_NCP_FLASH_HH -#define EWOMS_NCP_FLASH_HH +#ifndef OPM_NCP_FLASH_HH +#define OPM_NCP_FLASH_HH #include #include diff --git a/opm/material/eos/pengrobinson.hh b/opm/material/eos/pengrobinson.hh index 3324ba2fd..1b58f27af 100644 --- a/opm/material/eos/pengrobinson.hh +++ b/opm/material/eos/pengrobinson.hh @@ -21,8 +21,8 @@ * * \copydoc Opm::PengRobinson */ -#ifndef EWOMS_PENG_ROBINSON_HH -#define EWOMS_PENG_ROBINSON_HH +#ifndef OPM_PENG_ROBINSON_HH +#define OPM_PENG_ROBINSON_HH #include #include @@ -30,6 +30,8 @@ #include #include +#include + #include namespace Opm { diff --git a/opm/material/eos/pengrobinsonmixture.hh b/opm/material/eos/pengrobinsonmixture.hh index dff79be69..c6d4e65e7 100644 --- a/opm/material/eos/pengrobinsonmixture.hh +++ b/opm/material/eos/pengrobinsonmixture.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::PengRobinsonMixture */ -#ifndef EWOMS_PENG_ROBINSON_MIXTURE_HH -#define EWOMS_PENG_ROBINSON_MIXTURE_HH +#ifndef OPM_PENG_ROBINSON_MIXTURE_HH +#define OPM_PENG_ROBINSON_MIXTURE_HH #include "pengrobinson.hh" diff --git a/opm/material/eos/pengrobinsonparams.hh b/opm/material/eos/pengrobinsonparams.hh index 09f8a3361..f6d130cbd 100644 --- a/opm/material/eos/pengrobinsonparams.hh +++ b/opm/material/eos/pengrobinsonparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::PengRobinsonParams */ -#ifndef EWOMS_PENG_ROBINSON_PARAMS_HH -#define EWOMS_PENG_ROBINSON_PARAMS_HH +#ifndef OPM_PENG_ROBINSON_PARAMS_HH +#define OPM_PENG_ROBINSON_PARAMS_HH #include diff --git a/opm/material/eos/pengrobinsonparamsmixture.hh b/opm/material/eos/pengrobinsonparamsmixture.hh index 64d5f8b01..bc4abacc4 100644 --- a/opm/material/eos/pengrobinsonparamsmixture.hh +++ b/opm/material/eos/pengrobinsonparamsmixture.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::PengRobinsonParamsMixture */ -#ifndef EWOMS_PENG_ROBINSON_PARAMS_MIXTURE_HH -#define EWOMS_PENG_ROBINSON_PARAMS_MIXTURE_HH +#ifndef OPM_PENG_ROBINSON_PARAMS_MIXTURE_HH +#define OPM_PENG_ROBINSON_PARAMS_MIXTURE_HH #include #include diff --git a/opm/material/fluidmatrixinteractions/2p/brookscorey.hh b/opm/material/fluidmatrixinteractions/2p/brookscorey.hh index e5b411b91..ee7b3a56e 100644 --- a/opm/material/fluidmatrixinteractions/2p/brookscorey.hh +++ b/opm/material/fluidmatrixinteractions/2p/brookscorey.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BrooksCorey */ -#ifndef EWOMS_BROOKS_COREY_HH -#define EWOMS_BROOKS_COREY_HH +#ifndef OPM_BROOKS_COREY_HH +#define OPM_BROOKS_COREY_HH #include "brookscoreyparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/brookscoreyparams.hh b/opm/material/fluidmatrixinteractions/2p/brookscoreyparams.hh index 988950051..30b849e20 100644 --- a/opm/material/fluidmatrixinteractions/2p/brookscoreyparams.hh +++ b/opm/material/fluidmatrixinteractions/2p/brookscoreyparams.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::BrooksCoreyParams */ -#ifndef EWOMS_BROOKS_COREY_PARAMS_HH -#define EWOMS_BROOKS_COREY_PARAMS_HH +#ifndef OPM_BROOKS_COREY_PARAMS_HH +#define OPM_BROOKS_COREY_PARAMS_HH #include diff --git a/opm/material/fluidmatrixinteractions/2p/efftoabslaw.hh b/opm/material/fluidmatrixinteractions/2p/efftoabslaw.hh index 9bf58d74d..9f2835fb0 100644 --- a/opm/material/fluidmatrixinteractions/2p/efftoabslaw.hh +++ b/opm/material/fluidmatrixinteractions/2p/efftoabslaw.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::EffToAbsLaw */ -#ifndef EWOMS_EFF_TO_ABS_LAW_HH -#define EWOMS_EFF_TO_ABS_LAW_HH +#ifndef OPM_EFF_TO_ABS_LAW_HH +#define OPM_EFF_TO_ABS_LAW_HH #include "efftoabslawparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/efftoabslawparams.hh b/opm/material/fluidmatrixinteractions/2p/efftoabslawparams.hh index f64a3e3f6..7f032c7da 100644 --- a/opm/material/fluidmatrixinteractions/2p/efftoabslawparams.hh +++ b/opm/material/fluidmatrixinteractions/2p/efftoabslawparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::EffToAbsLawParams */ -#ifndef EWOMS_EFF_TO_ABS_LAW_PARAMS_HH -#define EWOMS_EFF_TO_ABS_LAW_PARAMS_HH +#ifndef OPM_EFF_TO_ABS_LAW_PARAMS_HH +#define OPM_EFF_TO_ABS_LAW_PARAMS_HH namespace Opm { /*! diff --git a/opm/material/fluidmatrixinteractions/2p/linearmaterial.hh b/opm/material/fluidmatrixinteractions/2p/linearmaterial.hh index 2ae4e58f8..160159e76 100644 --- a/opm/material/fluidmatrixinteractions/2p/linearmaterial.hh +++ b/opm/material/fluidmatrixinteractions/2p/linearmaterial.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::LinearMaterial */ -#ifndef EWOMS_LINEAR_MATERIAL_HH -#define EWOMS_LINEAR_MATERIAL_HH +#ifndef OPM_LINEAR_MATERIAL_HH +#define OPM_LINEAR_MATERIAL_HH #include "linearmaterialparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/parkerlenhard.hh b/opm/material/fluidmatrixinteractions/2p/parkerlenhard.hh index 9b24b8069..d26d2903f 100644 --- a/opm/material/fluidmatrixinteractions/2p/parkerlenhard.hh +++ b/opm/material/fluidmatrixinteractions/2p/parkerlenhard.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::ParkerLenhard */ -#ifndef EWOMS_PARKER_LENHARD_HH -#define EWOMS_PARKER_LENHARD_HH +#ifndef OPM_PARKER_LENHARD_HH +#define OPM_PARKER_LENHARD_HH #include "parkerlenhardparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/parkerlenhardparams.hh b/opm/material/fluidmatrixinteractions/2p/parkerlenhardparams.hh index c42947777..2907fee64 100644 --- a/opm/material/fluidmatrixinteractions/2p/parkerlenhardparams.hh +++ b/opm/material/fluidmatrixinteractions/2p/parkerlenhardparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::ParkerLenhardParams */ -#ifndef EWOMS_PARKER_LENHARD_PARAMS_HH -#define EWOMS_PARKER_LENHARD_PARAMS_HH +#ifndef OPM_PARKER_LENHARD_PARAMS_HH +#define OPM_PARKER_LENHARD_PARAMS_HH #include diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh b/opm/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh index 748b716a1..131acab38 100644 --- a/opm/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh +++ b/opm/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::RegularizedBrooksCoreyParams */ -#ifndef EWOMS_REGULARIZED_BROOKS_COREY_PARAMS_HH -#define EWOMS_REGULARIZED_BROOKS_COREY_PARAMS_HH +#ifndef OPM_REGULARIZED_BROOKS_COREY_PARAMS_HH +#define OPM_REGULARIZED_BROOKS_COREY_PARAMS_HH #include "brookscoreyparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh b/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh index e36a3683a..e37b4b73e 100644 --- a/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh +++ b/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::RegularizedLinearMaterial */ -#ifndef EWOMS_REGULARIZED_LINEAR_MATERIAL_HH -#define EWOMS_REGULARIZED_LINEAR_MATERIAL_HH +#ifndef OPM_REGULARIZED_LINEAR_MATERIAL_HH +#define OPM_REGULARIZED_LINEAR_MATERIAL_HH #include "linearmaterial.hh" #include "regularizedlinearmaterialparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh index b10dbbccc..eb555e7a9 100644 --- a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh +++ b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::RegularizedVanGenuchten */ -#ifndef EWOMS_REGULARIZED_VAN_GENUCHTEN_HH -#define EWOMS_REGULARIZED_VAN_GENUCHTEN_HH +#ifndef OPM_REGULARIZED_VAN_GENUCHTEN_HH +#define OPM_REGULARIZED_VAN_GENUCHTEN_HH #include "vangenuchten.hh" #include "regularizedvangenuchtenparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh index 85f1f8c1f..3a48553de 100644 --- a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh +++ b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::RegularizedVanGenuchtenParams */ -#ifndef EWOMS_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH -#define EWOMS_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH +#ifndef OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH +#define OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH #include "vangenuchtenparams.hh" diff --git a/opm/material/fluidmatrixinteractions/2p/vangenuchten.hh b/opm/material/fluidmatrixinteractions/2p/vangenuchten.hh index 5c18d3de4..ca260f1a9 100644 --- a/opm/material/fluidmatrixinteractions/2p/vangenuchten.hh +++ b/opm/material/fluidmatrixinteractions/2p/vangenuchten.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::VanGenuchten */ -#ifndef EWOMS_VAN_GENUCHTEN_HH -#define EWOMS_VAN_GENUCHTEN_HH +#ifndef OPM_VAN_GENUCHTEN_HH +#define OPM_VAN_GENUCHTEN_HH #include "vangenuchtenparams.hh" diff --git a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh index 9620ec9d5..9239d7264 100644 --- a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh +++ b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::ThreePParkerVanGenuchten */ -#ifndef EWOMS_3P_PARKER_VAN_GENUCHTEN_HH -#define EWOMS_3P_PARKER_VAN_GENUCHTEN_HH +#ifndef OPM_3P_PARKER_VAN_GENUCHTEN_HH +#define OPM_3P_PARKER_VAN_GENUCHTEN_HH #include "3pparkervangenuchtenparams.hh" @@ -61,8 +61,7 @@ public: Sw = wetting phase saturation, or, sum of wetting phase saturations alpha : VanGenuchten-alpha - this function is just copied from MUFTE/pml/constrel3p3cni.c - that is why variable names do not yet fulfill eWoms rules, TODO Change */ + this function is copied from MUFTE/pml/constrel3p3cni.c */ Scalar r,Se,x,vg_m; Scalar pc,pc_prime,Se_regu; @@ -105,8 +104,7 @@ public: Sw = wetting phase saturation, or, sum of wetting phase saturations alpha : VanGenuchten-alpha - this function is just copied from MUFTE/pml/constrel3p3cni.c - that is why variable names do not yet fulfill eWoms rules, TODO Change */ + this function is just copied from MUFTE/pml/constrel3p3cni.c */ Scalar r,Se,x,vg_m; Scalar pc,pc_prime,Se_regu; @@ -148,8 +146,7 @@ public: /* St = sum of wetting (liquid) phase saturations alpha : VanGenuchten-alpha - this function is just copied from MUFTE/pml/constrel3p3cni.c - that is why variable names do not yet fulfill eWoms rules, TODO Change */ + this function is just copied from MUFTE/pml/constrel3p3cni.c */ Scalar r,Se,x,vg_m; Scalar pc,pc_prime,Se_regu; diff --git a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh index 6ed0246ae..a414e1d28 100644 --- a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh +++ b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh @@ -22,8 +22,8 @@ * \file * \copydoc Opm::ParkerVanGen3PParams */ -#ifndef EWOMS_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH -#define EWOMS_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH +#ifndef OPM_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH +#define OPM_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH #include diff --git a/opm/material/fluidmatrixinteractions/mp/2padapter.hh b/opm/material/fluidmatrixinteractions/mp/2padapter.hh index d59fa2e44..e44cfe880 100644 --- a/opm/material/fluidmatrixinteractions/mp/2padapter.hh +++ b/opm/material/fluidmatrixinteractions/mp/2padapter.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::TwoPAdapter */ -#ifndef EWOMS_MP_2P_ADAPTER_HH -#define EWOMS_MP_2P_ADAPTER_HH +#ifndef OPM_MP_2P_ADAPTER_HH +#define OPM_MP_2P_ADAPTER_HH #include diff --git a/opm/material/fluidmatrixinteractions/mp/3padapter.hh b/opm/material/fluidmatrixinteractions/mp/3padapter.hh index c491e9d19..877b712a6 100644 --- a/opm/material/fluidmatrixinteractions/mp/3padapter.hh +++ b/opm/material/fluidmatrixinteractions/mp/3padapter.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::ThreePAdapter */ -#ifndef EWOMS_MP_3P_ADAPTER_HH -#define EWOMS_MP_3P_ADAPTER_HH +#ifndef OPM_MP_3P_ADAPTER_HH +#define OPM_MP_3P_ADAPTER_HH #include diff --git a/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh b/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh index ee567c716..233796873 100644 --- a/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh +++ b/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::MpLinearMaterial */ -#ifndef EWOMS_MP_LINEAR_MATERIAL_HH -#define EWOMS_MP_LINEAR_MATERIAL_HH +#ifndef OPM_MP_LINEAR_MATERIAL_HH +#define OPM_MP_LINEAR_MATERIAL_HH #include "mplinearmaterialparams.hh" diff --git a/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh b/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh index fac4b2dd2..e7008cab7 100644 --- a/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh +++ b/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::NullMaterialLaw */ -#ifndef EWOMS_NULL_MATERIAL_LAW_HH -#define EWOMS_NULL_MATERIAL_LAW_HH +#ifndef OPM_NULL_MATERIAL_LAW_HH +#define OPM_NULL_MATERIAL_LAW_HH #include "nullmateriallawparams.hh" diff --git a/opm/material/fluidmatrixinteractions/mp/nullmateriallawparams.hh b/opm/material/fluidmatrixinteractions/mp/nullmateriallawparams.hh index 17dfbac4a..26f0c7c1b 100644 --- a/opm/material/fluidmatrixinteractions/mp/nullmateriallawparams.hh +++ b/opm/material/fluidmatrixinteractions/mp/nullmateriallawparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::NullMaterialLawParams */ -#ifndef EWOMS_NULL_MATERIAL_LAW_PARAMS_HH -#define EWOMS_NULL_MATERIAL_LAW_PARAMS_HH +#ifndef OPM_NULL_MATERIAL_LAW_PARAMS_HH +#define OPM_NULL_MATERIAL_LAW_PARAMS_HH namespace Opm { /*! diff --git a/opm/material/fluidstates/compositionalfluidstate.hh b/opm/material/fluidstates/compositionalfluidstate.hh index f583ef189..ef5b8a582 100644 --- a/opm/material/fluidstates/compositionalfluidstate.hh +++ b/opm/material/fluidstates/compositionalfluidstate.hh @@ -23,8 +23,8 @@ * multi-phase, multi-component fluid system assuming * thermodynamic equilibrium. */ -#ifndef EWOMS_COMPOSITIONAL_FLUID_STATE_HH -#define EWOMS_COMPOSITIONAL_FLUID_STATE_HH +#ifndef OPM_COMPOSITIONAL_FLUID_STATE_HH +#define OPM_COMPOSITIONAL_FLUID_STATE_HH #include "modularfluidstate.hh" diff --git a/opm/material/fluidstates/fluidstatecompositionmodules.hh b/opm/material/fluidstates/fluidstatecompositionmodules.hh index ff586b9ce..0618ffd8a 100644 --- a/opm/material/fluidstates/fluidstatecompositionmodules.hh +++ b/opm/material/fluidstates/fluidstatecompositionmodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent composition. */ -#ifndef EWOMS_FLUID_STATE_COMPOSITION_MODULES_HH -#define EWOMS_FLUID_STATE_COMPOSITION_MODULES_HH +#ifndef OPM_FLUID_STATE_COMPOSITION_MODULES_HH +#define OPM_FLUID_STATE_COMPOSITION_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstatedensitymodules.hh b/opm/material/fluidstates/fluidstatedensitymodules.hh index 843b548af..b6c8d1719 100644 --- a/opm/material/fluidstates/fluidstatedensitymodules.hh +++ b/opm/material/fluidstates/fluidstatedensitymodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent density. */ -#ifndef EWOMS_FLUID_STATE_DENSITY_MODULES_HH -#define EWOMS_FLUID_STATE_DENSITY_MODULES_HH +#ifndef OPM_FLUID_STATE_DENSITY_MODULES_HH +#define OPM_FLUID_STATE_DENSITY_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstateenthalpymodules.hh b/opm/material/fluidstates/fluidstateenthalpymodules.hh index 6ae855fb9..1fae17aca 100644 --- a/opm/material/fluidstates/fluidstateenthalpymodules.hh +++ b/opm/material/fluidstates/fluidstateenthalpymodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent enthalpy. */ -#ifndef EWOMS_FLUID_STATE_ENTHALPY_MODULES_HH -#define EWOMS_FLUID_STATE_ENTHALPY_MODULES_HH +#ifndef OPM_FLUID_STATE_ENTHALPY_MODULES_HH +#define OPM_FLUID_STATE_ENTHALPY_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstatefugacitymodules.hh b/opm/material/fluidstates/fluidstatefugacitymodules.hh index dff2211ea..5b228efdd 100644 --- a/opm/material/fluidstates/fluidstatefugacitymodules.hh +++ b/opm/material/fluidstates/fluidstatefugacitymodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent fugacity/chemical potential. */ -#ifndef EWOMS_FLUID_STATE_FUGACITY_MODULES_HH -#define EWOMS_FLUID_STATE_FUGACITY_MODULES_HH +#ifndef OPM_FLUID_STATE_FUGACITY_MODULES_HH +#define OPM_FLUID_STATE_FUGACITY_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstatepressuremodules.hh b/opm/material/fluidstates/fluidstatepressuremodules.hh index 21f94183d..e12102074 100644 --- a/opm/material/fluidstates/fluidstatepressuremodules.hh +++ b/opm/material/fluidstates/fluidstatepressuremodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent pressure. */ -#ifndef EWOMS_FLUID_STATE_PRESSURE_MODULES_HH -#define EWOMS_FLUID_STATE_PRESSURE_MODULES_HH +#ifndef OPM_FLUID_STATE_PRESSURE_MODULES_HH +#define OPM_FLUID_STATE_PRESSURE_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstatesaturationmodules.hh b/opm/material/fluidstates/fluidstatesaturationmodules.hh index 903d81942..09c32176c 100644 --- a/opm/material/fluidstates/fluidstatesaturationmodules.hh +++ b/opm/material/fluidstates/fluidstatesaturationmodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent saturation. */ -#ifndef EWOMS_FLUID_STATE_SATURATION_MODULES_HH -#define EWOMS_FLUID_STATE_SATURATION_MODULES_HH +#ifndef OPM_FLUID_STATE_SATURATION_MODULES_HH +#define OPM_FLUID_STATE_SATURATION_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstatetemperaturemodules.hh b/opm/material/fluidstates/fluidstatetemperaturemodules.hh index 1efee1d39..1a9873f0c 100644 --- a/opm/material/fluidstates/fluidstatetemperaturemodules.hh +++ b/opm/material/fluidstates/fluidstatetemperaturemodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent temperature. */ -#ifndef EWOMS_FLUID_STATE_TEMPERATURE_MODULES_HH -#define EWOMS_FLUID_STATE_TEMPERATURE_MODULES_HH +#ifndef OPM_FLUID_STATE_TEMPERATURE_MODULES_HH +#define OPM_FLUID_STATE_TEMPERATURE_MODULES_HH #include diff --git a/opm/material/fluidstates/fluidstateviscositymodules.hh b/opm/material/fluidstates/fluidstateviscositymodules.hh index fdeb8d164..dfbfa7d60 100644 --- a/opm/material/fluidstates/fluidstateviscositymodules.hh +++ b/opm/material/fluidstates/fluidstateviscositymodules.hh @@ -21,8 +21,8 @@ * * \brief Modules for the ModularFluidState which represent viscosity. */ -#ifndef EWOMS_FLUID_STATE_VISCOSITY_MODULES_HH -#define EWOMS_FLUID_STATE_VISCOSITY_MODULES_HH +#ifndef OPM_FLUID_STATE_VISCOSITY_MODULES_HH +#define OPM_FLUID_STATE_VISCOSITY_MODULES_HH #include diff --git a/opm/material/fluidstates/immisciblefluidstate.hh b/opm/material/fluidstates/immisciblefluidstate.hh index 49768bd25..e3b1c0331 100644 --- a/opm/material/fluidstates/immisciblefluidstate.hh +++ b/opm/material/fluidstates/immisciblefluidstate.hh @@ -23,8 +23,8 @@ * multi-phase, multi-component fluid system assuming * thermodynamic equilibrium. */ -#ifndef EWOMS_IMMISCIBLE_FLUID_STATE_HH -#define EWOMS_IMMISCIBLE_FLUID_STATE_HH +#ifndef OPM_IMMISCIBLE_FLUID_STATE_HH +#define OPM_IMMISCIBLE_FLUID_STATE_HH #include "modularfluidstate.hh" diff --git a/opm/material/fluidstates/modularfluidstate.hh b/opm/material/fluidstates/modularfluidstate.hh index dab6bb714..3a57e50fc 100644 --- a/opm/material/fluidstates/modularfluidstate.hh +++ b/opm/material/fluidstates/modularfluidstate.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::ModularFluidState */ -#ifndef EWOMS_MODULAR_FLUID_STATE_HH -#define EWOMS_MODULAR_FLUID_STATE_HH +#ifndef OPM_MODULAR_FLUID_STATE_HH +#define OPM_MODULAR_FLUID_STATE_HH #include "fluidstatepressuremodules.hh" #include "fluidstatetemperaturemodules.hh" diff --git a/opm/material/fluidstates/nonequilibriumfluidstate.hh b/opm/material/fluidstates/nonequilibriumfluidstate.hh index b4aac1cf8..a0bcfc089 100644 --- a/opm/material/fluidstates/nonequilibriumfluidstate.hh +++ b/opm/material/fluidstates/nonequilibriumfluidstate.hh @@ -23,8 +23,8 @@ * multi-phase, multi-component fluid system _not_ assuming * thermodynamic equilibrium. */ -#ifndef EWOMS_NON_EQUILIBRIUM_FLUID_STATE_HH -#define EWOMS_NON_EQUILIBRIUM_FLUID_STATE_HH +#ifndef OPM_NON_EQUILIBRIUM_FLUID_STATE_HH +#define OPM_NON_EQUILIBRIUM_FLUID_STATE_HH #include "modularfluidstate.hh" diff --git a/opm/material/fluidstates/pressureoverlayfluidstate.hh b/opm/material/fluidstates/pressureoverlayfluidstate.hh index 6b4a63129..498d306f5 100644 --- a/opm/material/fluidstates/pressureoverlayfluidstate.hh +++ b/opm/material/fluidstates/pressureoverlayfluidstate.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::PressureOverlayFluidState */ -#ifndef EWOMS_PRESSURE_OVERLAY_FLUID_STATE_HH -#define EWOMS_PRESSURE_OVERLAY_FLUID_STATE_HH +#ifndef OPM_PRESSURE_OVERLAY_FLUID_STATE_HH +#define OPM_PRESSURE_OVERLAY_FLUID_STATE_HH #include diff --git a/opm/material/fluidstates/saturationoverlayfluidstate.hh b/opm/material/fluidstates/saturationoverlayfluidstate.hh index a075b88b8..d17fc9e9f 100644 --- a/opm/material/fluidstates/saturationoverlayfluidstate.hh +++ b/opm/material/fluidstates/saturationoverlayfluidstate.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::SaturationOverlayFluidState */ -#ifndef EWOMS_SATURATION_OVERLAY_FLUID_STATE_HH -#define EWOMS_SATURATION_OVERLAY_FLUID_STATE_HH +#ifndef OPM_SATURATION_OVERLAY_FLUID_STATE_HH +#define OPM_SATURATION_OVERLAY_FLUID_STATE_HH #include diff --git a/opm/material/fluidstates/temperatureoverlayfluidstate.hh b/opm/material/fluidstates/temperatureoverlayfluidstate.hh index baf3d15c6..d5f921502 100644 --- a/opm/material/fluidstates/temperatureoverlayfluidstate.hh +++ b/opm/material/fluidstates/temperatureoverlayfluidstate.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::TemperatureOverlayFluidState */ -#ifndef EWOMS_TEMPERATURE_OVERLAY_FLUID_STATE_HH -#define EWOMS_TEMPERATURE_OVERLAY_FLUID_STATE_HH +#ifndef OPM_TEMPERATURE_OVERLAY_FLUID_STATE_HH +#define OPM_TEMPERATURE_OVERLAY_FLUID_STATE_HH #include diff --git a/opm/material/fluidsystems/1pfluidsystem.hh b/opm/material/fluidsystems/1pfluidsystem.hh index 0b54537eb..c26ca02f8 100644 --- a/opm/material/fluidsystems/1pfluidsystem.hh +++ b/opm/material/fluidsystems/1pfluidsystem.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidSystems::OneP */ -#ifndef EWOMS_1P_FLUIDSYSTEM_HH -#define EWOMS_1P_FLUIDSYSTEM_HH +#ifndef OPM_1P_FLUIDSYSTEM_HH +#define OPM_1P_FLUIDSYSTEM_HH #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/2pimmisciblefluidsystem.hh b/opm/material/fluidsystems/2pimmisciblefluidsystem.hh index 66af5c5c5..ea21511d9 100644 --- a/opm/material/fluidsystems/2pimmisciblefluidsystem.hh +++ b/opm/material/fluidsystems/2pimmisciblefluidsystem.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::FluidSystems::TwoPImmiscible */ -#ifndef EWOMS_2P_IMMISCIBLE_FLUID_SYSTEM_HH -#define EWOMS_2P_IMMISCIBLE_FLUID_SYSTEM_HH +#ifndef OPM_2P_IMMISCIBLE_FLUID_SYSTEM_HH +#define OPM_2P_IMMISCIBLE_FLUID_SYSTEM_HH #include #include diff --git a/opm/material/fluidsystems/basefluidsystem.hh b/opm/material/fluidsystems/basefluidsystem.hh index ab0fa4958..58d5dd630 100644 --- a/opm/material/fluidsystems/basefluidsystem.hh +++ b/opm/material/fluidsystems/basefluidsystem.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::BaseFluidSystem */ -#ifndef EWOMS_BASE_FLUID_SYSTEM_HH -#define EWOMS_BASE_FLUID_SYSTEM_HH +#ifndef OPM_BASE_FLUID_SYSTEM_HH +#define OPM_BASE_FLUID_SYSTEM_HH #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/blackoilfluidsystem.hh b/opm/material/fluidsystems/blackoilfluidsystem.hh index 2c3501266..4e0b7951b 100644 --- a/opm/material/fluidsystems/blackoilfluidsystem.hh +++ b/opm/material/fluidsystems/blackoilfluidsystem.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidSystems::BlackOil */ -#ifndef EWOMS_BLACK_OIL_FLUID_SYSTEM_HH -#define EWOMS_BLACK_OIL_FLUID_SYSTEM_HH +#ifndef OPM_BLACK_OIL_FLUID_SYSTEM_HH +#define OPM_BLACK_OIL_FLUID_SYSTEM_HH #include #include diff --git a/opm/material/fluidsystems/brineco2fluidsystem.hh b/opm/material/fluidsystems/brineco2fluidsystem.hh index 2faabaa1b..8c4d39a96 100644 --- a/opm/material/fluidsystems/brineco2fluidsystem.hh +++ b/opm/material/fluidsystems/brineco2fluidsystem.hh @@ -21,8 +21,8 @@ * * \copydoc Opm::FluidSystems::BrineCO2 */ -#ifndef EWOMS_BRINE_CO2_SYSTEM_HH -#define EWOMS_BRINE_CO2_SYSTEM_HH +#ifndef OPM_BRINE_CO2_SYSTEM_HH +#define OPM_BRINE_CO2_SYSTEM_HH #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/gasphase.hh b/opm/material/fluidsystems/gasphase.hh index 84d5ac811..be009e17f 100644 --- a/opm/material/fluidsystems/gasphase.hh +++ b/opm/material/fluidsystems/gasphase.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::GasPhase */ -#ifndef EWOMS_GAS_PHASE_HH -#define EWOMS_GAS_PHASE_HH +#ifndef OPM_GAS_PHASE_HH +#define OPM_GAS_PHASE_HH namespace Opm { diff --git a/opm/material/fluidsystems/h2oairfluidsystem.hh b/opm/material/fluidsystems/h2oairfluidsystem.hh index 0317b24fb..374375f0e 100644 --- a/opm/material/fluidsystems/h2oairfluidsystem.hh +++ b/opm/material/fluidsystems/h2oairfluidsystem.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::FluidSystems::H2OAir */ -#ifndef EWOMS_H2O_AIR_SYSTEM_HH -#define EWOMS_H2O_AIR_SYSTEM_HH +#ifndef OPM_H2O_AIR_SYSTEM_HH +#define OPM_H2O_AIR_SYSTEM_HH #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh b/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh index 136739ecd..cfbcf5396 100644 --- a/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh +++ b/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidSystems::H2OAirMesitylene */ -#ifndef EWOMS_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH -#define EWOMS_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH +#ifndef OPM_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH +#define OPM_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/h2oairxylenefluidsystem.hh b/opm/material/fluidsystems/h2oairxylenefluidsystem.hh index bf6e66c08..6a5c6b35f 100644 --- a/opm/material/fluidsystems/h2oairxylenefluidsystem.hh +++ b/opm/material/fluidsystems/h2oairxylenefluidsystem.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::FluidSystems::H2OAirXylene */ -#ifndef EWOMS_H2O_AIR_XYLENE_FLUID_SYSTEM_HH -#define EWOMS_H2O_AIR_XYLENE_FLUID_SYSTEM_HH +#ifndef OPM_H2O_AIR_XYLENE_FLUID_SYSTEM_HH +#define OPM_H2O_AIR_XYLENE_FLUID_SYSTEM_HH #include #include diff --git a/opm/material/fluidsystems/h2on2fluidsystem.hh b/opm/material/fluidsystems/h2on2fluidsystem.hh index 1a5d43997..93c98d35f 100644 --- a/opm/material/fluidsystems/h2on2fluidsystem.hh +++ b/opm/material/fluidsystems/h2on2fluidsystem.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidSystems::H2ON2 */ -#ifndef EWOMS_H2O_N2_FLUID_SYSTEM_HH -#define EWOMS_H2O_N2_FLUID_SYSTEM_HH +#ifndef OPM_H2O_N2_FLUID_SYSTEM_HH +#define OPM_H2O_N2_FLUID_SYSTEM_HH #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh b/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh index 6da3e8099..0c84cc16d 100644 --- a/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh +++ b/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::FluidSystems::H2ON2LiquidPhase */ -#ifndef EWOMS_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH -#define EWOMS_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH +#ifndef OPM_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH +#define OPM_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/liquidphase.hh b/opm/material/fluidsystems/liquidphase.hh index 74fe4f06b..24cace9d6 100644 --- a/opm/material/fluidsystems/liquidphase.hh +++ b/opm/material/fluidsystems/liquidphase.hh @@ -21,8 +21,8 @@ * \file * \copydoc Opm::LiquidPhase */ -#ifndef EWOMS_LIQUID_PHASE_HH -#define EWOMS_LIQUID_PHASE_HH +#ifndef OPM_LIQUID_PHASE_HH +#define OPM_LIQUID_PHASE_HH namespace Opm { diff --git a/opm/material/fluidsystems/nullparametercache.hh b/opm/material/fluidsystems/nullparametercache.hh index bafb6d71d..7e780a21e 100644 --- a/opm/material/fluidsystems/nullparametercache.hh +++ b/opm/material/fluidsystems/nullparametercache.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::NullParameterCache */ -#ifndef EWOMS_NULL_PARAMETER_CACHE_HH -#define EWOMS_NULL_PARAMETER_CACHE_HH +#ifndef OPM_NULL_PARAMETER_CACHE_HH +#define OPM_NULL_PARAMETER_CACHE_HH #include "parametercachebase.hh" diff --git a/opm/material/fluidsystems/parametercachebase.hh b/opm/material/fluidsystems/parametercachebase.hh index e4895e381..63ccb0da8 100644 --- a/opm/material/fluidsystems/parametercachebase.hh +++ b/opm/material/fluidsystems/parametercachebase.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::ParameterCacheBase */ -#ifndef EWOMS_PARAMETER_CACHE_BASE_HH -#define EWOMS_PARAMETER_CACHE_BASE_HH +#ifndef OPM_PARAMETER_CACHE_BASE_HH +#define OPM_PARAMETER_CACHE_BASE_HH namespace Opm { diff --git a/opm/material/fluidsystems/spe5fluidsystem.hh b/opm/material/fluidsystems/spe5fluidsystem.hh index a5e9bea6a..babe29ab6 100644 --- a/opm/material/fluidsystems/spe5fluidsystem.hh +++ b/opm/material/fluidsystems/spe5fluidsystem.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidSystems::Spe5 */ -#ifndef EWOMS_SPE5_FLUID_SYSTEM_HH -#define EWOMS_SPE5_FLUID_SYSTEM_HH +#ifndef OPM_SPE5_FLUID_SYSTEM_HH +#define OPM_SPE5_FLUID_SYSTEM_HH #include "basefluidsystem.hh" #include "spe5parametercache.hh" diff --git a/opm/material/fluidsystems/spe5parametercache.hh b/opm/material/fluidsystems/spe5parametercache.hh index 56635ccd0..555ca9c5d 100644 --- a/opm/material/fluidsystems/spe5parametercache.hh +++ b/opm/material/fluidsystems/spe5parametercache.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::Spe5ParameterCache */ -#ifndef EWOMS_SPE5_PARAMETER_CACHE_HH -#define EWOMS_SPE5_PARAMETER_CACHE_HH +#ifndef OPM_SPE5_PARAMETER_CACHE_HH +#define OPM_SPE5_PARAMETER_CACHE_HH #include diff --git a/opm/material/heatconduction/dummyheatconductionlaw.hh b/opm/material/heatconduction/dummyheatconductionlaw.hh index b0eeccb11..16a4c6a08 100644 --- a/opm/material/heatconduction/dummyheatconductionlaw.hh +++ b/opm/material/heatconduction/dummyheatconductionlaw.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::DummyHeatConductionLaw */ -#ifndef EWOMS_DUMMY_HEATCONDUCTION_LAW_HH -#define EWOMS_DUMMY_HEATCONDUCTION_LAW_HH +#ifndef OPM_DUMMY_HEATCONDUCTION_LAW_HH +#define OPM_DUMMY_HEATCONDUCTION_LAW_HH #include diff --git a/opm/material/heatconduction/fluidconduction.hh b/opm/material/heatconduction/fluidconduction.hh index 8793e1f41..963230db0 100644 --- a/opm/material/heatconduction/fluidconduction.hh +++ b/opm/material/heatconduction/fluidconduction.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidHeatConduction */ -#ifndef EWOMS_FLUID_HEAT_CONDUCTION_HH -#define EWOMS_FLUID_HEAT_CONDUCTION_HH +#ifndef OPM_FLUID_HEAT_CONDUCTION_HH +#define OPM_FLUID_HEAT_CONDUCTION_HH #include "fluidconductionparams.hh" diff --git a/opm/material/heatconduction/fluidconductionparams.hh b/opm/material/heatconduction/fluidconductionparams.hh index a6b84bd6a..cf92016c7 100644 --- a/opm/material/heatconduction/fluidconductionparams.hh +++ b/opm/material/heatconduction/fluidconductionparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::FluidHeatConductionParams */ -#ifndef EWOMS_FLUID_HEAT_CONDUCTION_PARAMS_HH -#define EWOMS_FLUID_HEAT_CONDUCTION_PARAMS_HH +#ifndef OPM_FLUID_HEAT_CONDUCTION_PARAMS_HH +#define OPM_FLUID_HEAT_CONDUCTION_PARAMS_HH namespace Opm { /*! diff --git a/opm/material/heatconduction/somerton.hh b/opm/material/heatconduction/somerton.hh index 264d49700..1a50a33ea 100644 --- a/opm/material/heatconduction/somerton.hh +++ b/opm/material/heatconduction/somerton.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::Somerton */ -#ifndef EWOMS_SOMERTON_HH -#define EWOMS_SOMERTON_HH +#ifndef OPM_SOMERTON_HH +#define OPM_SOMERTON_HH #include "somertonparams.hh" diff --git a/opm/material/heatconduction/somertonparams.hh b/opm/material/heatconduction/somertonparams.hh index 314e5dc9f..8927e4047 100644 --- a/opm/material/heatconduction/somertonparams.hh +++ b/opm/material/heatconduction/somertonparams.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::SomertonParams */ -#ifndef EWOMS_SOMERTON_PARAMS_HH -#define EWOMS_SOMERTON_PARAMS_HH +#ifndef OPM_SOMERTON_PARAMS_HH +#define OPM_SOMERTON_PARAMS_HH #include diff --git a/opm/material/idealgas.hh b/opm/material/idealgas.hh index 05e7284f1..ecff9c21a 100644 --- a/opm/material/idealgas.hh +++ b/opm/material/idealgas.hh @@ -20,8 +20,8 @@ * \file * \copydoc Opm::IdealGas */ -#ifndef EWOMS_IDEAL_GAS_HH -#define EWOMS_IDEAL_GAS_HH +#ifndef OPM_IDEAL_GAS_HH +#define OPM_IDEAL_GAS_HH #include diff --git a/tests/material/CMakeLists.txt b/tests/material/CMakeLists.txt new file mode 100644 index 000000000..4edc99707 --- /dev/null +++ b/tests/material/CMakeLists.txt @@ -0,0 +1,6 @@ +add_subdirectory("fluidsystems") +add_subdirectory("immiscibleflash") +add_subdirectory("ncpflash") +add_subdirectory("pengrobinson") +add_subdirectory("tabulation") + diff --git a/tests/material/fluidsystems/CMakeLists.txt b/tests/material/fluidsystems/CMakeLists.txt new file mode 100644 index 000000000..be8a4aeb6 --- /dev/null +++ b/tests/material/fluidsystems/CMakeLists.txt @@ -0,0 +1,8 @@ +# set the CMake module include path +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules") + +# eWoms unit-testing infrastructure +include(EwomsAddTest) + +EwomsAddTest(test_fluidsystems DRIVER_ARGS --plain test_fluidsystems) + diff --git a/tests/material/fluidsystems/checkfluidsystem.hh b/tests/material/fluidsystems/checkfluidsystem.hh new file mode 100644 index 000000000..9cb642478 --- /dev/null +++ b/tests/material/fluidsystems/checkfluidsystem.hh @@ -0,0 +1,328 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2011-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc checkFluidSystem + */ +#ifndef OPM_CHECK_FLUIDSYSTEM_HH +#define OPM_CHECK_FLUIDSYSTEM_HH + +// include all fluid systems in opm-material +#include +#include +#include +#include +#include +#include +#include +#include + +// include all fluid states +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +/*! + * \brief This is a fluid state which makes sure that only the quantities + * allowed are accessed. + */ +template > +class HairSplittingFluidState + : protected BaseFluidState +{ +public: + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + + HairSplittingFluidState() + { + // initially, do not allow anything + allowTemperature(false); + allowPressure(false); + allowComposition(false); + allowDensity(false); + + // do not allow accessing any phase + restrictToPhase(1000); + } + + void allowTemperature(bool yesno) + { allowTemperature_ = yesno; } + + void allowPressure(bool yesno) + { allowPressure_ = yesno; } + + void allowComposition(bool yesno) + { allowComposition_ = yesno; } + + void allowDensity(bool yesno) + { allowDensity_ = yesno; } + + void restrictToPhase(int phaseIdx) + { restrictPhaseIdx_ = phaseIdx; } + + Scalar temperature(int phaseIdx) const + { + assert(allowTemperature_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::temperature(phaseIdx); + return 1e100; + } + + Scalar pressure(int phaseIdx) const + { + assert(allowPressure_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::pressure(phaseIdx); + return 1e100; + } + + Scalar moleFraction(int phaseIdx, int compIdx) const + { + assert(allowComposition_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::moleFraction(phaseIdx, compIdx); + return 1e100; + } + + Scalar massFraction(int phaseIdx, int compIdx) const + { + assert(allowComposition_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::massFraction(phaseIdx, compIdx); + return 1e100; + } + + Scalar averageMolarMass(int phaseIdx) const + { + assert(allowComposition_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::averageMolarMass(phaseIdx); + return 1e100; + } + + Scalar molarity(int phaseIdx, int compIdx) const + { + assert(allowDensity_ && allowComposition_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::molarity(phaseIdx, compIdx); + return 1e100; + } + + Scalar molarDensity(int phaseIdx) const + { + assert(allowDensity_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::molarDensity(phaseIdx); + return 1e100; + } + + Scalar molarVolume(int phaseIdx) const + { + assert(allowDensity_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::molarVolume(phaseIdx); + return 1e100; + } + + Scalar density(int phaseIdx) const + { + assert(allowDensity_); + assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); + DUNE_UNUSED Scalar tmp = BaseFluidState::density(phaseIdx); + return 1e100; + } + + Scalar saturation(int phaseIdx) const + { + assert(false); + DUNE_UNUSED Scalar tmp = BaseFluidState::saturation(phaseIdx); + return 1e100; + } + + Scalar fugacity(int phaseIdx, int compIdx) const + { + assert(false); + DUNE_UNUSED Scalar tmp = BaseFluidState::fugacity(phaseIdx, compIdx); + return 1e100; + } + + Scalar fugacityCoefficient(int phaseIdx, int compIdx) const + { + assert(false); + DUNE_UNUSED Scalar tmp = BaseFluidState::fugacityCoefficient(phaseIdx, compIdx); + return 1e100; + } + + Scalar enthalpy(int phaseIdx) const + { + assert(false); + DUNE_UNUSED Scalar tmp = BaseFluidState::enthalpy(phaseIdx); + return 1e100; + } + + Scalar internalEnergy(int phaseIdx) const + { + assert(false); + DUNE_UNUSED Scalar tmp = BaseFluidState::internalEnergy(phaseIdx); + return 1e100; + } + + Scalar viscosity(int phaseIdx) const + { + assert(false); + DUNE_UNUSED Scalar tmp = BaseFluidState::viscosity(phaseIdx); + return 1e100; + } + +private: + bool allowSaturation_; + bool allowTemperature_; + bool allowPressure_; + bool allowComposition_; + bool allowDensity_; + int restrictPhaseIdx_; +}; + +template +void checkFluidState(const BaseFluidState &fs) +{ + // fluid states must be copy-able + BaseFluidState tmpFs(fs); + tmpFs = fs; + + // a fluid state must provide a checkDefined() method + fs.checkDefined(); + + // make sure the fluid state provides all mandatory methods + while (false) { + Scalar DUNE_UNUSED val; + + val = fs.temperature(/*phaseIdx=*/0); + val = fs.pressure(/*phaseIdx=*/0); + val = fs.moleFraction(/*phaseIdx=*/0, /*compIdx=*/0); + val = fs.massFraction(/*phaseIdx=*/0, /*compIdx=*/0); + val = fs.averageMolarMass(/*phaseIdx=*/0); + val = fs.molarity(/*phaseIdx=*/0, /*compIdx=*/0); + val = fs.molarDensity(/*phaseIdx=*/0); + val = fs.molarVolume(/*phaseIdx=*/0); + val = fs.density(/*phaseIdx=*/0); + val = fs.saturation(/*phaseIdx=*/0); + val = fs.fugacity(/*phaseIdx=*/0, /*compIdx=*/0); + val = fs.fugacityCoefficient(/*phaseIdx=*/0, /*compIdx=*/0); + val = fs.enthalpy(/*phaseIdx=*/0); + val = fs.internalEnergy(/*phaseIdx=*/0); + val = fs.viscosity(/*phaseIdx=*/0); + }; +} + +/*! + * \brief Checks whether a fluid system adheres to the specification. + */ +template +void checkFluidSystem() +{ + std::cout << "Testing fluid system '" << Dune::className() << "'\n"; + + // make sure the fluid system provides the number of phases and + // the number of components + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + + HairSplittingFluidState fs; + fs.allowTemperature(true); + fs.allowPressure(true); + fs.allowComposition(true); + fs.restrictToPhase(-1); + + // check whether the parameter cache adheres to the API + typedef typename FluidSystem::ParameterCache PC; + PC paramCache; + try { paramCache.updateAll(fs); } catch (...) {}; + try { paramCache.updateAll(fs, /*except=*/PC::None); } catch (...) {}; + try { paramCache.updateAll(fs, /*except=*/PC::Temperature | PC::Pressure | PC::Composition); } catch (...) {}; + try { paramCache.updateAllPressures(fs); } catch (...) {}; + + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + fs.restrictToPhase(phaseIdx); + try { paramCache.updatePhase(fs, phaseIdx); } catch (...) {}; + try { paramCache.updatePhase(fs, phaseIdx, /*except=*/PC::None); } catch (...) {}; + try { paramCache.updatePhase(fs, phaseIdx, /*except=*/PC::Temperature | PC::Pressure | PC::Composition); } catch (...) {}; + try { paramCache.updateTemperature(fs, phaseIdx); } catch (...) {}; + try { paramCache.updatePressure(fs, phaseIdx); } catch (...) {}; + try { paramCache.updateComposition(fs, phaseIdx); } catch (...) {}; + try { paramCache.updateSingleMoleFraction(fs, phaseIdx, /*compIdx=*/0); } catch (...) {}; + } + + // some value to make sure the return values of the fluid system + // are convertible to scalars + Scalar DUNE_UNUSED val; + + // actually check the fluid system API + try { FluidSystem::init(); } catch (...) {}; + for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fs.restrictToPhase(phaseIdx); + fs.allowPressure(FluidSystem::isCompressible(phaseIdx)); + fs.allowComposition(true); + fs.allowDensity(false); + try { val = FluidSystem::density(fs, paramCache, phaseIdx); } catch (...) {}; + + fs.allowPressure(true); + fs.allowDensity(true); + try { val = FluidSystem::viscosity(fs, paramCache, phaseIdx); } catch (...) {}; + try { val = FluidSystem::enthalpy(fs, paramCache, phaseIdx); } catch (...) {}; + try { val = FluidSystem::heatCapacity(fs, paramCache, phaseIdx); } catch (...) {}; + try { val = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); } catch (...) {}; + + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { + fs.allowComposition(!FluidSystem::isIdealMixture(phaseIdx)); + try { val = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; + fs.allowComposition(true); + try { val = FluidSystem::diffusionCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; + } + + } + + // test for phaseName(), isLiquid() and isIdealGas() + for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + std::string DUNE_UNUSED name = FluidSystem::phaseName(phaseIdx); + bool DUNE_UNUSED bVal = FluidSystem::isLiquid(phaseIdx); + bVal = FluidSystem::isIdealGas(phaseIdx); + } + + // test for componentName() + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { + val = FluidSystem::molarMass(compIdx); + std::string DUNE_UNUSED name = FluidSystem::componentName(compIdx); + } + + std::cout << "----------------------------------\n"; +} + +#endif diff --git a/tests/material/fluidsystems/test_fluidsystems.cc b/tests/material/fluidsystems/test_fluidsystems.cc new file mode 100644 index 000000000..1843505ed --- /dev/null +++ b/tests/material/fluidsystems/test_fluidsystems.cc @@ -0,0 +1,176 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2011-2013 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief This test makes sure that the programming interface is + * observed by all fluid systems + */ +#include "config.h" + +#include "checkfluidsystem.hh" + +#include +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) +#include +#else +#include +#endif + +// include all fluid systems in opm-material +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// include all fluid states +#include +#include +#include +#include +#include +#include + +// include the tables for CO2 which are delivered with opm-material by +// default +#include +namespace Opm { +namespace FluidSystemsTest { +#include +} } + +int main(int argc, char **argv) +{ + typedef double Scalar; + typedef Opm::H2O H2O; + typedef Opm::N2 N2; + + typedef Opm::LiquidPhase Liquid; + typedef Opm::GasPhase Gas; + + // initialize MPI, finalize is done automatically on exit + Dune::MPIHelper::instance(argc, argv); + + // check all fluid states + { + typedef Opm::FluidSystems::H2ON2 FluidSystem; + + // CompositionalFluidState + { Opm::CompositionalFluidState fs; + checkFluidState(fs); } + + // NonEquilibriumFluidState + { Opm::NonEquilibriumFluidState fs; + checkFluidState(fs); } + + // ImmiscibleFluidState + { Opm::ImmiscibleFluidState fs; + checkFluidState(fs); } + + typedef Opm::CompositionalFluidState BaseFluidState; + BaseFluidState baseFs; + + // TemperatureOverlayFluidState + { Opm::TemperatureOverlayFluidState fs(baseFs); + checkFluidState(fs); } + + // PressureOverlayFluidState + { Opm::PressureOverlayFluidState fs(baseFs); + checkFluidState(fs); } + + // SaturationOverlayFluidState + { Opm::SaturationOverlayFluidState fs(baseFs); + checkFluidState(fs); } + } + + // black-oil + { typedef Opm::FluidSystems::BlackOil FluidSystem; + if (false) checkFluidSystem(); } + + // Brine -- CO2 + { typedef Opm::FluidSystems::BrineCO2 FluidSystem; + checkFluidSystem(); } + + // H2O -- N2 + { typedef Opm::FluidSystems::H2ON2 FluidSystem; + checkFluidSystem(); } + + { typedef Opm::FluidSystems::H2ON2 FluidSystem; + checkFluidSystem(); } + + // H2O -- N2 -- liquid phase + { typedef Opm::FluidSystems::H2ON2LiquidPhase FluidSystem; + checkFluidSystem(); } + + { typedef Opm::FluidSystems::H2ON2LiquidPhase FluidSystem; + checkFluidSystem(); } + + // H2O -- Air + { typedef Opm::SimpleH2O H2O; + const bool enableComplexRelations=false; + typedef Opm::FluidSystems::H2OAir FluidSystem; + checkFluidSystem(); } + + { typedef Opm::SimpleH2O H2O; + const bool enableComplexRelations=true; + typedef Opm::FluidSystems::H2OAir FluidSystem; + checkFluidSystem(); } + + { typedef Opm::H2O H2O; + const bool enableComplexRelations=false; + typedef Opm::FluidSystems::H2OAir FluidSystem; + checkFluidSystem(); } + + { typedef Opm::H2O H2O; + const bool enableComplexRelations=true; + typedef Opm::FluidSystems::H2OAir FluidSystem; + checkFluidSystem(); } + + // H2O -- Air -- Mesitylene + { typedef Opm::FluidSystems::H2OAirMesitylene FluidSystem; + checkFluidSystem(); } + + // H2O -- Air -- Xylene + { typedef Opm::FluidSystems::H2OAirXylene FluidSystem; + checkFluidSystem(); } + + // 2p-immiscible + { typedef Opm::FluidSystems::TwoPImmiscible FluidSystem; + checkFluidSystem(); } + + { typedef Opm::FluidSystems::TwoPImmiscible FluidSystem; + checkFluidSystem(); } + + { typedef Opm::FluidSystems::TwoPImmiscible FluidSystem; + checkFluidSystem(); } + + // 1p + { typedef Opm::FluidSystems::OneP FluidSystem; + checkFluidSystem(); } + + { typedef Opm::FluidSystems::OneP FluidSystem; + checkFluidSystem(); } + + return 0; +} diff --git a/tests/material/immiscibleflash/CMakeLists.txt b/tests/material/immiscibleflash/CMakeLists.txt new file mode 100644 index 000000000..2543a8706 --- /dev/null +++ b/tests/material/immiscibleflash/CMakeLists.txt @@ -0,0 +1,7 @@ +# set the CMake module include path +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules") + +# eWoms unit-testing infrastructure +include(EwomsAddTest) + +EwomsAddTest(test_immiscibleflash DRIVER_ARGS --plain test_immiscibleflash) diff --git a/tests/material/immiscibleflash/test_immiscibleflash.cc b/tests/material/immiscibleflash/test_immiscibleflash.cc new file mode 100644 index 000000000..ea3cc7b8c --- /dev/null +++ b/tests/material/immiscibleflash/test_immiscibleflash.cc @@ -0,0 +1,263 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2011-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief This is a program to test the flash calculation which uses + * non-linear complementarity problems (NCP) + * + * A flash calculation determines the pressures, saturations and + * composition of all phases given the total mass (or, as in this case + * the total number of moles) in a given amount of pore space. + */ +#include "config.h" + +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include + +template +void checkSame(const FluidState &fsRef, const FluidState &fsFlash) +{ + enum { numPhases = FluidState::numPhases }; + enum { numComponents = FluidState::numComponents }; + + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Scalar error; + + // check the pressures + error = 1 - fsRef.pressure(phaseIdx)/fsFlash.pressure(phaseIdx); + if (std::abs(error) > 1e-6) { + std::cout << "pressure error phase " << phaseIdx << ": " + << fsFlash.pressure(phaseIdx) << " flash vs " + << fsRef.pressure(phaseIdx) << " reference" + << " error=" << error << "\n"; + } + + // check the saturations + error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx); + if (std::abs(error) > 1e-6) + std::cout << "saturation error phase " << phaseIdx << ": " + << fsFlash.saturation(phaseIdx) << " flash vs " + << fsRef.saturation(phaseIdx) << " reference" + << " error=" << error << "\n"; + + // check the compositions + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { + error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx); + if (std::abs(error) > 1e-6) + std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": " + << fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs " + << fsRef.moleFraction(phaseIdx, compIdx) << " reference" + << " error=" << error << "\n"; + } + } +} + +template +void checkImmiscibleFlash(const FluidState &fsRef, + typename MaterialLaw::Params &matParams) +{ + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + typedef Dune::FieldVector ComponentVector; + + // calculate the total amount of stuff in the reference fluid + // phase + ComponentVector globalMolarities(0.0); + for (int compIdx = 0; compIdx < numComponents; ++compIdx) { + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + globalMolarities[compIdx] += + fsRef.saturation(phaseIdx)*fsRef.molarity(phaseIdx, compIdx); + } + } + + // initialize the fluid state for the flash calculation + typedef Opm::ImmiscibleFlash ImmiscibleFlash; + FluidState fsFlash; + + fsFlash.setTemperature(fsRef.temperature(/*phaseIdx=*/0)); + + // run the flash calculation + typename FluidSystem::ParameterCache paramCache; + ImmiscibleFlash::guessInitial(fsFlash, paramCache, globalMolarities); + ImmiscibleFlash::template solve(fsFlash, paramCache, matParams, globalMolarities); + + // compare the "flashed" fluid state with the reference one + checkSame(fsRef, fsFlash); +} + + +template +void completeReferenceFluidState(FluidState &fs, + typename MaterialLaw::Params &matParams, + int refPhaseIdx) +{ + enum { numPhases = FluidSystem::numPhases }; + typedef Dune::FieldVector PhaseVector; + + int otherPhaseIdx = 1 - refPhaseIdx; + + // calculate the other saturation + fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx)); + + // calulate the capillary pressure + PhaseVector pC; + MaterialLaw::capillaryPressures(pC, matParams, fs); + fs.setPressure(otherPhaseIdx, + fs.pressure(refPhaseIdx) + + (pC[otherPhaseIdx] - pC[refPhaseIdx])); + + // set all phase densities + typename FluidSystem::ParameterCache paramCache; + paramCache.updateAll(fs); + for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx); + fs.setDensity(phaseIdx, rho); + } +} + + +int main() +{ + typedef double Scalar; + typedef Opm::FluidSystems::H2ON2 FluidSystem; + typedef Opm::ImmiscibleFluidState ImmiscibleFluidState; + + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { lPhaseIdx = FluidSystem::lPhaseIdx }; + enum { gPhaseIdx = FluidSystem::gPhaseIdx }; + + enum { H2OIdx = FluidSystem::H2OIdx }; + enum { N2Idx = FluidSystem::N2Idx }; + + typedef Opm::RegularizedBrooksCorey EffMaterialLaw; + typedef Opm::EffToAbsLaw TwoPMaterialLaw; + typedef Opm::TwoPAdapter MaterialLaw; + typedef MaterialLaw::Params MaterialLawParams; + + Scalar T = 273.15 + 25; + + // initialize the tables of the fluid system + Scalar Tmin = T - 1.0; + Scalar Tmax = T + 1.0; + int nT = 3; + + Scalar pmin = 0.0; + Scalar pmax = 1.25 * 2e6; + int np = 100; + + FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np); + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + matParams.setSwr(0.0); + matParams.setSnr(0.0); + matParams.setPe(0); + matParams.setLambda(2.0); + + ImmiscibleFluidState fsRef; + + // create an fluid state which is consistent + + // set the fluid temperatures + fsRef.setTemperature(T); + + //////////////// + // only liquid + //////////////// + std::cout << "testing single-phase liquid\n"; + + // set liquid saturation and pressure + fsRef.setSaturation(lPhaseIdx, 1.0); + fsRef.setPressure(lPhaseIdx, 1e6); + + // set the remaining parameters of the reference fluid state + completeReferenceFluidState(fsRef, matParams, lPhaseIdx); + + // check the flash calculation + checkImmiscibleFlash(fsRef, matParams); + + //////////////// + // only gas + //////////////// + std::cout << "testing single-phase gas\n"; + + // set gas saturation and pressure + fsRef.setSaturation(gPhaseIdx, 1.0); + fsRef.setPressure(gPhaseIdx, 1e6); + + // set the remaining parameters of the reference fluid state + completeReferenceFluidState(fsRef, matParams, gPhaseIdx); + + // check the flash calculation + checkImmiscibleFlash(fsRef, matParams); + + //////////////// + // both phases + //////////////// + std::cout << "testing two-phase\n"; + + // set liquid saturation and pressure + fsRef.setSaturation(lPhaseIdx, 0.5); + fsRef.setPressure(lPhaseIdx, 1e6); + + // set the remaining parameters of the reference fluid state + completeReferenceFluidState(fsRef, matParams, lPhaseIdx); + + // check the flash calculation + checkImmiscibleFlash(fsRef, matParams); + + //////////////// + // with capillary pressure + //////////////// + std::cout << "testing two-phase with capillary pressure\n"; + + MaterialLawParams matParams2; + matParams2.setSwr(0.0); + matParams2.setSnr(0.0); + matParams2.setPe(1e3); + matParams2.setLambda(2.0); + + // set liquid saturation + fsRef.setSaturation(lPhaseIdx, 0.5); + + // set pressure of the liquid phase + fsRef.setPressure(lPhaseIdx, 1e6); + + // set the remaining parameters of the reference fluid state + completeReferenceFluidState(fsRef, matParams2, lPhaseIdx); + + // check the flash calculation + checkImmiscibleFlash(fsRef, matParams2); + + return 0; +} diff --git a/tests/material/ncpflash/CMakeLists.txt b/tests/material/ncpflash/CMakeLists.txt new file mode 100644 index 000000000..d668d668a --- /dev/null +++ b/tests/material/ncpflash/CMakeLists.txt @@ -0,0 +1,7 @@ +# set the CMake module include path +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules") + +# eWoms unit-testing infrastructure +include(EwomsAddTest) + +EwomsAddTest(test_ncpflash DRIVER_ARGS --plain test_ncpflash) diff --git a/tests/material/ncpflash/test_ncpflash.cc b/tests/material/ncpflash/test_ncpflash.cc new file mode 100644 index 000000000..21c6ee19a --- /dev/null +++ b/tests/material/ncpflash/test_ncpflash.cc @@ -0,0 +1,295 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2011-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief This is a program to test the flash calculation which uses + * non-linear complementarity problems (NCP) + * + * A flash calculation determines the pressures, saturations and + * composition of all phases given the total mass (or, as in this case + * the total number of moles) in a given amount of pore space. + */ +#include "config.h" + +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include + +template +void checkSame(const FluidState &fsRef, const FluidState &fsFlash) +{ + enum { numPhases = FluidState::numPhases }; + enum { numComponents = FluidState::numComponents }; + + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Scalar error; + + // check the pressures + error = 1 - fsRef.pressure(phaseIdx)/fsFlash.pressure(phaseIdx); + if (std::abs(error) > 1e-6) { + std::cout << "pressure error phase " << phaseIdx << ": " + << fsFlash.pressure(phaseIdx) << " flash vs " + << fsRef.pressure(phaseIdx) << " reference" + << " error=" << error << "\n"; + } + + // check the saturations + error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx); + if (std::abs(error) > 1e-6) + std::cout << "saturation error phase " << phaseIdx << ": " + << fsFlash.saturation(phaseIdx) << " flash vs " + << fsRef.saturation(phaseIdx) << " reference" + << " error=" << error << "\n"; + + // check the compositions + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { + error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx); + if (std::abs(error) > 1e-6) + std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": " + << fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs " + << fsRef.moleFraction(phaseIdx, compIdx) << " reference" + << " error=" << error << "\n"; + } + } +} + +template +void checkNcpFlash(const FluidState &fsRef, + typename MaterialLaw::Params &matParams) +{ + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + typedef Dune::FieldVector ComponentVector; + + // calculate the total amount of stuff in the reference fluid + // phase + ComponentVector globalMolarities(0.0); + for (int compIdx = 0; compIdx < numComponents; ++compIdx) { + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + globalMolarities[compIdx] += + fsRef.saturation(phaseIdx)*fsRef.molarity(phaseIdx, compIdx); + } + } + + // initialize the fluid state for the flash calculation + typedef Opm::NcpFlash NcpFlash; + FluidState fsFlash; + + fsFlash.setTemperature(fsRef.temperature(/*phaseIdx=*/0)); + + // run the flash calculation + typename FluidSystem::ParameterCache paramCache; + NcpFlash::guessInitial(fsFlash, paramCache, globalMolarities); + NcpFlash::template solve(fsFlash, paramCache, matParams, globalMolarities); + + // compare the "flashed" fluid state with the reference one + checkSame(fsRef, fsFlash); +} + + +template +void completeReferenceFluidState(FluidState &fs, + typename MaterialLaw::Params &matParams, + int refPhaseIdx) +{ + enum { numPhases = FluidSystem::numPhases }; + + typedef Opm::ComputeFromReferencePhase ComputeFromReferencePhase; + typedef Dune::FieldVector PhaseVector; + + int otherPhaseIdx = 1 - refPhaseIdx; + + // calculate the other saturation + fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx)); + + // calulate the capillary pressure + PhaseVector pC; + MaterialLaw::capillaryPressures(pC, matParams, fs); + fs.setPressure(otherPhaseIdx, + fs.pressure(refPhaseIdx) + + (pC[otherPhaseIdx] - pC[refPhaseIdx])); + + // make the fluid state consistent with local thermodynamic + // equilibrium + typename FluidSystem::ParameterCache paramCache; + ComputeFromReferencePhase::solve(fs, + paramCache, + refPhaseIdx, + /*setViscosity=*/false, + /*setEnthalpy=*/false); +} + + +int main() +{ + typedef double Scalar; + typedef Opm::FluidSystems::H2ON2 FluidSystem; + typedef Opm::CompositionalFluidState CompositionalFluidState; + + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { lPhaseIdx = FluidSystem::lPhaseIdx }; + enum { gPhaseIdx = FluidSystem::gPhaseIdx }; + + enum { H2OIdx = FluidSystem::H2OIdx }; + enum { N2Idx = FluidSystem::N2Idx }; + + typedef Opm::RegularizedBrooksCorey EffMaterialLaw; + typedef Opm::EffToAbsLaw TwoPMaterialLaw; + typedef Opm::TwoPAdapter MaterialLaw; + typedef MaterialLaw::Params MaterialLawParams; + + Scalar T = 273.15 + 25; + + // initialize the tables of the fluid system + Scalar Tmin = T - 1.0; + Scalar Tmax = T + 1.0; + int nT = 3; + + Scalar pmin = 0.0; + Scalar pmax = 1.25 * 2e6; + int np = 100; + + FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np); + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + matParams.setSwr(0.0); + matParams.setSnr(0.0); + matParams.setPe(0); + matParams.setLambda(2.0); + + CompositionalFluidState fsRef; + + // create an fluid state which is consistent + + // set the fluid temperatures + fsRef.setTemperature(T); + + //////////////// + // only liquid + //////////////// + std::cout << "testing single-phase liquid\n"; + + // set liquid saturation + fsRef.setSaturation(lPhaseIdx, 1.0); + + // set pressure of the liquid phase + fsRef.setPressure(lPhaseIdx, 2e5); + + // set the liquid composition to pure water + fsRef.setMoleFraction(lPhaseIdx, N2Idx, 0.0); + fsRef.setMoleFraction(lPhaseIdx, H2OIdx, 1.0 - fsRef.moleFraction(lPhaseIdx, N2Idx)); + + // "complete" the fluid state + completeReferenceFluidState(fsRef, matParams, lPhaseIdx); + + // check the flash calculation + checkNcpFlash(fsRef, matParams); + + //////////////// + // only gas + //////////////// + std::cout << "testing single-phase gas\n"; + // set gas saturation + fsRef.setSaturation(gPhaseIdx, 1.0); + + // set pressure of the gas phase + fsRef.setPressure(gPhaseIdx, 1e6); + + // set the gas composition to 99.9% nitrogen and 0.1% water + fsRef.setMoleFraction(gPhaseIdx, N2Idx, 0.999); + fsRef.setMoleFraction(gPhaseIdx, H2OIdx, 0.001); + + // "complete" the fluid state + completeReferenceFluidState(fsRef, matParams, gPhaseIdx); + + // check the flash calculation + checkNcpFlash(fsRef, matParams); + + //////////////// + // both phases + //////////////// + std::cout << "testing two-phase\n"; + + // set saturations + fsRef.setSaturation(lPhaseIdx, 0.5); + fsRef.setSaturation(gPhaseIdx, 0.5); + + // set pressures + fsRef.setPressure(lPhaseIdx, 1e6); + fsRef.setPressure(gPhaseIdx, 1e6); + + FluidSystem::ParameterCache paramCache; + typedef Opm::MiscibleMultiPhaseComposition MiscibleMultiPhaseComposition; + MiscibleMultiPhaseComposition::solve(fsRef, paramCache, + /*setViscosity=*/false, + /*setEnthalpy=*/false); + + // check the flash calculation + checkNcpFlash(fsRef, matParams); + + //////////////// + // with capillary pressure + //////////////// + + MaterialLawParams matParams2; + matParams2.setSwr(0.0); + matParams2.setSnr(0.0); + matParams2.setPe(1e3); + matParams2.setLambda(2.0); + + // set gas saturation + fsRef.setSaturation(gPhaseIdx, 0.5); + fsRef.setSaturation(lPhaseIdx, 0.5); + + // set pressure of the liquid phase + fsRef.setPressure(lPhaseIdx, 1e6); + + // calulate the capillary pressure + typedef Dune::FieldVector PhaseVector; + PhaseVector pC; + MaterialLaw::capillaryPressures(pC, matParams2, fsRef); + fsRef.setPressure(gPhaseIdx, + fsRef.pressure(lPhaseIdx) + + (pC[gPhaseIdx] - pC[lPhaseIdx])); + + typedef Opm::MiscibleMultiPhaseComposition MiscibleMultiPhaseComposition; + MiscibleMultiPhaseComposition::solve(fsRef, paramCache, + /*setViscosity=*/false, + /*setEnthalpy=*/false); + + + // check the flash calculation + checkNcpFlash(fsRef, matParams2); + + return 0; +} diff --git a/tests/material/pengrobinson/CMakeLists.txt b/tests/material/pengrobinson/CMakeLists.txt new file mode 100644 index 000000000..e41e3469e --- /dev/null +++ b/tests/material/pengrobinson/CMakeLists.txt @@ -0,0 +1,8 @@ +# set the CMake module include path +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules") + +# eWoms unit-testing infrastructure +include(EwomsAddTest) + +EwomsAddTest(test_pengrobinson DRIVER_ARGS --plain test_pengrobinson) + diff --git a/tests/material/pengrobinson/test_pengrobinson.cc b/tests/material/pengrobinson/test_pengrobinson.cc new file mode 100644 index 000000000..71f443b7c --- /dev/null +++ b/tests/material/pengrobinson/test_pengrobinson.cc @@ -0,0 +1,337 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief This is test for the SPE5 fluid system (which uses the + * Peng-Robinson EOS) and the NCP flash solver. + */ +#include "config.h" + +#include +#include +#include +#include +#include + +template +void guessInitial(FluidState &fluidState, + int phaseIdx) +{ + if (phaseIdx == FluidSystem::gPhaseIdx) { + fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.74785); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.0121364); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.00606028); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.00268136); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.000204256); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 8.78291e-06); + } + else if (phaseIdx == FluidSystem::oPhaseIdx) { + fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.50); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.03); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.07); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.20); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.15); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 0.05); + } + else { + assert(phaseIdx == FluidSystem::wPhaseIdx); + } +} + +template +Scalar bringOilToSurface(FluidState &surfaceFluidState, Scalar alpha, const FluidState &reservoirFluidState, bool guessInitial) +{ + enum { + numPhases = FluidSystem::numPhases, + wPhaseIdx = FluidSystem::wPhaseIdx, + gPhaseIdx = FluidSystem::gPhaseIdx, + oPhaseIdx = FluidSystem::oPhaseIdx, + + numComponents = FluidSystem::numComponents + }; + + typedef Opm::NcpFlash Flash; + typedef Opm::MpLinearMaterial MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + typedef Dune::FieldVector ComponentVector; + + const Scalar refPressure = 1.0135e5; // [Pa] + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + matParams.setPcMinSat(phaseIdx, 0.0); + matParams.setPcMaxSat(phaseIdx, 0.0); + } + + // retieve the global volumetric component molarities + surfaceFluidState.setTemperature(273.15 + 20); + + ComponentVector molarities; + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) + molarities[compIdx] = reservoirFluidState.molarity(oPhaseIdx, compIdx); + + if (guessInitial) { + // we start at a fluid state with reservoir oil. + for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { + surfaceFluidState.setMoleFraction(phaseIdx, + compIdx, + reservoirFluidState.moleFraction(phaseIdx, compIdx)); + } + surfaceFluidState.setDensity(phaseIdx, reservoirFluidState.density(phaseIdx)); + surfaceFluidState.setPressure(phaseIdx, reservoirFluidState.pressure(phaseIdx)); + surfaceFluidState.setSaturation(phaseIdx, 0.0); + } + surfaceFluidState.setSaturation(oPhaseIdx, 1.0); + } + + typename FluidSystem::ParameterCache paramCache; + paramCache.updateAll(surfaceFluidState); + + // increase volume until we are at surface pressure. use the + // newton method for this + ComponentVector tmpMolarities; + for (int i = 0;; ++i) { + if (i >= 20) + DUNE_THROW(Opm::NumericalProblem, + "Newton method did not converge after 20 iterations"); + + // calculate the deviation from the standard pressure + tmpMolarities = molarities; + tmpMolarities /= alpha; + Flash::template solve(surfaceFluidState, paramCache, matParams, tmpMolarities); + Scalar f = surfaceFluidState.pressure(gPhaseIdx) - refPressure; + + // calculate the derivative of the deviation from the standard + // pressure + Scalar eps = alpha*1e-10; + tmpMolarities = molarities; + tmpMolarities /= alpha + eps; + Flash::template solve(surfaceFluidState, paramCache, matParams, tmpMolarities); + Scalar fStar = surfaceFluidState.pressure(gPhaseIdx) - refPressure; + Scalar fPrime = (fStar - f)/eps; + + // newton update + Scalar delta = f/fPrime; + alpha -= delta; + if (std::abs(delta) < std::abs(alpha)*1e-9) { + break; + } + } + + // calculate the final result + tmpMolarities = molarities; + tmpMolarities /= alpha; + Flash::template solve(surfaceFluidState, paramCache, matParams, tmpMolarities); + return alpha; +} + +int main(int argc, char** argv) +{ + typedef double Scalar; + typedef Opm::FluidSystems::Spe5 FluidSystem; + + enum { + numPhases = FluidSystem::numPhases, + wPhaseIdx = FluidSystem::wPhaseIdx, + gPhaseIdx = FluidSystem::gPhaseIdx, + oPhaseIdx = FluidSystem::oPhaseIdx, + + numComponents = FluidSystem::numComponents, + H2OIdx = FluidSystem::H2OIdx, + C1Idx = FluidSystem::C1Idx, + C3Idx = FluidSystem::C3Idx, + C6Idx = FluidSystem::C6Idx, + C10Idx = FluidSystem::C10Idx, + C15Idx = FluidSystem::C15Idx, + C20Idx = FluidSystem::C20Idx + }; + + typedef Opm::NcpFlash Flash; + typedef Dune::FieldVector ComponentVector; + typedef Opm::CompositionalFluidState FluidState; + + typedef Opm::MpLinearMaterial MaterialLaw; + typedef MaterialLaw::Params MaterialLawParams; + + typedef FluidSystem::ParameterCache ParameterCache; + + //////////// + // Initialize the fluid system and create the capillary pressure + // parameters + //////////// + Scalar T = 273.15 + 20; // 20 deg Celsius + FluidSystem::init(/*minTemperature=*/T - 1, + /*maxTemperature=*/T + 1, + /*minPressure=*/1.0e4, + /*maxTemperature=*/40.0e6); + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + matParams.setPcMinSat(phaseIdx, 0.0); + matParams.setPcMaxSat(phaseIdx, 0.0); + } + + //////////// + // Create a fluid state + //////////// + FluidState fluidState; + ParameterCache paramCache; + + // temperature + fluidState.setTemperature(T); + + // oil pressure + fluidState.setPressure(oPhaseIdx, 4000 * 6894.7573); // 4000 PSI + + // oil saturation + fluidState.setSaturation(oPhaseIdx, 1.0); + + // composition: SPE-5 reservoir oil + fluidState.setMoleFraction(oPhaseIdx, H2OIdx, 0.0); + fluidState.setMoleFraction(oPhaseIdx, C1Idx, 0.50); + fluidState.setMoleFraction(oPhaseIdx, C3Idx, 0.03); + fluidState.setMoleFraction(oPhaseIdx, C6Idx, 0.07); + fluidState.setMoleFraction(oPhaseIdx, C10Idx, 0.20); + fluidState.setMoleFraction(oPhaseIdx, C15Idx, 0.15); + fluidState.setMoleFraction(oPhaseIdx, C20Idx, 0.05); + + // set the fugacities in the oil phase + paramCache.updatePhase(fluidState, oPhaseIdx); + ComponentVector fugVec; + for (int compIdx = 0; compIdx < numComponents; ++compIdx) { + Scalar Phi = + FluidSystem::fugacityCoefficient(fluidState, paramCache, oPhaseIdx, compIdx); + fluidState.setFugacityCoefficient(oPhaseIdx, compIdx, Phi); + fugVec[compIdx] = fluidState.fugacity(oPhaseIdx, compIdx); + + fluidState.setMoleFraction(gPhaseIdx, compIdx, fluidState.moleFraction(oPhaseIdx, compIdx)); + } + +/* + fluidState.setPressure(gPhaseIdx, fluidState.pressure(oPhaseIdx)); // 4000 PSI + { + Scalar TMin = fluidState.temperature(oPhaseIdx)/2; + Scalar TMax = fluidState.temperature(oPhaseIdx)*2; + int n = 1000; + for (int i = 0; i < n; ++i) { + Scalar T = TMin + (TMax - TMin)*i/(n - 1); + + fluidState.setTemperature(T); + paramCache.updatePhase(fluidState, oPhaseIdx); + paramCache.updatePhase(fluidState, gPhaseIdx); + + Scalar rhoO = FluidSystem::density(fluidState, paramCache, oPhaseIdx); + Scalar rhoG = FluidSystem::density(fluidState, paramCache, gPhaseIdx); + fluidState.setDensity(oPhaseIdx, rhoO); + fluidState.setDensity(gPhaseIdx, rhoG); + std::cout << T << " " + << fluidState.density(oPhaseIdx) << " " + << fluidState.density(gPhaseIdx) << " " + << paramCache.gasPhaseParams().a() << " " + << paramCache.gasPhaseParams().b() << " " + << "\n"; + } + exit(0); + } +*/ + + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (phaseIdx != oPhaseIdx) { + fluidState.setSaturation(phaseIdx, 0.0); + fluidState.setPressure(phaseIdx, fluidState.pressure(oPhaseIdx)); + } + + // initial guess for the composition + guessInitial(fluidState, phaseIdx); + } + + typedef Opm::ComputeFromReferencePhase CFRP; + CFRP::solve(fluidState, + paramCache, + /*refPhaseIdx=*/oPhaseIdx, + /*setViscosity=*/false, + /*setEnthalpy=*/false); + + //////////// + // Calculate the total molarities of the components + //////////// + ComponentVector molarities; + for (int compIdx = 0; compIdx < numComponents; ++ compIdx) + molarities[compIdx] = fluidState.saturation(oPhaseIdx)*fluidState.molarity(oPhaseIdx, compIdx); + + //////////// + // Gradually increase the volume for and calculate the gas + // formation factor, oil formation volume factor and gas formation + // volume factor. + //////////// + + FluidState flashFluidState, surfaceFluidState; + flashFluidState.assign(fluidState); + //Flash::guessInitial(flashFluidState, paramCache, molarities); + Flash::solve(flashFluidState, paramCache, matParams, molarities); + + Scalar surfaceAlpha = 1; + surfaceAlpha = bringOilToSurface(surfaceFluidState, surfaceAlpha, flashFluidState, /*guessInitial=*/true); + Scalar rho_gRef = surfaceFluidState.density(gPhaseIdx); + Scalar rho_oRef = surfaceFluidState.density(oPhaseIdx); + + std::cout << "alpha[-] p[Pa] S_g[-] rho_o[kg/m^3] rho_g[kg/m^3] [kg/mol] [kg/mol] R_s[m^3/m^3] B_g[-] B_o[-]\n"; + int n = 1000; + for (int i = 0; i < n; ++i) { + Scalar minAlpha = 0.98; + Scalar maxAlpha = 5.0; + + // ratio between the original and the current volume + Scalar alpha = minAlpha + (maxAlpha - minAlpha)*i/(n - 1); + + // increasing the volume means decreasing the molartity + ComponentVector curMolarities = molarities; + curMolarities /= alpha; + + // "flash" the modified reservoir oil + Flash::solve(flashFluidState, paramCache, matParams, curMolarities); + + surfaceAlpha = bringOilToSurface(surfaceFluidState, + surfaceAlpha, + flashFluidState, + /*guessInitial=*/false); + std::cout << alpha << " " + << flashFluidState.pressure(oPhaseIdx) << " " + << flashFluidState.saturation(gPhaseIdx) << " " + << flashFluidState.density(oPhaseIdx) << " " + << flashFluidState.density(gPhaseIdx) << " " + << flashFluidState.averageMolarMass(oPhaseIdx) << " " + << flashFluidState.averageMolarMass(gPhaseIdx) << " " + << surfaceFluidState.saturation(gPhaseIdx)*surfaceAlpha << " " + << rho_gRef/flashFluidState.density(gPhaseIdx) << " " + << rho_oRef/flashFluidState.density(oPhaseIdx) << " " + << "\n"; + } + + std::cout << "reference density oil [kg/m^3]: " << rho_oRef << "\n"; + std::cout << "reference density gas [kg/m^3]: " << rho_gRef << "\n"; + + return 0; +} diff --git a/tests/material/tabulation/CMakeLists.txt b/tests/material/tabulation/CMakeLists.txt new file mode 100644 index 000000000..d2c2aac42 --- /dev/null +++ b/tests/material/tabulation/CMakeLists.txt @@ -0,0 +1,8 @@ +# set the CMake module include path +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules") + +# eWoms unit-testing infrastructure +include(EwomsAddTest) + +EwomsAddTest(test_tabulation DRIVER_ARGS --plain test_tabulation) + diff --git a/tests/material/tabulation/test_tabulation.cc b/tests/material/tabulation/test_tabulation.cc new file mode 100644 index 000000000..f18495590 --- /dev/null +++ b/tests/material/tabulation/test_tabulation.cc @@ -0,0 +1,113 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief This is a program to test the tabulation class for of + * individual components. + * + * It either prints "success" or "error", it does not do anything + * else. + */ +#include "config.h" + +#include +#include + +bool success; + +template +void isSame(const char *str, Scalar v, Scalar vRef, Scalar tol=1e-3) +{ + if (std::abs( (v - vRef)/vRef ) > tol) { + std::cout << "error for \"" << str << "\": " << (v - vRef)/vRef*100 << "% difference (tolerance: " << tol*100 << "%)\n"; + success = false; + //exit(1); + } +} + +int main() +{ + typedef double Scalar; + typedef Opm::H2O IapwsH2O; + typedef Opm::TabulatedComponent TabulatedH2O; + + Scalar tempMin = 274.15; + Scalar tempMax = 622.15; + int nTemp = static_cast(tempMax - tempMin) * 6/8; + + Scalar pMin = 10.00; + Scalar pMax = IapwsH2O::vaporPressure(tempMax*1.1); + int nPress = 200; + + std::cout << "Creating tabulation with " << nTemp*nPress << " entries per quantity\n"; + TabulatedH2O::init(tempMin, tempMax, nTemp, + pMin, pMax, nPress); + + std::cout << "Checking tabulation\n"; + success = true; + int m = nTemp*3; + int n = nPress*3; + for (int i = 0; i < m; ++i) { + Scalar T = tempMin + (tempMax - tempMin)*Scalar(i)/m; + + if (i % std::max(1, m/1000) == 0) { + std::cout << Scalar(i)/m*100 << "% done \r"; + std::cout.flush(); + } + + isSame("vaporPressure", + TabulatedH2O::vaporPressure(T), + IapwsH2O::vaporPressure(T), + 1e-3); + for (int j = 0; j < n; ++j) { + Scalar p = pMin + (pMax - pMin)*Scalar(j)/n; + if (p < IapwsH2O::vaporPressure(T) * 1.001) { + Scalar tol = 1e-3; + if (p > IapwsH2O::vaporPressure(T)) + tol = 1e-2; + Scalar rho = IapwsH2O::gasDensity(T,p); + //isSame("Iapws::gasPressure", IapwsH2O::gasPressure(T,rho), p, 1e-6); + //isSame("gasPressure", TabulatedH2O::gasPressure(T,rho), p, 2e-2); + isSame("gasEnthalpy", TabulatedH2O::gasEnthalpy(T,p), IapwsH2O::gasEnthalpy(T,p), tol); + isSame("gasInternalEnergy", TabulatedH2O::gasInternalEnergy(T,p), IapwsH2O::gasInternalEnergy(T,p), tol); + isSame("gasDensity", TabulatedH2O::gasDensity(T,p), rho, tol); + isSame("gasViscosity", TabulatedH2O::gasViscosity(T,p), IapwsH2O::gasViscosity(T,p), tol); + } + + if (p > IapwsH2O::vaporPressure(T) / 1.001) { + Scalar tol = 1e-3; + if (p < IapwsH2O::vaporPressure(T)) + tol = 1e-2; + Scalar rho = IapwsH2O::liquidDensity(T,p); + //isSame("Iapws::liquidPressure", IapwsH2O::liquidPressure(T,rho), p, 1e-6); + //isSame("liquidPressure", TabulatedH2O::liquidPressure(T,rho), p, 2e-2); + isSame("liquidEnthalpy", TabulatedH2O::liquidEnthalpy(T,p), IapwsH2O::liquidEnthalpy(T,p), tol); + isSame("liquidInternalEnergy", TabulatedH2O::liquidInternalEnergy(T,p), IapwsH2O::liquidInternalEnergy(T,p), tol); + isSame("liquidDensity", TabulatedH2O::liquidDensity(T,p), rho, tol); + isSame("liquidViscosity", TabulatedH2O::liquidViscosity(T,p), IapwsH2O::liquidViscosity(T,p), tol); + } + } + //std::cerr << "\n"; + } + + if (success) + std::cout << "\nsuccess\n"; + return 0; +}