From 7bd162270957acab6bcecaa014df77fa871a0d8a Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 20 Sep 2013 15:04:27 +0200 Subject: [PATCH] use OPM instead of dune infrastructure wherever possible basically the only Dune thing which is still used are the FieldVector and FieldMatrix classes used by some constraint solvers. Until something similar goes into opm-core, opm-material must depend on dune-common... --- opm/common/exceptions.hh | 77 -- opm/common/fixedlengthspline_.hh | 412 ----------- opm/common/math.hh | 392 ---------- opm/common/spline.hh | 532 -------------- opm/common/splinecommon_.hh | 685 ------------------ opm/common/variablelengthspline_.hh | 479 ------------ .../binarycoefficients/air_mesitylene.hh | 7 +- opm/material/binarycoefficients/air_xylene.hh | 7 +- .../binarycoefficients/fullermethod.hh | 13 +- opm/material/binarycoefficients/h2o_co2.hh | 5 +- opm/material/components/air.hh | 8 +- opm/material/components/co2.hh | 11 +- opm/material/components/component.hh | 52 +- opm/material/components/h2o.hh | 33 +- opm/material/components/mesitylene.hh | 4 +- opm/material/components/simpleh2o.hh | 2 +- opm/material/components/tabulatedcomponent.hh | 37 +- opm/material/components/xylene.hh | 8 +- .../compositionfromfugacities.hh | 17 +- .../computefromreferencephase.hh | 5 +- .../constraintsolvers/immiscibleflash.hh | 13 +- .../misciblemultiphasecomposition.hh | 9 +- opm/material/constraintsolvers/ncpflash.hh | 19 +- .../dynamictabulated2dfunction.hh | 5 +- opm/material/eos/pengrobinson.hh | 13 +- opm/material/eos/pengrobinsonparams.hh | 2 +- .../2p/brookscoreyparams.hh | 2 +- .../2p/regularizedbrookscorey.hh | 2 +- .../2p/regularizedlinearmaterial.hh | 2 +- .../2p/regularizedvangenuchten.hh | 2 +- .../3p/3pparkervangenuchten.hh | 8 +- .../3p/3pparkervangenuchtenparams.hh | 4 +- .../fluidmatrixinteractions/mp/3padapter.hh | 12 +- .../mp/mplinearmaterial.hh | 7 +- .../mp/nullmateriallaw.hh | 5 +- .../fluidstates/compositionalfluidstate.hh | 2 +- .../fluidstatecompositionmodules.hh | 15 +- .../fluidstates/fluidstatedensitymodules.hh | 9 +- .../fluidstates/fluidstateenthalpymodules.hh | 11 +- .../fluidstates/fluidstatefugacitymodules.hh | 11 +- .../fluidstates/fluidstatepressuremodules.hh | 9 +- .../fluidstatesaturationmodules.hh | 9 +- .../fluidstatetemperaturemodules.hh | 9 +- .../fluidstates/fluidstateviscositymodules.hh | 9 +- .../fluidstates/immisciblefluidstate.hh | 2 +- opm/material/fluidstates/modularfluidstate.hh | 2 +- .../fluidstates/nonequilibriumfluidstate.hh | 2 +- .../fluidstates/pressureoverlayfluidstate.hh | 2 +- .../saturationoverlayfluidstate.hh | 2 +- .../temperatureoverlayfluidstate.hh | 2 +- opm/material/fluidsystems/1pfluidsystem.hh | 2 - .../fluidsystems/2pimmisciblefluidsystem.hh | 3 - opm/material/fluidsystems/basefluidsystem.hh | 49 +- .../fluidsystems/blackoilfluidsystem.hh | 19 +- .../fluidsystems/brineco2fluidsystem.hh | 4 +- .../fluidsystems/h2oairfluidsystem.hh | 15 +- .../h2oairmesitylenefluidsystem.hh | 14 +- .../fluidsystems/h2oairxylenefluidsystem.hh | 12 +- opm/material/fluidsystems/h2on2fluidsystem.hh | 5 +- .../h2on2liquidphasefluidsystem.hh | 5 +- opm/material/fluidsystems/spe5fluidsystem.hh | 4 +- .../fluidsystems/spe5parametercache.hh | 8 +- .../heatconduction/dummyheatconductionlaw.hh | 11 +- .../heatconduction/fluidconduction.hh | 2 +- opm/material/heatconduction/somerton.hh | 2 +- .../statictabulated2dfunction.hh | 7 +- .../tabulated2dfunction.hh | 5 +- opm/{common => material}/valgrind.hh | 14 +- .../material/fluidsystems/checkfluidsystem.hh | 45 +- .../fluidsystems/test_fluidsystems.cc | 12 +- .../pengrobinson/test_pengrobinson.cc | 2 +- 71 files changed, 326 insertions(+), 2901 deletions(-) delete mode 100644 opm/common/exceptions.hh delete mode 100644 opm/common/fixedlengthspline_.hh delete mode 100644 opm/common/math.hh delete mode 100644 opm/common/spline.hh delete mode 100644 opm/common/splinecommon_.hh delete mode 100644 opm/common/variablelengthspline_.hh rename opm/{common => material}/dynamictabulated2dfunction.hh (97%) rename opm/{common => material}/statictabulated2dfunction.hh (95%) rename opm/{common => material}/tabulated2dfunction.hh (97%) rename opm/{common => material}/valgrind.hh (94%) diff --git a/opm/common/exceptions.hh b/opm/common/exceptions.hh deleted file mode 100644 index 3004455fb..000000000 --- a/opm/common/exceptions.hh +++ /dev/null @@ -1,77 +0,0 @@ -// -*- 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 deleted file mode 100644 index 8e090580c..000000000 --- a/opm/common/fixedlengthspline_.hh +++ /dev/null @@ -1,412 +0,0 @@ -// -*- 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 deleted file mode 100644 index 55b2376d2..000000000 --- a/opm/common/math.hh +++ /dev/null @@ -1,392 +0,0 @@ -// -*- 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 deleted file mode 100644 index 60f6d25ed..000000000 --- a/opm/common/spline.hh +++ /dev/null @@ -1,532 +0,0 @@ -// -*- 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 deleted file mode 100644 index 11f79e15b..000000000 --- a/opm/common/splinecommon_.hh +++ /dev/null @@ -1,685 +0,0 @@ -// -*- 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/variablelengthspline_.hh b/opm/common/variablelengthspline_.hh deleted file mode 100644 index 6e743e253..000000000 --- a/opm/common/variablelengthspline_.hh +++ /dev/null @@ -1,479 +0,0 @@ -// -*- 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 db9c166a9..73df039a2 100644 --- a/opm/material/binarycoefficients/air_mesitylene.hh +++ b/opm/material/binarycoefficients/air_mesitylene.hh @@ -43,9 +43,7 @@ public: */ template static Scalar henry(Scalar temperature) - { DUNE_THROW(Dune::NotImplemented, - "Henry coefficient of air in mesitylene"); - } + { OPM_THROW(std::runtime_error, "Not implemented: Henry coefficient of air in mesitylene"); } /*! * \brief Binary diffusion coefficent [m^2/s] for air and mesitylene. @@ -96,8 +94,7 @@ public: */ template static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, - "Binary liquid diffusion coefficients of air and mesitylene"); + { OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and mesitylene"); } }; diff --git a/opm/material/binarycoefficients/air_xylene.hh b/opm/material/binarycoefficients/air_xylene.hh index 82a19c7e4..f1e8ecce6 100644 --- a/opm/material/binarycoefficients/air_xylene.hh +++ b/opm/material/binarycoefficients/air_xylene.hh @@ -43,9 +43,7 @@ public: */ template static Scalar henry(Scalar temperature) - { DUNE_THROW(Dune::NotImplemented, - "Henry coefficient of air in xylene"); - } + { OPM_THROW(std::runtime_error, "Not implemented: Henry coefficient of air in xylene"); } /*! * \brief Binary diffusion coefficent [m^2/s] for air and xylene. @@ -96,8 +94,7 @@ public: */ template static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, - "Binary liquid diffusion coefficients of air and xylene"); + { OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and xylene"); } }; diff --git a/opm/material/binarycoefficients/fullermethod.hh b/opm/material/binarycoefficients/fullermethod.hh index 958fb1b6f..3086c7df0 100644 --- a/opm/material/binarycoefficients/fullermethod.hh +++ b/opm/material/binarycoefficients/fullermethod.hh @@ -23,12 +23,13 @@ #ifndef OPM_FULLERMETHOD_HH #define OPM_FULLERMETHOD_HH -#include +#include + +#include + +namespace Opm { +namespace BinaryCoeff { -namespace Opm -{ -namespace BinaryCoeff -{ /*! * \ingroup Binarycoefficients * \brief Estimate binary diffusion coefficents \f$\mathrm{[m^2/s]}\f$ in gases according to @@ -53,7 +54,7 @@ inline Scalar fullerMethod(const Scalar *M, // molar masses [g/mol] const Scalar pressure) // [Pa] { // "effective" molar mass in [g/m^3] - Scalar Mab = harmonicMean(M[0], M[1]); + Scalar Mab = Opm::utils::harmonicAverage(M[0], M[1]); // Fuller's method Scalar tmp = std::pow(SigmaNu[0], 1./3) + std::pow(SigmaNu[1], 1./3); diff --git a/opm/material/binarycoefficients/h2o_co2.hh b/opm/material/binarycoefficients/h2o_co2.hh index fc0bde485..8c8113bd9 100644 --- a/opm/material/binarycoefficients/h2o_co2.hh +++ b/opm/material/binarycoefficients/h2o_co2.hh @@ -87,10 +87,7 @@ public: */ template static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure) - { - DUNE_THROW(Dune::NotImplemented, - "Binary liquid diffusion coefficients of CO2 and CH4"); - } + { OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of CO2 and CH4"); } }; } diff --git a/opm/material/components/air.hh b/opm/material/components/air.hh index 4705fa0e9..6c6333adb 100644 --- a/opm/material/components/air.hh +++ b/opm/material/components/air.hh @@ -24,10 +24,12 @@ #ifndef OPM_AIR_HH #define OPM_AIR_HH -#include #include #include +#include +#include + namespace Opm { /*! @@ -156,8 +158,8 @@ public: Scalar r; if(temperature < 273.15 || temperature > 660.) { - DUNE_THROW(Dune::NotImplemented, - "ConstrelAir: Temperature out of range at "); + OPM_THROW(MaterialLawProblem, + "ConstrelAir: Temperature out of range at "); } r = 1.496*1.E-6*std::pow(temperature,1.5)/(temperature+120.); return (r); diff --git a/opm/material/components/co2.hh b/opm/material/components/co2.hh index 9fb4cbbad..d1a9ea849 100644 --- a/opm/material/components/co2.hh +++ b/opm/material/components/co2.hh @@ -24,7 +24,8 @@ #ifndef OPM_CO2_HH #define OPM_CO2_HH -#include +#include +#include #include #include #include @@ -163,8 +164,8 @@ public: #ifndef NDEBUG if ((temperature < criticalTemperature() or pressure < criticalPressure()) and !warningPrinted) { - Dune::dwarn << "The tables used for the CO2 exhibit subcritical as well as critical values: " - << "Double-check your results if you cross the critical-subcritical line during the simulation!\n"; + std::cout << "The tables used for the CO2 exhibit subcritical as well as critical values: " + << "Double-check your results if you cross the critical-subcritical line during the simulation!\n"; warningPrinted=true; } #endif @@ -192,8 +193,8 @@ public: #ifndef NDEBUG if ((temperature < criticalTemperature() or pressure < criticalPressure()) and !warningPrinted) { - Dune::dwarn << "The tables used for the CO2 exhibit subcritical as well as critical values: " - << "Double-check your results if you cross the critical-subcritical line during the simulation!\n"; + std::cout << "The tables used for the CO2 exhibit subcritical as well as critical values: " + << "Double-check your results if you cross the critical-subcritical line during the simulation!\n"; warningPrinted=true; } #endif diff --git a/opm/material/components/component.hh b/opm/material/components/component.hh index cc64d274b..d1578e44e 100644 --- a/opm/material/components/component.hh +++ b/opm/material/components/component.hh @@ -25,10 +25,10 @@ #ifndef OPM_COMPONENT_HH #define OPM_COMPONENT_HH -#include +#include +#include -namespace Opm -{ +namespace Opm { /*! * \ingroup Components @@ -57,61 +57,61 @@ public: */ static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp, Scalar pressMin, Scalar pressMax, unsigned nPress) - { Dune::dwarn << "No init routine defined - make sure that this is not necessary!" << std::endl; } + { } /*! * \brief Returns true iff the gas phase is assumed to be compressible */ static bool gasIsCompressible() - { DUNE_THROW(Dune::NotImplemented, "Component::gasIsCompressible()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasIsCompressible()"); } /*! * \brief Returns true iff the gas phase is assumed to be ideal */ static bool gasIsIdeal() - { DUNE_THROW(Dune::NotImplemented, "Component::gasIsIdeal()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasIsIdeal()"); } /*! * \brief Returns true iff the liquid phase is assumed to be compressible */ static bool liquidIsCompressible() - { DUNE_THROW(Dune::NotImplemented, "Component::liquidIsCompressible()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidIsCompressible()"); } /*! * \brief A human readable name for the component. */ static const char *name() - { DUNE_THROW(Dune::NotImplemented, "Component::name()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::name()"); } /*! * \brief The molar mass in \f$\mathrm{[kg]}\f$ of the component. */ static Scalar molarMass() - { DUNE_THROW(Dune::NotImplemented, "Component::molarMass()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::molarMass()"); } /*! * \brief Returns the critical temperature in \f$\mathrm{[K]}\f$ of the component. */ static Scalar criticalTemperature() - { DUNE_THROW(Dune::NotImplemented, "Component::criticalTemperature()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::criticalTemperature()"); } /*! * \brief Returns the critical pressure in \f$\mathrm{[Pa]}\f$ of the component. */ static Scalar criticalPressure() - { DUNE_THROW(Dune::NotImplemented, "Component::criticalPressure()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::criticalPressure()"); } /*! * \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point. */ static Scalar tripleTemperature() - { DUNE_THROW(Dune::NotImplemented, "Component::tripleTemperature()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::tripleTemperature()"); } /*! * \brief Returns the pressure in \f$\mathrm{[Pa]}\f$ at the component's triple point. */ static Scalar triplePressure() - { DUNE_THROW(Dune::NotImplemented, "Component::triplePressure()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::triplePressure()"); } /*! * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given @@ -120,7 +120,7 @@ public: * \param T temperature of the component in \f$\mathrm{[K]}\f$ */ static Scalar vaporPressure(Scalar T) - { DUNE_THROW(Dune::NotImplemented, "Component::vaporPressure()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::vaporPressure()"); } /*! * \brief The density in \f$\mathrm{[kg/m^3]}\f$ of the component at a given pressure in \f$\mathrm{[Pa]}\f$ and temperature in \f$\mathrm{[K]}\f$. @@ -129,7 +129,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static Scalar gasDensity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasDensity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasDensity()"); } /*! * \brief The density \f$\mathrm{[kg/m^3]}\f$ of the liquid component at a given pressure in \f$\mathrm{[Pa]}\f$ and temperature in \f$\mathrm{[K]}\f$. @@ -138,7 +138,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static Scalar liquidDensity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidDensity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidDensity()"); } /*! * \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of the pure component in gas. @@ -147,7 +147,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasEnthalpy()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasEnthalpy()"); } /*! * \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of the pure component in liquid. @@ -156,7 +156,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidEnthalpy()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidEnthalpy()"); } /*! * \brief Specific internal energy \f$\mathrm{[J/kg]}\f$ of the pure component in gas. @@ -165,7 +165,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasInternalEnergy()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasInternalEnergy()"); } /*! * \brief Specific internal energy \f$\mathrm{[J/kg]}\f$ of pure the pure component in liquid. @@ -174,7 +174,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidInternalEnergy()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidInternalEnergy()"); } /*! * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of the pure component at a given pressure in \f$\mathrm{[Pa]}\f$ and @@ -184,7 +184,7 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static Scalar gasViscosity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasViscosity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasViscosity()"); } /*! * \brief The dynamic liquid viscosity \f$\mathrm{[Pa*s]}\f$ of the pure component. @@ -193,31 +193,31 @@ public: * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ static Scalar liquidViscosity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidViscosity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidViscosity()"); } /*! * \brief Thermal conductivity of the component [W/(m^2 K/m)] as a gas. */ static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasThermalConductivity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasThermalConductivity()"); } /*! * \brief Thermal conductivity of the component [W/(m^2 K/m)] as a liquid. */ static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidThermalConductivity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidThermalConductivity()"); } /*! * \brief Specific isobaric heat capacity of the component [J/kg] as a gas. */ static Scalar gasHeatCapacity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasHeatCapacity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::gasHeatCapacity()"); } /*! * \brief Specific isobaric heat capacity of the component [J/kg] as a liquid. */ static Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidHeatCapacity()"); } + { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidHeatCapacity()"); } }; } // end namepace diff --git a/opm/material/components/h2o.hh b/opm/material/components/h2o.hh index 71a58c4fe..03b58b426 100644 --- a/opm/material/components/h2o.hh +++ b/opm/material/components/h2o.hh @@ -28,8 +28,9 @@ #include #include -#include -#include +#include +#include +#include #include "component.hh" @@ -172,7 +173,7 @@ public: { if (!Region2::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Enthalpy of steam is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -222,7 +223,7 @@ public: { if (!Region1::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Enthalpy of water is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -263,7 +264,7 @@ public: { if (!Region2::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Heat capacity of steam is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -299,7 +300,7 @@ public: { if (!Region1::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "heat Capacity of water is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -332,7 +333,7 @@ public: { if (!Region1::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Internal Energy of water is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -388,7 +389,7 @@ public: { if (!Region2::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Internal Energy of steam is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -462,7 +463,7 @@ public: { if (!Region1::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Heat capacity of water is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -495,7 +496,7 @@ public: { if (!Region2::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Heat capacity of steam is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -541,7 +542,7 @@ public: { if (!Region2::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Density of steam is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -662,7 +663,7 @@ public: { if (!Region1::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Density of water is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -761,7 +762,7 @@ public: { if (!Region2::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Viscosity of steam is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); } @@ -785,7 +786,7 @@ public: { if (!Region1::isValid(temperature, pressure)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Viscosity of water is only implemented for temperatures below 623.15K and " "pressures below 100MPa. (T = " << temperature << ", p=" << pressure); }; @@ -818,7 +819,7 @@ public: || (pressure <= 150e6 && ((523.15 < temperature) || (temperature > 673.15)) ) || (pressure <= 100e6 && ((673.15 < temperature) || (temperature > 1073.15)) ) ) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Liquid thermal conductivity of H2O for " << "T="< 673.15)) ) || (pressure <= 100e6 && ((673.15 < temperature) || (temperature > 1073.15)) ) ) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Gas thermal conductivity of H2O for " " T="< #include -#include +#include +#include namespace Opm { @@ -109,7 +110,7 @@ public: Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; try { vaporPressure_[iT] = RawComponent::vaporPressure(temperature); } - catch (Dune::NotImplemented) { vaporPressure_[iT] = NaN; } + catch (std::exception) { vaporPressure_[iT] = NaN; } catch (NumericalProblem e) { vaporPressure_[iT] = NaN; }; Scalar pgMax = maxGasPressure_(iT); @@ -122,23 +123,23 @@ public: unsigned i = iT + iP*nTemp_; try { gasEnthalpy_[i] = RawComponent::gasEnthalpy(temperature, pressure); } - catch (Dune::NotImplemented) { gasEnthalpy_[i] = NaN; } + catch (std::exception) { gasEnthalpy_[i] = NaN; } catch (NumericalProblem) { gasEnthalpy_[i] = NaN; }; try { gasHeatCapacity_[i] = RawComponent::gasHeatCapacity(temperature, pressure); } - catch (Dune::NotImplemented) { gasHeatCapacity_[i] = NaN; } + catch (std::exception) { gasHeatCapacity_[i] = NaN; } catch (NumericalProblem) { gasHeatCapacity_[i] = NaN; }; try { gasDensity_[i] = RawComponent::gasDensity(temperature, pressure); } - catch (Dune::NotImplemented) { gasDensity_[i] = NaN; } + catch (std::exception) { gasDensity_[i] = NaN; } catch (NumericalProblem) { gasDensity_[i] = NaN; }; try { gasViscosity_[i] = RawComponent::gasViscosity(temperature, pressure); } - catch (Dune::NotImplemented) { gasViscosity_[i] = NaN; } + catch (std::exception) { gasViscosity_[i] = NaN; } catch (NumericalProblem) { gasViscosity_[i] = NaN; }; try { gasThermalConductivity_[i] = RawComponent::gasThermalConductivity(temperature, pressure); } - catch (Dune::NotImplemented) { gasThermalConductivity_[i] = NaN; } + catch (std::exception) { gasThermalConductivity_[i] = NaN; } catch (NumericalProblem) { gasThermalConductivity_[i] = NaN; }; }; @@ -150,23 +151,23 @@ public: unsigned i = iT + iP*nTemp_; try { liquidEnthalpy_[i] = RawComponent::liquidEnthalpy(temperature, pressure); } - catch (Dune::NotImplemented) { liquidEnthalpy_[i] = NaN; } + catch (std::exception) { liquidEnthalpy_[i] = NaN; } catch (NumericalProblem) { liquidEnthalpy_[i] = NaN; }; try { liquidHeatCapacity_[i] = RawComponent::liquidHeatCapacity(temperature, pressure); } - catch (Dune::NotImplemented) { liquidHeatCapacity_[i] = NaN; } + catch (std::exception) { liquidHeatCapacity_[i] = NaN; } catch (NumericalProblem) { liquidHeatCapacity_[i] = NaN; }; try { liquidDensity_[i] = RawComponent::liquidDensity(temperature, pressure); } - catch (Dune::NotImplemented) { liquidDensity_[i] = NaN; } + catch (std::exception) { liquidDensity_[i] = NaN; } catch (NumericalProblem) { liquidDensity_[i] = NaN; }; try { liquidViscosity_[i] = RawComponent::liquidViscosity(temperature, pressure); } - catch (Dune::NotImplemented) { liquidViscosity_[i] = NaN; } + catch (std::exception) { liquidViscosity_[i] = NaN; } catch (NumericalProblem) { liquidViscosity_[i] = NaN; }; try { liquidThermalConductivity_[i] = RawComponent::liquidThermalConductivity(temperature, pressure); } - catch (Dune::NotImplemented) { liquidThermalConductivity_[i] = NaN; } + catch (std::exception) { liquidThermalConductivity_[i] = NaN; } catch (NumericalProblem) { liquidThermalConductivity_[i] = NaN; }; } } @@ -595,13 +596,13 @@ private: #if 0 && !defined NDEBUG if(!(0 <= alphaT && alphaT <= 1.0)) - DUNE_THROW(NumericalProblem, "Temperature out of range: " + OPM_THROW(NumericalProblem, "Temperature out of range: " << "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]"); if(!(0 <= alphaP1 && alphaP1 <= 1.0)) - DUNE_THROW(NumericalProblem, "First liquid pressure out of range: " + OPM_THROW(NumericalProblem, "First liquid pressure out of range: " << "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T)) << ", " << maxLiquidPressure_(tempIdx_(T)) << "]"); if(!(0 <= alphaP2 && alphaP2 <= 1.0)) - DUNE_THROW(NumericalProblem, "Second liquid pressure out of range: " + OPM_THROW(NumericalProblem, "Second liquid pressure out of range: " << "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T) + 1) << ", " << maxLiquidPressure_(tempIdx_(T) + 1) << "]"); #endif @@ -633,13 +634,13 @@ private: #if 0 && !defined NDEBUG if(!(0 <= alphaT && alphaT <= 1.0)) - DUNE_THROW(NumericalProblem, "Temperature out of range: " + OPM_THROW(NumericalProblem, "Temperature out of range: " << "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]"); if(!(0 <= alphaP1 && alphaP1 <= 1.0)) - DUNE_THROW(NumericalProblem, "First gas pressure out of range: " + OPM_THROW(NumericalProblem, "First gas pressure out of range: " << "p=" << p << " range: [" << minGasPressure_(tempIdx_(T)) << ", " << maxGasPressure_(tempIdx_(T)) << "]"); if(!(0 <= alphaP2 && alphaP2 <= 1.0)) - DUNE_THROW(NumericalProblem, "Second gas pressure out of range: " + OPM_THROW(NumericalProblem, "Second gas pressure out of range: " << "p=" << p << " range: [" << minGasPressure_(tempIdx_(T) + 1) << ", " << maxGasPressure_(tempIdx_(T) + 1) << "]"); #endif diff --git a/opm/material/components/xylene.hh b/opm/material/components/xylene.hh index d87c5cf25..42643ba8f 100644 --- a/opm/material/components/xylene.hh +++ b/opm/material/components/xylene.hh @@ -73,17 +73,13 @@ public: * \brief Returns the temperature \f$\mathrm{[K]}\f$ at xylene's triple point. */ static const Scalar tripleTemperature() - { - DUNE_THROW(Dune::NotImplemented, "tripleTemperature for xylene"); - } + { OPM_THROW(std::runtime_error, "Not implemented: tripleTemperature for xylene"); } /*! * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at xylene's triple point. */ static const Scalar triplePressure() - { - DUNE_THROW(Dune::NotImplemented, "triplePressure for xylene"); - } + { OPM_THROW(std::runtime_error, "Not implemented: triplePressure for xylene"); } /*! * \brief The saturation vapor pressure in \f$\mathrm{[Pa]}\f$ of pure xylene diff --git a/opm/material/constraintsolvers/compositionfromfugacities.hh b/opm/material/constraintsolvers/compositionfromfugacities.hh index 61bac0d8b..b83bb4de6 100644 --- a/opm/material/constraintsolvers/compositionfromfugacities.hh +++ b/opm/material/constraintsolvers/compositionfromfugacities.hh @@ -26,8 +26,9 @@ #include #include -#include -#include +#include +#include +#include #include @@ -165,12 +166,12 @@ public: } } - DUNE_THROW(NumericalProblem, - "Calculating the " << FluidSystem::phaseName(phaseIdx) - << "Phase composition failed. Initial {x} = {" - << xInit - << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx) - << ", T = " << fluidState.temperature(phaseIdx)); + OPM_THROW(Opm::NumericalProblem, + "Calculating the " << FluidSystem::phaseName(phaseIdx) + << "Phase composition failed. Initial {x} = {" + << xInit + << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx) + << ", T = " << fluidState.temperature(phaseIdx)); } diff --git a/opm/material/constraintsolvers/computefromreferencephase.hh b/opm/material/constraintsolvers/computefromreferencephase.hh index 8bb312902..97d3b58df 100644 --- a/opm/material/constraintsolvers/computefromreferencephase.hh +++ b/opm/material/constraintsolvers/computefromreferencephase.hh @@ -25,8 +25,9 @@ #include -#include -#include +#include +#include +#include #include diff --git a/opm/material/constraintsolvers/immiscibleflash.hh b/opm/material/constraintsolvers/immiscibleflash.hh index c2e3b2d3c..7f74f17f3 100644 --- a/opm/material/constraintsolvers/immiscibleflash.hh +++ b/opm/material/constraintsolvers/immiscibleflash.hh @@ -23,8 +23,9 @@ #ifndef OPM_IMMISCIBLE_FLASH_HH #define OPM_IMMISCIBLE_FLASH_HH -#include -#include +#include +#include +#include #include #include @@ -214,10 +215,10 @@ public: std::cout << "\n"; */ - DUNE_THROW(NumericalProblem, - "Flash calculation failed." - " {c_alpha^kappa} = {" << globalMolarities << "}, T = " - << fluidState.temperature(/*phaseIdx=*/0)); + OPM_THROW(Opm::NumericalProblem, + "Flash calculation failed." + " {c_alpha^kappa} = {" << globalMolarities << "}, T = " + << fluidState.temperature(/*phaseIdx=*/0)); } diff --git a/opm/material/constraintsolvers/misciblemultiphasecomposition.hh b/opm/material/constraintsolvers/misciblemultiphasecomposition.hh index 924256fe9..b44223707 100644 --- a/opm/material/constraintsolvers/misciblemultiphasecomposition.hh +++ b/opm/material/constraintsolvers/misciblemultiphasecomposition.hh @@ -26,8 +26,9 @@ #include #include -#include -#include +#include +#include +#include namespace Opm { @@ -242,8 +243,8 @@ public: M.solve(x, b); } catch (const Dune::FMatrixError &e) { - DUNE_THROW(NumericalProblem, - "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << NumericalProblem(e.what()) << "; M="< +#include namespace Opm { /*! diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh b/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh index e37b4b73e..7c31c0cad 100644 --- a/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh +++ b/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh @@ -28,7 +28,7 @@ #include "linearmaterial.hh" #include "regularizedlinearmaterialparams.hh" -#include +#include namespace Opm { /*! diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh index eb555e7a9..7a92f09d7 100644 --- a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh +++ b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh @@ -28,7 +28,7 @@ #include "vangenuchten.hh" #include "regularizedvangenuchtenparams.hh" -#include +#include #include diff --git a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh index 9239d7264..bace198e1 100644 --- a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh +++ b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh @@ -52,7 +52,7 @@ public: */ static Scalar pC(const Params ¶ms, Scalar Sw) { - DUNE_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pCGN, pCNW, and pcGW"); + OPM_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pCGN, pCNW, and pcGW"); } static Scalar pCGW(const Params ¶ms, Scalar Sw) @@ -208,7 +208,7 @@ public: */ static Scalar Sw(const Params ¶ms, Scalar pC) { - DUNE_THROW(Dune::NotImplemented, "Sw(pc) for three phases not implemented! Do it yourself!"); + OPM_THROW(Dune::NotImplemented, "Sw(pc) for three phases not implemented! Do it yourself!"); } /*! @@ -218,7 +218,7 @@ public: */ static Scalar dpC_dSw(const Params ¶ms, Scalar Sw) { - DUNE_THROW(Dune::NotImplemented, "dpC/dSw for three phases not implemented! Do it yourself!"); + OPM_THROW(Dune::NotImplemented, "dpC/dSw for three phases not implemented! Do it yourself!"); } /*! @@ -227,7 +227,7 @@ public: */ static Scalar dSw_dpC(const Params ¶ms, Scalar pC) { - DUNE_THROW(Dune::NotImplemented, "dSw/dpC for three phases not implemented! Do it yourself!"); + OPM_THROW(Dune::NotImplemented, "dSw/dpC for three phases not implemented! Do it yourself!"); } /*! diff --git a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh index a414e1d28..6eb0e8a59 100644 --- a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh +++ b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh @@ -27,7 +27,7 @@ #include -#include +#include namespace Opm { /*! @@ -128,7 +128,7 @@ public: return Sgr_; break; }; - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } /*! diff --git a/opm/material/fluidmatrixinteractions/mp/3padapter.hh b/opm/material/fluidmatrixinteractions/mp/3padapter.hh index 877b712a6..4f916cf05 100644 --- a/opm/material/fluidmatrixinteractions/mp/3padapter.hh +++ b/opm/material/fluidmatrixinteractions/mp/3padapter.hh @@ -20,12 +20,12 @@ * \file * \copydoc Opm::ThreePAdapter */ -#ifndef OPM_MP_3P_ADAPTER_HH -#define OPM_MP_3P_ADAPTER_HH +#ifndef OPM_3P_ADAPTER_HPP +#define OPM_3P_ADAPTER_HPP -#include - -#include +#include +#include +#include #include @@ -80,7 +80,7 @@ public: static void saturations(ContainerT &values, const Params ¶ms, const FluidState &fluidState) - { DUNE_THROW(Dune::NotImplemented, "Inverse capillary pressure curves"); } + { OPM_THROW(std::runtime_error, "Not implemented: Inverse capillary pressure curves"); } /*! * \brief The relative permeability of all phases. diff --git a/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh b/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh index 233796873..3834b44f3 100644 --- a/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh +++ b/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh @@ -25,9 +25,10 @@ #include "mplinearmaterialparams.hh" -#include +#include -#include +#include +#include #include @@ -93,7 +94,7 @@ public: const Params ¶ms, const FluidState &state) { - DUNE_THROW(Dune::NotImplemented, "MpLinearMaterial::saturations()"); + OPM_THROW(std::runtime_error, "Not implemented: MpLinearMaterial::saturations()"); } /*! diff --git a/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh b/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh index e7008cab7..fcbd6e5e9 100644 --- a/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh +++ b/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh @@ -25,7 +25,8 @@ #include "nullmateriallawparams.hh" -#include +#include +#include #include @@ -69,7 +70,7 @@ public: const Params ¶ms, const FluidState &state) { - DUNE_THROW(Dune::NotImplemented, "MpLinearMaterial::saturations()"); + OPM_THROW(std::runtime_error, "Not implemented: MpLinearMaterial::saturations()"); } /*! diff --git a/opm/material/fluidstates/compositionalfluidstate.hh b/opm/material/fluidstates/compositionalfluidstate.hh index ef5b8a582..cd3ab54e4 100644 --- a/opm/material/fluidstates/compositionalfluidstate.hh +++ b/opm/material/fluidstates/compositionalfluidstate.hh @@ -28,7 +28,7 @@ #include "modularfluidstate.hh" -#include +#include #include namespace Opm { diff --git a/opm/material/fluidstates/fluidstatecompositionmodules.hh b/opm/material/fluidstates/fluidstatecompositionmodules.hh index 0618ffd8a..dbd463b3e 100644 --- a/opm/material/fluidstates/fluidstatecompositionmodules.hh +++ b/opm/material/fluidstates/fluidstatecompositionmodules.hh @@ -24,9 +24,10 @@ #ifndef OPM_FLUID_STATE_COMPOSITION_MODULES_HH #define OPM_FLUID_STATE_COMPOSITION_MODULES_HH -#include +#include -#include +#include +#include #include #include @@ -238,7 +239,7 @@ protected: /*! * \brief Module for the modular fluid state which does not store the - * compositions but throws Dune::InvalidState instead. + * compositions but throws std::logic_error instead. */ template +#include -#include +#include +#include #include @@ -105,7 +106,7 @@ protected: /*! * \brief Module for the modular fluid state which does not the - * densitys but throws Dune::InvalidState instead. + * densitys but throws std::logic_error instead. */ template +#include -#include +#include +#include #include @@ -99,7 +100,7 @@ protected: /*! * \brief Module for the modular fluid state which does not store the - * enthalpies but throws Dune::InvalidState instead. + * enthalpies but throws std::logic_error instead. */ template +#include -#include +#include +#include #include #include @@ -173,7 +174,7 @@ protected: /*! * \brief Module for the modular fluid state which does not store the - * fugacitys but throws Dune::InvalidState instead. + * fugacitys but throws std::logic_error instead. */ template +#include -#include +#include +#include #include @@ -91,7 +92,7 @@ protected: /*! * \brief Module for the modular fluid state which does not the - * pressures but throws Dune::InvalidState instead. + * pressures but throws std::logic_error instead. */ template +#include -#include +#include +#include #include @@ -90,7 +91,7 @@ protected: /*! * \brief Module for the modular fluid state which does not the - * saturations but throws Dune::InvalidState instead. + * saturations but throws std::logic_error instead. */ template +#include -#include +#include +#include #include #include @@ -149,7 +150,7 @@ protected: /*! * \brief Module for the modular fluid state which does not the - * temperatures but throws Dune::InvalidState instead. + * temperatures but throws std::logic_error instead. */ template +#include -#include +#include +#include #include @@ -90,7 +91,7 @@ protected: /*! * \brief Module for the modular fluid state which does not the - * viscositys but throws Dune::InvalidState instead. + * viscositys but throws std::logic_error instead. */ template +#include #include namespace Opm { diff --git a/opm/material/fluidstates/modularfluidstate.hh b/opm/material/fluidstates/modularfluidstate.hh index 3a57e50fc..faf7d2640 100644 --- a/opm/material/fluidstates/modularfluidstate.hh +++ b/opm/material/fluidstates/modularfluidstate.hh @@ -32,7 +32,7 @@ #include "fluidstateviscositymodules.hh" #include "fluidstateenthalpymodules.hh" -#include +#include #include namespace Opm { diff --git a/opm/material/fluidstates/nonequilibriumfluidstate.hh b/opm/material/fluidstates/nonequilibriumfluidstate.hh index a0bcfc089..b30a23925 100644 --- a/opm/material/fluidstates/nonequilibriumfluidstate.hh +++ b/opm/material/fluidstates/nonequilibriumfluidstate.hh @@ -28,7 +28,7 @@ #include "modularfluidstate.hh" -#include +#include #include namespace Opm { diff --git a/opm/material/fluidstates/pressureoverlayfluidstate.hh b/opm/material/fluidstates/pressureoverlayfluidstate.hh index 498d306f5..77ab18e75 100644 --- a/opm/material/fluidstates/pressureoverlayfluidstate.hh +++ b/opm/material/fluidstates/pressureoverlayfluidstate.hh @@ -23,7 +23,7 @@ #ifndef OPM_PRESSURE_OVERLAY_FLUID_STATE_HH #define OPM_PRESSURE_OVERLAY_FLUID_STATE_HH -#include +#include #include diff --git a/opm/material/fluidstates/saturationoverlayfluidstate.hh b/opm/material/fluidstates/saturationoverlayfluidstate.hh index d17fc9e9f..b83bf7168 100644 --- a/opm/material/fluidstates/saturationoverlayfluidstate.hh +++ b/opm/material/fluidstates/saturationoverlayfluidstate.hh @@ -23,7 +23,7 @@ #ifndef OPM_SATURATION_OVERLAY_FLUID_STATE_HH #define OPM_SATURATION_OVERLAY_FLUID_STATE_HH -#include +#include #include diff --git a/opm/material/fluidstates/temperatureoverlayfluidstate.hh b/opm/material/fluidstates/temperatureoverlayfluidstate.hh index d5f921502..5084ce53d 100644 --- a/opm/material/fluidstates/temperatureoverlayfluidstate.hh +++ b/opm/material/fluidstates/temperatureoverlayfluidstate.hh @@ -23,7 +23,7 @@ #ifndef OPM_TEMPERATURE_OVERLAY_FLUID_STATE_HH #define OPM_TEMPERATURE_OVERLAY_FLUID_STATE_HH -#include +#include namespace Opm { diff --git a/opm/material/fluidsystems/1pfluidsystem.hh b/opm/material/fluidsystems/1pfluidsystem.hh index 59db3992f..c01cecfc1 100644 --- a/opm/material/fluidsystems/1pfluidsystem.hh +++ b/opm/material/fluidsystems/1pfluidsystem.hh @@ -33,8 +33,6 @@ #include #include -#include - #include #include diff --git a/opm/material/fluidsystems/2pimmisciblefluidsystem.hh b/opm/material/fluidsystems/2pimmisciblefluidsystem.hh index f7e6d1d25..e72e7676f 100644 --- a/opm/material/fluidsystems/2pimmisciblefluidsystem.hh +++ b/opm/material/fluidsystems/2pimmisciblefluidsystem.hh @@ -31,9 +31,6 @@ #include #include -#include - - #include "basefluidsystem.hh" #include "nullparametercache.hh" diff --git a/opm/material/fluidsystems/basefluidsystem.hh b/opm/material/fluidsystems/basefluidsystem.hh index 58d5dd630..a7421cbce 100644 --- a/opm/material/fluidsystems/basefluidsystem.hh +++ b/opm/material/fluidsystems/basefluidsystem.hh @@ -25,9 +25,10 @@ #include "nullparametercache.hh" -#include +#include +#include -#include +#include namespace Opm { @@ -62,8 +63,8 @@ public: */ static char *phaseName(int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a phaseName() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a phaseName() method!"); } /*! @@ -73,8 +74,8 @@ public: */ static bool isLiquid(int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a isLiquid() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isLiquid() method!"); } /*! @@ -93,8 +94,8 @@ public: */ static bool isIdealMixture(int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a isIdealMixture() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isIdealMixture() method!"); } /*! @@ -108,8 +109,8 @@ public: */ static bool isCompressible(int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a isCompressible() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isCompressible() method!"); } /*! @@ -120,8 +121,8 @@ public: */ static bool isIdealGas(int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a isIdealGas() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isIdealGas() method!"); } /*! @@ -131,8 +132,8 @@ public: */ static const char *componentName(int compIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a componentName() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a componentName() method!"); } /*! @@ -142,8 +143,8 @@ public: */ static Scalar molarMass(int compIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a molarMass() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a molarMass() method!"); } /*! @@ -163,8 +164,8 @@ public: const ParameterCache ¶mCache, int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, - "The fluid system '" << Dune::className() << "' does not provide a density() method!"); + OPM_THROW(std::runtime_error, + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a density() method!"); } /*! @@ -187,7 +188,7 @@ public: int phaseIdx, int compIdx) { - DUNE_THROW(Dune::NotImplemented, "The fluid system '" << Dune::className() << "' does not provide a fugacityCoefficient() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a fugacityCoefficient() method!"); } /*! @@ -201,7 +202,7 @@ public: const ParameterCache ¶mCache, int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, "The fluid system '" << Dune::className() << "' does not provide a viscosity() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a viscosity() method!"); } /*! @@ -227,7 +228,7 @@ public: int phaseIdx, int compIdx) { - DUNE_THROW(Dune::NotImplemented, "The fluid system '" << Dune::className() << "' does not provide a diffusionCoefficient() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a diffusionCoefficient() method!"); } /*! @@ -242,7 +243,7 @@ public: const ParameterCache ¶mCache, int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, "The fluid system '" << Dune::className() << "' does not provide an enthalpy() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide an enthalpy() method!"); } /*! @@ -256,7 +257,7 @@ public: const ParameterCache ¶mCache, int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, "The fluid system '" << Dune::className() << "' does not provide a thermalConductivity() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a thermalConductivity() method!"); } /*! @@ -270,7 +271,7 @@ public: const ParameterCache ¶mCache, int phaseIdx) { - DUNE_THROW(Dune::NotImplemented, "The fluid system '" << Dune::className() << "' does not provide a heatCapacity() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a heatCapacity() method!"); } }; diff --git a/opm/material/fluidsystems/blackoilfluidsystem.hh b/opm/material/fluidsystems/blackoilfluidsystem.hh index 9a5d4ef46..dbdb08b5f 100644 --- a/opm/material/fluidsystems/blackoilfluidsystem.hh +++ b/opm/material/fluidsystems/blackoilfluidsystem.hh @@ -23,8 +23,9 @@ #ifndef OPM_BLACK_OIL_FLUID_SYSTEM_HH #define OPM_BLACK_OIL_FLUID_SYSTEM_HH -#include -#include +#include +#include +#include #include #include @@ -47,7 +48,7 @@ template class BlackOil : public BaseFluidSystem > { - typedef Opm::Spline Spline; + typedef Opm::Spline Spline; typedef std::vector > SplineSamplingPoints; public: @@ -72,7 +73,7 @@ public: * \copydoc BaseFluidSystem::init * * \attention For this fluid system, this method just throws a - * Dune::InvalidStateException as there is no + * std::logic_error as there is no * way to generically calculate the required black oil * parameters. Instead of this method, use * \code @@ -83,7 +84,7 @@ public: */ static void init() { - DUNE_THROW(Dune::InvalidStateException, "No generic init() method for this fluid system. The black-oil fluid system must be initialized with:\n" + OPM_THROW(std::logic_error, "No generic init() method for this fluid system. The black-oil fluid system must be initialized with:\n" << " FluidSystem::initBegin()\n" << " // set black oil parameters\n" << " FluidSystem::initEnd()\n"); @@ -399,7 +400,7 @@ public: } } - DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::fugacityCoefficient @@ -419,7 +420,7 @@ public: case oPhaseIdx: return fugCoefficientInOil(compIdx, p); } - DUNE_THROW(Dune::InvalidStateException, "Unhandled phase or component index"); + OPM_THROW(std::logic_error, "Unhandled phase or component index"); } //! \copydoc BaseFluidSystem::viscosity @@ -438,7 +439,7 @@ public: case gPhaseIdx: return gasViscosity_(p); } - DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx); } /*! @@ -593,7 +594,7 @@ public: return pSat; } - DUNE_THROW(NumericalProblem, "Could find the oil saturation pressure for X_o^g = " << X_oG); + OPM_THROW(NumericalProblem, "Could find the oil saturation pressure for X_o^g = " << X_oG); } diff --git a/opm/material/fluidsystems/brineco2fluidsystem.hh b/opm/material/fluidsystems/brineco2fluidsystem.hh index 8c4d39a96..3add1edb0 100644 --- a/opm/material/fluidsystems/brineco2fluidsystem.hh +++ b/opm/material/fluidsystems/brineco2fluidsystem.hh @@ -463,12 +463,12 @@ private: Valgrind::CheckDefined(xlCO2); if(T < 273.15) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Liquid density for Brine and CO2 is only " "defined above 273.15K (is " << T << "K)"); } if(pl >= 2.5e8) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Liquid density for Brine and CO2 is only " "defined below 250MPa (is " << pl << "Pa)"); } diff --git a/opm/material/fluidsystems/h2oairfluidsystem.hh b/opm/material/fluidsystems/h2oairfluidsystem.hh index aac059f6e..aeab24ac6 100644 --- a/opm/material/fluidsystems/h2oairfluidsystem.hh +++ b/opm/material/fluidsystems/h2oairfluidsystem.hh @@ -32,8 +32,9 @@ #include #include #include -#include -#include +#include +#include +#include #include #include @@ -85,7 +86,7 @@ public: case lPhaseIdx: return "liquid"; case gPhaseIdx: return "gas"; }; - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::isLiquid @@ -146,7 +147,7 @@ public: case H2OIdx: return H2O::name(); case AirIdx: return Air::name(); }; - DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + OPM_THROW(std::logic_error, "Invalid component index " << compIdx); } //! \copydoc BaseFluidSystem::molarMass @@ -320,7 +321,7 @@ public: H2O::gasDensity(T, partialPressureH2O) + Air::gasDensity(T, partialPressureAir); } - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::viscosity @@ -386,7 +387,7 @@ public: return muResult; } } - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::fugacityCoefficient @@ -459,7 +460,7 @@ public: fluidState.massFraction(gPhaseIdx, AirIdx); return result; } - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::thermalConductivity diff --git a/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh b/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh index 7f204465f..87ab59ab9 100644 --- a/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh +++ b/opm/material/fluidsystems/h2oairmesitylenefluidsystem.hh @@ -168,7 +168,7 @@ public: case nPhaseIdx: return "n"; case gPhaseIdx: return "g";; }; - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::componentName @@ -179,7 +179,7 @@ public: case airIdx: return Air::name(); case NAPLIdx: return NAPL::name(); }; - DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + OPM_THROW(std::logic_error, "Invalid component index " << compIdx); } //! \copydoc BaseFluidSystem::molarMass @@ -193,7 +193,7 @@ public: : (compIdx == NAPLIdx) ? NAPL::molarMass() : -1e10; - //DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + //OPM_THROW(std::logic_error, "Invalid component index " << compIdx); } //! \copydoc BaseFluidSystem::density @@ -335,7 +335,7 @@ public: if (compIdx==NAPLIdx) return (1 - xgw)/(xga/diffAW + xgc/diffWC); else if (compIdx==H2OIdx) return (1 - xgc)/(xgw/diffWC + xga/diffAC); - else if (compIdx==airIdx) DUNE_THROW(Dune::InvalidStateException, + else if (compIdx==airIdx) OPM_THROW(std::logic_error, "Diffusivity of air in the gas phase " "is constraint by sum of diffusive fluxes = 0 !\n"); } @@ -356,13 +356,13 @@ public: diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl); return diffCont; case H2OIdx: - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "Diffusivity of water in the water phase " "is constraint by sum of diffusive fluxes = 0 !\n"); }; } else if (phaseIdx==nPhaseIdx) { - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "Diffusion coefficients of " "substances in liquid phase are undefined!\n"); } @@ -441,7 +441,7 @@ public: return result; } - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::thermalConductivity diff --git a/opm/material/fluidsystems/h2oairxylenefluidsystem.hh b/opm/material/fluidsystems/h2oairxylenefluidsystem.hh index 20f2edf85..690822892 100644 --- a/opm/material/fluidsystems/h2oairxylenefluidsystem.hh +++ b/opm/material/fluidsystems/h2oairxylenefluidsystem.hh @@ -130,7 +130,7 @@ public: case nPhaseIdx: return "n"; case gPhaseIdx: return "g";; }; - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::componentName @@ -141,7 +141,7 @@ public: case airIdx: return Air::name(); case NAPLIdx: return NAPL::name(); }; - DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + OPM_THROW(std::logic_error, "Invalid component index " << compIdx); } //! \copydoc BaseFluidSystem::molarMass @@ -292,7 +292,7 @@ public: if (compIdx==NAPLIdx) return (1.- xgw)/(xga/diffAW + xgc/diffWC); else if (compIdx==H2OIdx) return (1.- xgc)/(xgw/diffWC + xga/diffAC); - else if (compIdx==airIdx) DUNE_THROW(Dune::InvalidStateException, + else if (compIdx==airIdx) OPM_THROW(std::logic_error, "Diffusivity of air in the gas phase " "is constraint by sum of diffusive fluxes = 0 !\n"); } else if (phaseIdx==wPhaseIdx){ @@ -312,13 +312,13 @@ public: diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl); return diffCont; case H2OIdx: - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "Diffusivity of water in the water phase " "is constraint by sum of diffusive fluxes = 0 !\n"); }; } else if (phaseIdx==nPhaseIdx) { - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "Diffusion coefficients of " "substances in liquid phase are undefined!\n"); } @@ -394,7 +394,7 @@ public: return result; } - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } private: diff --git a/opm/material/fluidsystems/h2on2fluidsystem.hh b/opm/material/fluidsystems/h2on2fluidsystem.hh index fb0a44a72..11a2940b1 100644 --- a/opm/material/fluidsystems/h2on2fluidsystem.hh +++ b/opm/material/fluidsystems/h2on2fluidsystem.hh @@ -32,8 +32,9 @@ #include #include #include -#include -#include +#include +#include +#include #include #include diff --git a/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh b/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh index 3ab6108eb..2f1e393c9 100644 --- a/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh +++ b/opm/material/fluidsystems/h2on2liquidphasefluidsystem.hh @@ -33,8 +33,9 @@ #include #include #include -#include -#include +#include +#include +#include #include #include diff --git a/opm/material/fluidsystems/spe5fluidsystem.hh b/opm/material/fluidsystems/spe5fluidsystem.hh index 7d6762279..2a0ff41eb 100644 --- a/opm/material/fluidsystems/spe5fluidsystem.hh +++ b/opm/material/fluidsystems/spe5fluidsystem.hh @@ -26,7 +26,7 @@ #include "basefluidsystem.hh" #include "spe5parametercache.hh" -#include +#include #include #include @@ -424,7 +424,7 @@ protected: case C10Idx: return 1e10; case C15Idx: return 1e10; case C20Idx: return 1e10; - default: DUNE_THROW(Dune::InvalidStateException, "Unknown component index " << compIdx); + default: OPM_THROW(std::logic_error, "Unknown component index " << compIdx); } } }; diff --git a/opm/material/fluidsystems/spe5parametercache.hh b/opm/material/fluidsystems/spe5parametercache.hh index 555ca9c5d..0af95889a 100644 --- a/opm/material/fluidsystems/spe5parametercache.hh +++ b/opm/material/fluidsystems/spe5parametercache.hh @@ -110,7 +110,7 @@ public: case oPhaseIdx: return oilPhaseParams_.a(); case gPhaseIdx: return gasPhaseParams_.a(); default: - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "The a() parameter is only defined for " "oil and gas phases"); }; @@ -128,7 +128,7 @@ public: case oPhaseIdx: return oilPhaseParams_.b(); case gPhaseIdx: return gasPhaseParams_.b(); default: - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "The b() parameter is only defined for " "oil and gas phases"); }; @@ -149,7 +149,7 @@ public: case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a(); case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a(); default: - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "The a() parameter is only defined for " "oil and gas phases"); }; @@ -169,7 +169,7 @@ public: case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b(); case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b(); default: - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "The b() parameter is only defined for " "oil and gas phases"); }; diff --git a/opm/material/heatconduction/dummyheatconductionlaw.hh b/opm/material/heatconduction/dummyheatconductionlaw.hh index 16a4c6a08..df121e238 100644 --- a/opm/material/heatconduction/dummyheatconductionlaw.hh +++ b/opm/material/heatconduction/dummyheatconductionlaw.hh @@ -23,17 +23,18 @@ #ifndef OPM_DUMMY_HEATCONDUCTION_LAW_HH #define OPM_DUMMY_HEATCONDUCTION_LAW_HH -#include +#include +#include namespace Opm { /*! * \ingroup material * - * \brief Implements a dumm law for heat conduction to which isothermal models - * can fall back to + * \brief Implements a dummy law for heat conduction to which isothermal models + * can fall back to. * - * If any method of this law is called, it throws an excetion + * If any method of this law is called, it throws std::logic_error. */ template class DummyHeatConductionLaw @@ -52,7 +53,7 @@ public: static Scalar heatConductivity(const Params ¶ms, const FluidState &fluidState) { - DUNE_THROW(Dune::InvalidStateException, + OPM_THROW(std::logic_error, "No heat conduction law specified!"); } }; diff --git a/opm/material/heatconduction/fluidconduction.hh b/opm/material/heatconduction/fluidconduction.hh index 963230db0..8e16fae50 100644 --- a/opm/material/heatconduction/fluidconduction.hh +++ b/opm/material/heatconduction/fluidconduction.hh @@ -25,7 +25,7 @@ #include "fluidconductionparams.hh" -#include +#include #include namespace Opm diff --git a/opm/material/heatconduction/somerton.hh b/opm/material/heatconduction/somerton.hh index 1a50a33ea..89f73e52a 100644 --- a/opm/material/heatconduction/somerton.hh +++ b/opm/material/heatconduction/somerton.hh @@ -25,7 +25,7 @@ #include "somertonparams.hh" -#include +#include #include namespace Opm diff --git a/opm/common/statictabulated2dfunction.hh b/opm/material/statictabulated2dfunction.hh similarity index 95% rename from opm/common/statictabulated2dfunction.hh rename to opm/material/statictabulated2dfunction.hh index ba57abf12..e10c0924a 100644 --- a/opm/common/statictabulated2dfunction.hh +++ b/opm/material/statictabulated2dfunction.hh @@ -24,8 +24,9 @@ #ifndef OPM_STATIC_TABULATED_2D_FUNCTION_HH #define OPM_STATIC_TABULATED_2D_FUNCTION_HH -#include -#include +#include +#include +#include #include @@ -92,7 +93,7 @@ public: #if !defined NDEBUG if (i < 0 || i >= Traits::numX || j < 0 || j >= Traits::numY) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Attempt to access element (" << i << ", " << j << ") on a " << Traits::name << " table of size (" diff --git a/opm/common/tabulated2dfunction.hh b/opm/material/tabulated2dfunction.hh similarity index 97% rename from opm/common/tabulated2dfunction.hh rename to opm/material/tabulated2dfunction.hh index 48f7d064a..cf714b411 100644 --- a/opm/common/tabulated2dfunction.hh +++ b/opm/material/tabulated2dfunction.hh @@ -24,7 +24,8 @@ #ifndef OPM_TABULATED_2D_FUNCTION_HH #define OPM_TABULATED_2D_FUNCTION_HH -#include +#include +#include #include @@ -110,7 +111,7 @@ public: #ifndef NDEBUG if (!applies(x,y)) { - DUNE_THROW(NumericalProblem, + OPM_THROW(NumericalProblem, "Attempt to get tabulated value for (" << x << ", " << y << ") on a table of extend " diff --git a/opm/common/valgrind.hh b/opm/material/valgrind.hh similarity index 94% rename from opm/common/valgrind.hh rename to opm/material/valgrind.hh index 7a70d0fe0..16535f019 100644 --- a/opm/common/valgrind.hh +++ b/opm/material/valgrind.hh @@ -23,7 +23,7 @@ #ifndef OPM_VALGRIND_HH #define OPM_VALGRIND_HH -#include +#include #if ! HAVE_VALGRIND && ! defined(DOXYGEN) namespace Valgrind @@ -150,7 +150,7 @@ template inline void SetUndefined(const T &value) { #if !defined NDEBUG && HAVE_VALGRIND - auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(&value, sizeof(T)); + auto OPM_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(&value, sizeof(T)); #endif } @@ -176,7 +176,7 @@ 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)); + auto OPM_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(value, size*sizeof(T)); #endif } @@ -200,7 +200,7 @@ template inline void SetDefined(const T &value) { #if !defined NDEBUG && HAVE_VALGRIND - auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(&value, sizeof(T)); + auto OPM_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(&value, sizeof(T)); #endif } @@ -226,7 +226,7 @@ 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)); + auto OPM_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(value, n*sizeof(T)); #endif } @@ -250,7 +250,7 @@ template inline void SetNoAccess(const T &value) { #if !defined NDEBUG && HAVE_VALGRIND - auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(&value, sizeof(T)); + auto OPM_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(&value, sizeof(T)); #endif } @@ -276,7 +276,7 @@ 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)); + auto OPM_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(value, size*sizeof(T)); #endif } diff --git a/tests/material/fluidsystems/checkfluidsystem.hh b/tests/material/fluidsystems/checkfluidsystem.hh index 9cb642478..6eb173837 100644 --- a/tests/material/fluidsystems/checkfluidsystem.hh +++ b/tests/material/fluidsystems/checkfluidsystem.hh @@ -41,7 +41,8 @@ #include #include -#include +#include +#include #include #include @@ -91,7 +92,7 @@ public: { assert(allowTemperature_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::temperature(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::temperature(phaseIdx); return 1e100; } @@ -99,7 +100,7 @@ public: { assert(allowPressure_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::pressure(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::pressure(phaseIdx); return 1e100; } @@ -107,7 +108,7 @@ public: { assert(allowComposition_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::moleFraction(phaseIdx, compIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::moleFraction(phaseIdx, compIdx); return 1e100; } @@ -115,7 +116,7 @@ public: { assert(allowComposition_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::massFraction(phaseIdx, compIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::massFraction(phaseIdx, compIdx); return 1e100; } @@ -123,7 +124,7 @@ public: { assert(allowComposition_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::averageMolarMass(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::averageMolarMass(phaseIdx); return 1e100; } @@ -131,7 +132,7 @@ public: { assert(allowDensity_ && allowComposition_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::molarity(phaseIdx, compIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::molarity(phaseIdx, compIdx); return 1e100; } @@ -139,7 +140,7 @@ public: { assert(allowDensity_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::molarDensity(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::molarDensity(phaseIdx); return 1e100; } @@ -147,7 +148,7 @@ public: { assert(allowDensity_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::molarVolume(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::molarVolume(phaseIdx); return 1e100; } @@ -155,49 +156,49 @@ public: { assert(allowDensity_); assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx); - DUNE_UNUSED Scalar tmp = BaseFluidState::density(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::density(phaseIdx); return 1e100; } Scalar saturation(int phaseIdx) const { assert(false); - DUNE_UNUSED Scalar tmp = BaseFluidState::saturation(phaseIdx); + OPM_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); + OPM_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); + OPM_UNUSED Scalar tmp = BaseFluidState::fugacityCoefficient(phaseIdx, compIdx); return 1e100; } Scalar enthalpy(int phaseIdx) const { assert(false); - DUNE_UNUSED Scalar tmp = BaseFluidState::enthalpy(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::enthalpy(phaseIdx); return 1e100; } Scalar internalEnergy(int phaseIdx) const { assert(false); - DUNE_UNUSED Scalar tmp = BaseFluidState::internalEnergy(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::internalEnergy(phaseIdx); return 1e100; } Scalar viscosity(int phaseIdx) const { assert(false); - DUNE_UNUSED Scalar tmp = BaseFluidState::viscosity(phaseIdx); + OPM_UNUSED Scalar tmp = BaseFluidState::viscosity(phaseIdx); return 1e100; } @@ -222,7 +223,7 @@ void checkFluidState(const BaseFluidState &fs) // make sure the fluid state provides all mandatory methods while (false) { - Scalar DUNE_UNUSED val; + Scalar OPM_UNUSED val; val = fs.temperature(/*phaseIdx=*/0); val = fs.pressure(/*phaseIdx=*/0); @@ -248,7 +249,7 @@ void checkFluidState(const BaseFluidState &fs) template void checkFluidSystem() { - std::cout << "Testing fluid system '" << Dune::className() << "'\n"; + std::cout << "Testing fluid system '" << Opm::className() << "'\n"; // make sure the fluid system provides the number of phases and // the number of components @@ -282,7 +283,7 @@ void checkFluidSystem() // some value to make sure the return values of the fluid system // are convertible to scalars - Scalar DUNE_UNUSED val; + Scalar OPM_UNUSED val; // actually check the fluid system API try { FluidSystem::init(); } catch (...) {}; @@ -311,15 +312,15 @@ void checkFluidSystem() // 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); + std::string OPM_UNUSED name = FluidSystem::phaseName(phaseIdx); + bool OPM_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::string OPM_UNUSED name = FluidSystem::componentName(compIdx); } std::cout << "----------------------------------\n"; diff --git a/tests/material/fluidsystems/test_fluidsystems.cc b/tests/material/fluidsystems/test_fluidsystems.cc index 1843505ed..8bec37843 100644 --- a/tests/material/fluidsystems/test_fluidsystems.cc +++ b/tests/material/fluidsystems/test_fluidsystems.cc @@ -26,13 +26,6 @@ #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 @@ -54,7 +47,7 @@ // include the tables for CO2 which are delivered with opm-material by // default -#include +#include namespace Opm { namespace FluidSystemsTest { #include @@ -69,9 +62,6 @@ int main(int argc, char **argv) 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; diff --git a/tests/material/pengrobinson/test_pengrobinson.cc b/tests/material/pengrobinson/test_pengrobinson.cc index 71f443b7c..f2f436b86 100644 --- a/tests/material/pengrobinson/test_pengrobinson.cc +++ b/tests/material/pengrobinson/test_pengrobinson.cc @@ -113,7 +113,7 @@ Scalar bringOilToSurface(FluidState &surfaceFluidState, Scalar alpha, const Flui ComponentVector tmpMolarities; for (int i = 0;; ++i) { if (i >= 20) - DUNE_THROW(Opm::NumericalProblem, + OPM_THROW(Opm::NumericalProblem, "Newton method did not converge after 20 iterations"); // calculate the deviation from the standard pressure