OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #345 from andlaus/import_splines

Browse Source 
Import splines
This commit is contained in:
Atgeirr Flø Rasmussen 2013-09-19 10:26:02 -07:00
commit 86d09115b9
6 changed files with 3149 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
CMakeLists_files.cmake
View File

@ -139,6 +139,7 @@ list (APPEND MAIN_SOURCE_FILES
# originally generated with the command:
# find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
list (APPEND TEST_SOURCE_FILES
tests/test_spline.cpp
tests/test_dgbasis.cpp
tests/test_cartgrid.cpp
tests/test_cubic.cpp
@ -331,12 +332,16 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/Factory.hpp
opm/core/utility/MonotCubicInterpolator.hpp
opm/core/utility/NonuniformTableLinear.hpp
opm/core/utility/PolynomialUtils.hpp
opm/core/utility/RootFinders.hpp
opm/core/utility/SparseTable.hpp
opm/core/utility/SparseVector.hpp
opm/core/utility/Spline.hpp
opm/core/utility/StopWatch.hpp
opm/core/utility/TridiagonalMatrix.hpp
opm/core/utility/UniformTableLinear.hpp
opm/core/utility/Units.hpp
opm/core/utility/Unused.hpp
opm/core/utility/VelocityInterpolation.hpp
opm/core/utility/WachspressCoord.hpp
opm/core/utility/buildUniformMonotoneTable.hpp

310
opm/core/utility/PolynomialUtils.hpp Normal file
View File

@ -0,0 +1,310 @@
// -*- 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-2013 by Andreas Lauser *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Define some often used mathematical functions
*/
#ifndef OPM_MATH_HH
#define OPM_MATH_HH
#include <cmath>
#include <algorithm>
namespace Opm
{
/*!
* \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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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;
}
}
#endif

1703
opm/core/utility/Spline.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

849
opm/core/utility/TridiagonalMatrix.hpp Normal file
View File

@ -0,0 +1,849 @@
// -*- 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-2013 by Andreas Lauser *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTBILITY 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \copydetails Opm::TridiagonalMatrix
*/
#ifndef OPM_TRIDIAGONAL_MATRIX_HH
#define OPM_TRIDIAGONAL_MATRIX_HH
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <assert.h>
namespace Opm {
/*!
* \brief Provides a tridiagonal matrix that also supports non-zero
* entries in the upper right and lower left
*
* The entries in the lower left and upper right are supported to make
* implementing periodic systems easy.
*
* The API of this class is designed to be close to the one used by
* the DUNE matrix classes.
*/
template <class Scalar>
class TridiagonalMatrix
{
struct TridiagRow_ {
TridiagRow_(TridiagonalMatrix &m, size_t rowIdx)
: matrix_(m)
, rowIdx_(rowIdx)
{};
Scalar &operator[](size_t colIdx)
{ return matrix_.at(rowIdx_, colIdx); }
Scalar operator[](size_t colIdx) const
{ return matrix_.at(rowIdx_, colIdx); }
/*!
* \brief Prefix increment operator
*/
TridiagRow_ &operator++()
{ ++ rowIdx_; return *this; }
/*!
* \brief Prefix decrement operator
*/
TridiagRow_ &operator--()
{ -- rowIdx_; return *this; }
/*!
* \brief Comparision operator
*/
bool operator==(const TridiagRow_ &other) const
{ return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; }
/*!
* \brief Comparision operator
*/
bool operator!=(const TridiagRow_ &other) const
{ return !operator==(other); }
/*!
* \brief Dereference operator
*/
TridiagRow_ &operator*()
{ return *this; }
/*!
* \brief Return the row index of the this row.
*
* 0 is the first row.
*/
size_t index() const
{ return rowIdx_; }
private:
TridiagonalMatrix &matrix_;
mutable size_t rowIdx_;
};
public:
typedef Scalar FieldType;
typedef TridiagRow_ RowType;
typedef size_t SizeType;
typedef TridiagRow_ iterator;
typedef TridiagRow_ const_iterator;
explicit TridiagonalMatrix(int numRows = 0)
{
resize(numRows);
};
TridiagonalMatrix(int numRows, Scalar value)
{
resize(numRows);
this->operator=(value);
};
/*!
* \brief Copy constructor.
*/
TridiagonalMatrix(const TridiagonalMatrix &source)
{ *this = source; };
/*!
* \brief Return the number of rows/columns of the matrix.
*/
size_t size() const
{ return diag_[0].size(); }
/*!
* \brief Return the number of rows of the matrix.
*/
size_t rows() const
{ return size(); }
/*!
* \brief Return the number of columns of the matrix.
*/
size_t cols() const
{ return size(); }
/*!
* \brief Change the number of rows of the matrix.
*/
void resize(size_t n)
{
if (n == size())
return;
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx].resize(n);
}
/*!
* \brief Access an entry.
*/
Scalar &at(size_t rowIdx, size_t colIdx)
{
size_t n = size();
// special cases
if (n > 2) {
if (rowIdx == 0 && colIdx == n - 1)
return diag_[2][0];
if (rowIdx == n - 1 && colIdx == 0)
return diag_[0][n - 1];
}
int diagIdx = 1 + colIdx - rowIdx;
// make sure that the requested column is in range
assert(0 <= diagIdx && diagIdx < 3);
return diag_[diagIdx][colIdx];
}
/*!
* \brief Access an entry.
*/
Scalar at(size_t rowIdx, size_t colIdx) const
{
int n = size();
// special cases
if (rowIdx == 0 && colIdx == n - 1)
return diag_[2][0];
if (rowIdx == n - 1 && colIdx == 0)
return diag_[0][n - 1];
int diagIdx = 1 + colIdx - rowIdx;
// make sure that the requested column is in range
assert(0 <= diagIdx && diagIdx < 3);
return diag_[diagIdx][colIdx];
}
/*!
* \brief Assignment operator from another tridiagonal matrix.
*/
TridiagonalMatrix &operator=(const TridiagonalMatrix &source)
{
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx] = source.diag_[diagIdx];
return *this;
}
/*!
* \brief Assignment operator from a Scalar.
*/
TridiagonalMatrix &operator=(Scalar value)
{
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx].assign(size(), value);
return *this;
}
/*!
* \begin Iterator for the first row
*/
iterator begin()
{ return TridiagRow_(*this, 0); }
/*!
* \begin Const iterator for the first row
*/
const_iterator begin() const
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), 0); }
/*!
* \begin Const iterator for the next-to-last row
*/
const_iterator end() const
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), size()); }
/*!
* \brief Row access operator.
*/
TridiagRow_ operator[](size_t rowIdx)
{ return TridiagRow_(*this, rowIdx); }
/*!
* \brief Row access operator.
*/
const TridiagRow_ operator[](size_t rowIdx) const
{ return TridiagRow_(*this, rowIdx); }
/*!
* \brief Multiplication with a Scalar
*/
TridiagonalMatrix &operator*=(Scalar alpha)
{
int n = size();
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) {
for (int i = 0; i < n; ++i) {
diag_[diagIdx][i] *= alpha;
}
}
return *this;
}
/*!
* \brief Division by a Scalar
*/
TridiagonalMatrix &operator/=(Scalar alpha)
{
alpha = 1.0/alpha;
int n = size();
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) {
for (int i = 0; i < n; ++i) {
diag_[diagIdx][i] *= alpha;
}
}
return *this;
}
/*!
* \brief Subtraction operator
*/
TridiagonalMatrix &operator-=(const TridiagonalMatrix &other)
{ return axpy(-1.0, other); }
/*!
* \brief Addition operator
*/
TridiagonalMatrix &operator+=(const TridiagonalMatrix &other)
{ return axpy(1.0, other); }
/*!
* \brief Multiply and add the matrix entries of another
* tridiagonal matrix.
*
* This means that
* \code
* A.axpy(alpha, B)
* \endcode
* is equivalent to
* \code
* A += alpha*C
* \endcode
*/
TridiagonalMatrix &axpy(Scalar alpha, const TridiagonalMatrix &other)
{
assert(size() == other.size());
int n = size();
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
for (int i = 0; i < n; ++ i)
diag_[diagIdx][i] += alpha * other[diagIdx][i];
return *this;
}
/*!
* \brief Matrix-vector product
*
* This means that
* \code
* A.mv(x, y)
* \endcode
* is equivalent to
* \code
* y = A*x
* \endcode
*/
template<class Vector>
void mv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] =
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] =
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] =
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Additive matrix-vector product
*
* This means that
* \code
* A.umv(x, y)
* \endcode
* is equivalent to
* \code
* y += A*x
* \endcode
*/
template<class Vector>
void umv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] +=
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] +=
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Subtractive matrix-vector product
*
* This means that
* \code
* A.mmv(x, y)
* \endcode
* is equivalent to
* \code
* y -= A*x
* \endcode
*/
template<class Vector>
void mmv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] -=
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] -=
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] -=
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Scaled additive matrix-vector product
*
* This means that
* \code
* A.usmv(x, y)
* \endcode
* is equivalent to
* \code
* y += alpha*(A*x)
* \endcode
*/
template<class Vector>
void usmv(Scalar alpha, const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
alpha*(
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1]);
}
// rows 0 and n-1
dest[0] +=
alpha*(
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1]);
dest[n - 1] +=
alpha*(
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1]);
}
/*!
* \brief Transposed matrix-vector product
*
* This means that
* \code
* A.mtv(x, y)
* \endcode
* is equivalent to
* \code
* y = A^T*x
* \endcode
*/
template<class Vector>
void mtv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] =
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] =
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] =
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Transposed additive matrix-vector product
*
* This means that
* \code
* A.umtv(x, y)
* \endcode
* is equivalent to
* \code
* y += A^T*x
* \endcode
*/
template<class Vector>
void umtv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] +=
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] +=
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Transposed subtractive matrix-vector product
*
* This means that
* \code
* A.mmtv(x, y)
* \endcode
* is equivalent to
* \code
* y -= A^T*x
* \endcode
*/
template<class Vector>
void mmtv (const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] -=
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] -=
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] -=
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Transposed scaled additive matrix-vector product
*
* This means that
* \code
* A.umtv(alpha, x, y)
* \endcode
* is equivalent to
* \code
* y += alpha*A^T*x
* \endcode
*/
template<class Vector>
void usmtv(Scalar alpha, const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
alpha*(
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1]);
}
// rows 0 and n-1
dest[0] +=
alpha*(
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1]);
dest[n - 1] +=
alpha*(
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1]);
}
/*!
* \brief Calculate the frobenius norm
*
* i.e., the square root of the sum of all squared entries. This
* corresponds to the euclidean norm for vectors.
*/
Scalar frobeniusNorm() const
{ return std::sqrt(frobeniusNormSquared()); }
/*!
* \brief Calculate the squared frobenius norm
*
* i.e., the sum of all squared entries.
*/
Scalar frobeniusNormSquared() const
{
Scalar result = 0;
int n = size();
for (int i = 0; i < n; ++ i)
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
result += diag_[diagIdx][i];
return result;
}
/*!
* \brief Calculate the infinity norm
*
* i.e., the maximum of the sum of the absolute values of all rows.
*/
Scalar infinityNorm() const
{
Scalar result=0;
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
result = std::max(result,
std::abs(diag_[0][i - 1]) +
std::abs(diag_[1][i]) +
std::abs(diag_[2][i + 1]));
}
// rows 0 and n-1
result = std::max(result,
std::abs(diag_[1][0]) +
std::abs(diag_[2][1]) +
std::abs(diag_[2][0]));
result = std::max(result,
std::abs(diag_[0][n-1]) +
std::abs(diag_[0][n-2]) +
std::abs(diag_[1][n-2]));
return result;
}
/*!
* \brief Calculate the solution for a linear system of equations
*
* i.e., calculate x, so that it solves Ax = b, where A is a
* tridiagonal matrix.
*/
template <class XVector, class BVector>
void solve(XVector &x, const BVector &b) const
{
if (size() > 2 && diag_[2][0] != 0)
solveWithUpperRight_(x, b);
else
solveWithoutUpperRight_(x, b);
}
/*!
* \brief Print the matrix to a given output stream.
*/
void print(std::ostream &os = std::cout) const
{
int n = size();
// row 0
os << at(0, 0) << "\t"
<< at(0, 1) << "\t";
if (n > 3)
os << "\t";
if (n > 2)
os << at(0, n-1);
os << "\n";
// row 1 .. n - 2
for (int rowIdx = 1; rowIdx < n-1; ++rowIdx) {
if (rowIdx > 1)
os << "\t";
if (rowIdx == n - 2)
os << "\t";
os << at(rowIdx, rowIdx - 1) << "\t"
<< at(rowIdx, rowIdx) << "\t"
<< at(rowIdx, rowIdx + 1) << "\n";
}
// row n - 1
if (n > 2)
os << at(n-1, 0) << "\t";
if (n > 3)
os << "\t";
if (n > 4)
os << "\t";
os << at(n-1, n-2) << "\t"
<< at(n-1, n-1) << "\n";
}
private:
template <class XVector, class BVector>
void solveWithUpperRight_(XVector &x, const BVector &b) const
{
size_t n = size();
std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]), lastColumn(n);
std::vector<Scalar> bStar(n);
std::copy(b.begin(), b.end(), bStar.begin());
lastColumn[0] = upperDiag[0];
// forward elimination
for (size_t i = 1; i < n; ++i) {
Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
mainDiag[i] -= alpha * upperDiag[i];
bStar[i] -= alpha * bStar[i - 1];
};
// deal with the last row if the entry on the lower left is not zero
if (lowerDiag[n - 1] != 0.0 && size() > 2) {
Scalar lastRow = lowerDiag[n - 1];
for (size_t i = 0; i < n - 1; ++i) {
Scalar alpha = lastRow/mainDiag[i];
lastRow = - alpha*upperDiag[i + 1];
bStar[n - 1] -= alpha * bStar[i];
}
mainDiag[n-1] += lastRow;
}
// backward elimination
x[n - 1] = bStar[n - 1]/mainDiag[n-1];
for (int i = n - 2; i >= 0; --i) {
x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1] - x[n-1]*lastColumn[i])/mainDiag[i];
}
}
template <class XVector, class BVector>
void solveWithoutUpperRight_(XVector &x, const BVector &b) const
{
size_t n = size();
std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]);
std::vector<Scalar> bStar(n);
std::copy(b.begin(), b.end(), bStar.begin());
// forward elimination
for (size_t i = 1; i < n; ++i) {
Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
mainDiag[i] -= alpha * upperDiag[i];
bStar[i] -= alpha * bStar[i - 1];
};
// deal with the last row if the entry on the lower left is not zero
if (lowerDiag[n - 1] != 0.0 && size() > 2) {
Scalar lastRow = lowerDiag[n - 1];
for (size_t i = 0; i < n - 1; ++i) {
Scalar alpha = lastRow/mainDiag[i];
lastRow = - alpha*upperDiag[i + 1];
bStar[n - 1] -= alpha * bStar[i];
}
mainDiag[n-1] += lastRow;
}
// backward elimination
x[n - 1] = bStar[n - 1]/mainDiag[n-1];
for (int i = n - 2; i >= 0; --i) {
x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1])/mainDiag[i];
}
}
mutable std::vector<Scalar> diag_[3];
};
} // namespace Opm
template <class Scalar>
std::ostream &operator<<(std::ostream &os, const Opm::TridiagonalMatrix<Scalar> &mat)
{
mat.print(os);
return os;
}
#endif

28
opm/core/utility/Unused.hpp Normal file
View File

@ -0,0 +1,28 @@
// -*- 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-2013 by Andreas Lauser *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef OPM_CORE_UNUSED_HH
#define OPM_CORE_UNUSED_HH
#ifndef HAS_ATTRIBUTE_UNUSED
#define OPM_UNUSED
#else
#define OPM_UNUSED __attribute__((unused))
#endif
#endif

254
tests/test_spline.cpp Normal file
View File

@ -0,0 +1,254 @@
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief This is a program to test the polynomial spline interpolation.
*
* It just prints some function to stdout. You can look at the result
* using the following commands:
*
----------- snip -----------
./test_spline > 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 "Monotonical"
----------- snap -----------
*/
#include "config.h"
#include <opm/core/utility/Spline.hpp>
#define GCC_VERSION (__GNUC__ * 10000 \
+ __GNUC_MINOR__ * 100 \
+ __GNUC_PATCHLEVEL__)
template <class Spline>
void testCommon(const Spline &sp,
const double *x,
const double *y)
{
static double eps = 1e-10;
int n = sp.numSamples();
for (int i = 0; i < n; ++i) {
// sure that we hit all sampling points
double y0 = (i>0)?sp.eval(x[i]-eps):y[0];
double y1 = sp.eval(x[i]);
double y2 = (i<n-1)?sp.eval(x[i]+eps):y[n-1];
if (std::abs(y0 - y[i]) > 100*eps || std::abs(y2 - y[i]) > 100*eps)
OPM_THROW(std::runtime_error,
"Spline seems to be discontinuous at sampling point " << i << "!");
if (std::abs(y1 - y[i]) > eps)
OPM_THROW(std::runtime_error,
"Spline does not capture sampling point " << i << "!");
// make sure the derivative is continuous (assuming that the
// second derivative is smaller than 1000)
double d1 = sp.evalDerivative(x[i]);
double d0 = (i>0)?sp.evalDerivative(x[i]-eps):d1;
double d2 = (i<n-1)?sp.evalDerivative(x[i]+eps):d1;
if (std::abs(d1 - d0) > 1000*eps || std::abs(d2 - d0) > 1000*eps)
OPM_THROW(std::runtime_error,
"Spline seems to exhibit a discontinuous derivative at sampling point " << i << "!");
}
}
template <class Spline>
void testFull(const Spline &sp,
const double *x,
const double *y,
double m0,
double m1)
{
// test the common properties of splines
testCommon(sp, x, y);
static double eps = 1e-5;
int n = sp.numSamples();
// make sure the derivative at both end points is correct
double d0 = sp.evalDerivative(x[0]);
double d1 = sp.evalDerivative(x[n-1]);
if (std::abs(d0 - m0) > eps)
OPM_THROW(std::runtime_error,
"Invalid derivative at beginning of interval: is "
<< d0 << " ought to be " << m0);
if (std::abs(d1 - m1) > eps)
OPM_THROW(std::runtime_error,
"Invalid derivative at end of interval: is "
<< d1 << " ought to be " << m1);
}
template <class Spline>
void testNatural(const Spline &sp,
const double *x,
const double *y)
{
// test the common properties of splines
testCommon(sp, x, y);
static double eps = 1e-5;
int n = sp.numSamples();
// make sure the second derivatives at both end points are 0
double d0 = sp.evalDerivative(x[0]);
double d1 = sp.evalDerivative(x[0] + eps);
double d2 = sp.evalDerivative(x[n-1] - eps);
double d3 = sp.evalDerivative(x[n-1]);
if (std::abs(d1 - d0)/eps > 1000*eps)
OPM_THROW(std::runtime_error,
"Invalid second derivative at beginning of interval: is "
<< (d1 - d0)/eps << " ought to be 0");
if (std::abs(d3 - d2)/eps > 1000*eps)
OPM_THROW(std::runtime_error,
"Invalid second derivative at end of interval: is "
<< (d3 - d2)/eps << " ought to be 0");
}
void testAll()
{
double x[] = { 0, 5, 7.5, 8.75, 9.375 };
double y[] = { 10, 0, 10, 0, 10 };
double m0 = 10;
double m1 = -10;
double points[][2] =
{
{x[0], y[0]},
{x[1], y[1]},
{x[2], y[2]},
{x[3], y[3]},
{x[4], y[4]},
};
#if GCC_VERSION >= 40500
std::initializer_list<const std::pair<double, double> > pointsInitList =
{
{x[0], y[0]},
{x[1], y[1]},
{x[2], y[2]},
{x[3], y[3]},
{x[4], y[4]},
};
#endif
std::vector<double> xVec;
std::vector<double> yVec;
std::vector<double*> pointVec;
for (int i = 0; i < 5; ++i) {
xVec.push_back(x[i]);
yVec.push_back(y[i]);
pointVec.push_back(points[i]);
}
/////////
// test spline with two sampling points
/////////
// full spline
{ Opm::Spline<double> sp(x[0], x[1], y[0], y[1], m0, m1); sp.set(x[0],x[1],y[0],y[1],m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp(2, x, y, m0, m1); sp.setXYArrays(2, x, y, m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp(2, points, m0, m1); sp.setArrayOfPoints(2, points, m0, m1); testFull(sp, x, y, m0, m1); };
/////////
// test variable length splines
/////////
// full spline
{ Opm::Spline<double> sp(5, x, y, m0, m1); sp.setXYArrays(5,x,y,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp(xVec, yVec, m0, m1); sp.setXYContainers(xVec,yVec,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp; sp.setArrayOfPoints(5,points,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp; sp.setContainerOfPoints(pointVec,m0, m1); testFull(sp, x, y, m0, m1); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double> sp; sp.setContainerOfTuples(pointsInitList,m0, m1); testFull(sp, x, y, m0, m1); };
#endif
// natural spline
{ Opm::Spline<double> sp(5, x, y); sp.setXYArrays(5,x,y); testNatural(sp, x, y); };
{ Opm::Spline<double> sp(xVec, yVec); sp.setXYContainers(xVec,yVec); testNatural(sp, x, y); };
{ Opm::Spline<double> sp; sp.setArrayOfPoints(5,points); testNatural(sp, x, y); };
{ Opm::Spline<double> sp; sp.setContainerOfPoints(pointVec); testNatural(sp, x, y); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double> sp; sp.setContainerOfTuples(pointsInitList); testNatural(sp, x, y); };
#endif
}
void plot()
{
const int numSamples = 5;
const int n = numSamples - 1;
typedef std::array<double, numSamples> FV;
double x_[] = { 0, 5, 7.5, 8.75, 10 };
double y_[] = { 10, 0, 10, 0, 10 };
double m1 = 10;
double m2 = -10;
FV &xs = *reinterpret_cast<FV*>(x_);
FV &ys = *reinterpret_cast<FV*>(y_);
Opm::Spline<double> spFull(xs, ys, m1, m2);
Opm::Spline<double> spNatural(xs, ys);
Opm::Spline<double> spPeriodic(xs, ys, /*type=*/Opm::Spline<double>::Periodic);
Opm::Spline<double> spMonotonic(xs, ys, /*type=*/Opm::Spline<double>::Monotonic);
spFull.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
spNatural.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
spPeriodic.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
spMonotonic.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
std::cerr << "Spline is monotonic: " << spFull.monotonic(x_[0], x_[n]) << "\n";
}
int main(int argc, char** argv)
{
try {
testAll();
plot();
}
catch (const std::exception &e) {
std::cout << "Caught OPM exception: " << e.what() << "\n";
}
return 0;
}