From 976ab03f687c69ed6d0da448c70539a7330a7cbd Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 21 Jun 2018 12:14:17 +0200 Subject: [PATCH] fork ISTLSolver into a legacy and non-legacy version this is necessary because only the latter can use the property system. --- CMakeLists_files.cmake | 1 + opm/autodiff/BlackoilModelEbos.hpp | 4 +- opm/autodiff/ISTLSolver.hpp | 290 +----------- opm/autodiff/ISTLSolverEbos.hpp | 696 +++++++++++++++++++++++++++++ 4 files changed, 700 insertions(+), 291 deletions(-) create mode 100644 opm/autodiff/ISTLSolverEbos.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 8a9e6174e..939da47be 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -297,6 +297,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/GridInit.hpp opm/autodiff/ImpesTPFAAD.hpp opm/autodiff/ISTLSolver.hpp + opm/autodiff/ISTLSolverEbos.hpp opm/autodiff/IterationReport.hpp opm/autodiff/moduleVersion.hpp opm/autodiff/multiPhaseUpwind.hpp diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index fbb466abc..f0ddaf57a 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -47,7 +47,7 @@ #include #include -#include +#include #include #include @@ -124,7 +124,7 @@ namespace Opm { typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector BVector; - typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType; + typedef ISTLSolverEbos ISTLSolverType; //typedef typename SolutionVector :: value_type PrimaryVariables ; // --------- Public methods --------- diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index 7692dbbdd..2948c7054 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -20,6 +20,7 @@ #ifndef OPM_ISTLSOLVER_HEADER_INCLUDED #define OPM_ISTLSOLVER_HEADER_INCLUDED +#include #include #include #include @@ -41,297 +42,8 @@ #include -namespace Dune -{ -namespace FMatrixHelp { -//! invert 4x4 Matrix without changing the original matrix -template -static inline K invertMatrix (const FieldMatrix &matrix, FieldMatrix &inverse) -{ - inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] - - matrix[1][1] * matrix[2][3] * matrix[3][2] - - matrix[2][1] * matrix[1][2] * matrix[3][3] + - matrix[2][1] * matrix[1][3] * matrix[3][2] + - matrix[3][1] * matrix[1][2] * matrix[2][3] - - matrix[3][1] * matrix[1][3] * matrix[2][2]; - - inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] + - matrix[1][0] * matrix[2][3] * matrix[3][2] + - matrix[2][0] * matrix[1][2] * matrix[3][3] - - matrix[2][0] * matrix[1][3] * matrix[3][2] - - matrix[3][0] * matrix[1][2] * matrix[2][3] + - matrix[3][0] * matrix[1][3] * matrix[2][2]; - - inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] - - matrix[1][0] * matrix[2][3] * matrix[3][1] - - matrix[2][0] * matrix[1][1] * matrix[3][3] + - matrix[2][0] * matrix[1][3] * matrix[3][1] + - matrix[3][0] * matrix[1][1] * matrix[2][3] - - matrix[3][0] * matrix[1][3] * matrix[2][1]; - - inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] + - matrix[1][0] * matrix[2][2] * matrix[3][1] + - matrix[2][0] * matrix[1][1] * matrix[3][2] - - matrix[2][0] * matrix[1][2] * matrix[3][1] - - matrix[3][0] * matrix[1][1] * matrix[2][2] + - matrix[3][0] * matrix[1][2] * matrix[2][1]; - - inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] + - matrix[0][1] * matrix[2][3] * matrix[3][2] + - matrix[2][1] * matrix[0][2] * matrix[3][3] - - matrix[2][1] * matrix[0][3] * matrix[3][2] - - matrix[3][1] * matrix[0][2] * matrix[2][3] + - matrix[3][1] * matrix[0][3] * matrix[2][2]; - - inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] - - matrix[0][0] * matrix[2][3] * matrix[3][2] - - matrix[2][0] * matrix[0][2] * matrix[3][3] + - matrix[2][0] * matrix[0][3] * matrix[3][2] + - matrix[3][0] * matrix[0][2] * matrix[2][3] - - matrix[3][0] * matrix[0][3] * matrix[2][2]; - - inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] + - matrix[0][0] * matrix[2][3] * matrix[3][1] + - matrix[2][0] * matrix[0][1] * matrix[3][3] - - matrix[2][0] * matrix[0][3] * matrix[3][1] - - matrix[3][0] * matrix[0][1] * matrix[2][3] + - matrix[3][0] * matrix[0][3] * matrix[2][1]; - - inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] - - matrix[0][0] * matrix[2][2] * matrix[3][1] - - matrix[2][0] * matrix[0][1] * matrix[3][2] + - matrix[2][0] * matrix[0][2] * matrix[3][1] + - matrix[3][0] * matrix[0][1] * matrix[2][2] - - matrix[3][0] * matrix[0][2] * matrix[2][1]; - - inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] - - matrix[0][1] * matrix[1][3] * matrix[3][2] - - matrix[1][1] * matrix[0][2] * matrix[3][3] + - matrix[1][1] * matrix[0][3] * matrix[3][2] + - matrix[3][1] * matrix[0][2] * matrix[1][3] - - matrix[3][1] * matrix[0][3] * matrix[1][2]; - - inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] + - matrix[0][0] * matrix[1][3] * matrix[3][2] + - matrix[1][0] * matrix[0][2] * matrix[3][3] - - matrix[1][0] * matrix[0][3] * matrix[3][2] - - matrix[3][0] * matrix[0][2] * matrix[1][3] + - matrix[3][0] * matrix[0][3] * matrix[1][2]; - - inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] - - matrix[0][0] * matrix[1][3] * matrix[3][1] - - matrix[1][0] * matrix[0][1] * matrix[3][3] + - matrix[1][0] * matrix[0][3] * matrix[3][1] + - matrix[3][0] * matrix[0][1] * matrix[1][3] - - matrix[3][0] * matrix[0][3] * matrix[1][1]; - - inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] + - matrix[0][0] * matrix[1][2] * matrix[3][1] + - matrix[1][0] * matrix[0][1] * matrix[3][2] - - matrix[1][0] * matrix[0][2] * matrix[3][1] - - matrix[3][0] * matrix[0][1] * matrix[1][2] + - matrix[3][0] * matrix[0][2] * matrix[1][1]; - - inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] + - matrix[0][1] * matrix[1][3] * matrix[2][2] + - matrix[1][1] * matrix[0][2] * matrix[2][3] - - matrix[1][1] * matrix[0][3] * matrix[2][2] - - matrix[2][1] * matrix[0][2] * matrix[1][3] + - matrix[2][1] * matrix[0][3] * matrix[1][2]; - - inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] - - matrix[0][0] * matrix[1][3] * matrix[2][2] - - matrix[1][0] * matrix[0][2] * matrix[2][3] + - matrix[1][0] * matrix[0][3] * matrix[2][2] + - matrix[2][0] * matrix[0][2] * matrix[1][3] - - matrix[2][0] * matrix[0][3] * matrix[1][2]; - - inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] + - matrix[0][0] * matrix[1][3] * matrix[2][1] + - matrix[1][0] * matrix[0][1] * matrix[2][3] - - matrix[1][0] * matrix[0][3] * matrix[2][1] - - matrix[2][0] * matrix[0][1] * matrix[1][3] + - matrix[2][0] * matrix[0][3] * matrix[1][1]; - - inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] - - matrix[0][0] * matrix[1][2] * matrix[2][1] - - matrix[1][0] * matrix[0][1] * matrix[2][2] + - matrix[1][0] * matrix[0][2] * matrix[2][1] + - matrix[2][0] * matrix[0][1] * matrix[1][2] - - matrix[2][0] * matrix[0][2] * matrix[1][1]; - - K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0]; - - // return identity for singular or nearly singular matrices. - if (std::abs(det) < 1e-40) { - for (int i = 0; i < 4; ++i){ - inverse[i][i] = 1.0; - } - return 1.0; - } - K inv_det = 1.0 / det; - inverse *= inv_det; - - return det; -} -} // end FMatrixHelp - -namespace ISTLUtility { - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling matrix.invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - Dune::FMatrixPrecision::set_singular_limit(1.e-20); - matrix.invert(); -} - -} // end ISTLUtility - -template -class MatrixBlock : public Dune::FieldMatrix -{ -public: - typedef Dune::FieldMatrix BaseType; - - using BaseType :: operator= ; - using BaseType :: rows; - using BaseType :: cols; - explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {} - void invert() - { - ISTLUtility::invertMatrix( *this ); - } - const BaseType& asBase() const { return static_cast< const BaseType& > (*this); } - BaseType& asBase() { return static_cast< BaseType& > (*this); } -}; - -template -void -print_row (std::ostream& s, const MatrixBlock& A, - typename FieldMatrix::size_type I, - typename FieldMatrix::size_type J, - typename FieldMatrix::size_type therow, int width, - int precision) -{ - print_row(s, A.asBase(), I, J, therow, width, precision); -} - -template -K& firstmatrixelement (MatrixBlock& A) -{ - return firstmatrixelement( A.asBase() ); -} - - - -template -struct MatrixDimension< MatrixBlock< Scalar, n, m > > -: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType > -{ -}; - - -#if HAVE_UMFPACK - -/// \brief UMFPack specialization for MatrixBlock to make AMG happy -/// -/// Without this the empty default implementation would be used. -template -class UMFPack, A> > - : public UMFPack, A> > -{ - typedef UMFPack, A> > Base; - typedef BCRSMatrix, A> Matrix; - -public: - typedef BCRSMatrix, A> RealMatrix; - - UMFPack(const RealMatrix& matrix, int verbose, bool) - : Base(reinterpret_cast(matrix), verbose) - {} -}; -#endif - -#if HAVE_SUPERLU - -/// \brief SuperLU specialization for MatrixBlock to make AMG happy -/// -/// Without this the empty default implementation would be used. -template -class SuperLU, A> > - : public SuperLU, A> > -{ - typedef SuperLU, A> > Base; - typedef BCRSMatrix, A> Matrix; - -public: - typedef BCRSMatrix, A> RealMatrix; - - SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true) - : Base(reinterpret_cast(matrix), verbose, reuse) - {} -}; -#endif - - -} // end namespace Dune - namespace Opm { -namespace Detail -{ - //! calculates ret = A^T * B - template< class K, int m, int n, int p > - static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, n, m > &A, - const Dune::FieldMatrix< K, n, p > &B, - Dune::FieldMatrix< K, m, p > &ret ) - { - typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type; - - for( size_type i = 0; i < m; ++i ) - { - for( size_type j = 0; j < p; ++j ) - { - ret[ i ][ j ] = K( 0 ); - for( size_type k = 0; k < n; ++k ) - ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ]; - } - } - } -} /// This class solves the fully implicit black-oil system by /// solving the reduced system (after eliminating well variables) /// as a block-structured matrix (one block for all cell variables) for a fixed diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp new file mode 100644 index 000000000..f09a6030c --- /dev/null +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -0,0 +1,696 @@ +/* + Copyright 2016 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED +#define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(GlobalEqVector); +NEW_PROP_TAG(JacobianMatrix); +NEW_PROP_TAG(Indices); + +END_PROPERTIES + +namespace Dune +{ +namespace FMatrixHelp { +//! invert 4x4 Matrix without changing the original matrix +template +static inline K invertMatrix (const FieldMatrix &matrix, FieldMatrix &inverse) +{ + inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] - + matrix[1][1] * matrix[2][3] * matrix[3][2] - + matrix[2][1] * matrix[1][2] * matrix[3][3] + + matrix[2][1] * matrix[1][3] * matrix[3][2] + + matrix[3][1] * matrix[1][2] * matrix[2][3] - + matrix[3][1] * matrix[1][3] * matrix[2][2]; + + inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] + + matrix[1][0] * matrix[2][3] * matrix[3][2] + + matrix[2][0] * matrix[1][2] * matrix[3][3] - + matrix[2][0] * matrix[1][3] * matrix[3][2] - + matrix[3][0] * matrix[1][2] * matrix[2][3] + + matrix[3][0] * matrix[1][3] * matrix[2][2]; + + inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] - + matrix[1][0] * matrix[2][3] * matrix[3][1] - + matrix[2][0] * matrix[1][1] * matrix[3][3] + + matrix[2][0] * matrix[1][3] * matrix[3][1] + + matrix[3][0] * matrix[1][1] * matrix[2][3] - + matrix[3][0] * matrix[1][3] * matrix[2][1]; + + inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] + + matrix[1][0] * matrix[2][2] * matrix[3][1] + + matrix[2][0] * matrix[1][1] * matrix[3][2] - + matrix[2][0] * matrix[1][2] * matrix[3][1] - + matrix[3][0] * matrix[1][1] * matrix[2][2] + + matrix[3][0] * matrix[1][2] * matrix[2][1]; + + inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] + + matrix[0][1] * matrix[2][3] * matrix[3][2] + + matrix[2][1] * matrix[0][2] * matrix[3][3] - + matrix[2][1] * matrix[0][3] * matrix[3][2] - + matrix[3][1] * matrix[0][2] * matrix[2][3] + + matrix[3][1] * matrix[0][3] * matrix[2][2]; + + inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] - + matrix[0][0] * matrix[2][3] * matrix[3][2] - + matrix[2][0] * matrix[0][2] * matrix[3][3] + + matrix[2][0] * matrix[0][3] * matrix[3][2] + + matrix[3][0] * matrix[0][2] * matrix[2][3] - + matrix[3][0] * matrix[0][3] * matrix[2][2]; + + inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] + + matrix[0][0] * matrix[2][3] * matrix[3][1] + + matrix[2][0] * matrix[0][1] * matrix[3][3] - + matrix[2][0] * matrix[0][3] * matrix[3][1] - + matrix[3][0] * matrix[0][1] * matrix[2][3] + + matrix[3][0] * matrix[0][3] * matrix[2][1]; + + inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] - + matrix[0][0] * matrix[2][2] * matrix[3][1] - + matrix[2][0] * matrix[0][1] * matrix[3][2] + + matrix[2][0] * matrix[0][2] * matrix[3][1] + + matrix[3][0] * matrix[0][1] * matrix[2][2] - + matrix[3][0] * matrix[0][2] * matrix[2][1]; + + inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] - + matrix[0][1] * matrix[1][3] * matrix[3][2] - + matrix[1][1] * matrix[0][2] * matrix[3][3] + + matrix[1][1] * matrix[0][3] * matrix[3][2] + + matrix[3][1] * matrix[0][2] * matrix[1][3] - + matrix[3][1] * matrix[0][3] * matrix[1][2]; + + inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] + + matrix[0][0] * matrix[1][3] * matrix[3][2] + + matrix[1][0] * matrix[0][2] * matrix[3][3] - + matrix[1][0] * matrix[0][3] * matrix[3][2] - + matrix[3][0] * matrix[0][2] * matrix[1][3] + + matrix[3][0] * matrix[0][3] * matrix[1][2]; + + inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] - + matrix[0][0] * matrix[1][3] * matrix[3][1] - + matrix[1][0] * matrix[0][1] * matrix[3][3] + + matrix[1][0] * matrix[0][3] * matrix[3][1] + + matrix[3][0] * matrix[0][1] * matrix[1][3] - + matrix[3][0] * matrix[0][3] * matrix[1][1]; + + inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] + + matrix[0][0] * matrix[1][2] * matrix[3][1] + + matrix[1][0] * matrix[0][1] * matrix[3][2] - + matrix[1][0] * matrix[0][2] * matrix[3][1] - + matrix[3][0] * matrix[0][1] * matrix[1][2] + + matrix[3][0] * matrix[0][2] * matrix[1][1]; + + inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] + + matrix[0][1] * matrix[1][3] * matrix[2][2] + + matrix[1][1] * matrix[0][2] * matrix[2][3] - + matrix[1][1] * matrix[0][3] * matrix[2][2] - + matrix[2][1] * matrix[0][2] * matrix[1][3] + + matrix[2][1] * matrix[0][3] * matrix[1][2]; + + inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] - + matrix[0][0] * matrix[1][3] * matrix[2][2] - + matrix[1][0] * matrix[0][2] * matrix[2][3] + + matrix[1][0] * matrix[0][3] * matrix[2][2] + + matrix[2][0] * matrix[0][2] * matrix[1][3] - + matrix[2][0] * matrix[0][3] * matrix[1][2]; + + inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] + + matrix[0][0] * matrix[1][3] * matrix[2][1] + + matrix[1][0] * matrix[0][1] * matrix[2][3] - + matrix[1][0] * matrix[0][3] * matrix[2][1] - + matrix[2][0] * matrix[0][1] * matrix[1][3] + + matrix[2][0] * matrix[0][3] * matrix[1][1]; + + inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] - + matrix[0][0] * matrix[1][2] * matrix[2][1] - + matrix[1][0] * matrix[0][1] * matrix[2][2] + + matrix[1][0] * matrix[0][2] * matrix[2][1] + + matrix[2][0] * matrix[0][1] * matrix[1][2] - + matrix[2][0] * matrix[0][2] * matrix[1][1]; + + K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0]; + + // return identity for singular or nearly singular matrices. + if (std::abs(det) < 1e-40) { + for (int i = 0; i < 4; ++i){ + inverse[i][i] = 1.0; + } + return 1.0; + } + K inv_det = 1.0 / det; + inverse *= inv_det; + + return det; +} +} // end FMatrixHelp + +namespace ISTLUtility { + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling matrix.invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + Dune::FMatrixPrecision::set_singular_limit(1.e-20); + matrix.invert(); +} + +} // end ISTLUtility + +template +class MatrixBlock : public Dune::FieldMatrix +{ +public: + typedef Dune::FieldMatrix BaseType; + + using BaseType :: operator= ; + using BaseType :: rows; + using BaseType :: cols; + explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {} + void invert() + { + ISTLUtility::invertMatrix( *this ); + } + const BaseType& asBase() const { return static_cast< const BaseType& > (*this); } + BaseType& asBase() { return static_cast< BaseType& > (*this); } +}; + +template +void +print_row (std::ostream& s, const MatrixBlock& A, + typename FieldMatrix::size_type I, + typename FieldMatrix::size_type J, + typename FieldMatrix::size_type therow, int width, + int precision) +{ + print_row(s, A.asBase(), I, J, therow, width, precision); +} + +template +K& firstmatrixelement (MatrixBlock& A) +{ + return firstmatrixelement( A.asBase() ); +} + + + +template +struct MatrixDimension< MatrixBlock< Scalar, n, m > > +: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType > +{ +}; + + +#if HAVE_UMFPACK + +/// \brief UMFPack specialization for MatrixBlock to make AMG happy +/// +/// Without this the empty default implementation would be used. +template +class UMFPack, A> > + : public UMFPack, A> > +{ + typedef UMFPack, A> > Base; + typedef BCRSMatrix, A> Matrix; + +public: + typedef BCRSMatrix, A> RealMatrix; + + UMFPack(const RealMatrix& matrix, int verbose, bool) + : Base(reinterpret_cast(matrix), verbose) + {} +}; +#endif + +#if HAVE_SUPERLU + +/// \brief SuperLU specialization for MatrixBlock to make AMG happy +/// +/// Without this the empty default implementation would be used. +template +class SuperLU, A> > + : public SuperLU, A> > +{ + typedef SuperLU, A> > Base; + typedef BCRSMatrix, A> Matrix; + +public: + typedef BCRSMatrix, A> RealMatrix; + + SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true) + : Base(reinterpret_cast(matrix), verbose, reuse) + {} +}; +#endif + + +} // end namespace Dune + +namespace Opm +{ +namespace Detail +{ + //! calculates ret = A^T * B + template< class K, int m, int n, int p > + static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, n, m > &A, + const Dune::FieldMatrix< K, n, p > &B, + Dune::FieldMatrix< K, m, p > &ret ) + { + typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type; + + for( size_type i = 0; i < m; ++i ) + { + for( size_type j = 0; j < p; ++j ) + { + ret[ i ][ j ] = K( 0 ); + for( size_type k = 0; k < n; ++k ) + ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ]; + } + } + } +} + /// This class solves the fully implicit black-oil system by + /// solving the reduced system (after eliminating well variables) + /// as a block-structured matrix (one block for all cell variables) for a fixed + /// number of cell variables np . + /// \tparam MatrixBlockType The type of the matrix block used. + /// \tparam VectorBlockType The type of the vector block used. + /// \tparam pressureIndex The index of the pressure component in the vector + /// vector block. It is used to guide the AMG coarsening. + /// Default is zero. + template + class ISTLSolverEbos : public NewtonIterationBlackoilInterface + { + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + enum { pressureIndex = Indices::pressureSwitchIdx }; + + public: + typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; + + typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector; + /// Construct a system solver. + /// \param[in] param parameters controlling the behaviour of the linear solvers + /// \param[in] parallelInformation In the case of a parallel run + /// with dune-istl the information about the parallelization. + ISTLSolverEbos(const NewtonIterationBlackoilInterleavedParameters& param, + const boost::any& parallelInformation_arg=boost::any()) + : iterations_( 0 ), + parallelInformation_(parallelInformation_arg), + isIORank_(isIORank(parallelInformation_arg)), + parameters_( param ) + { + } + + /// Construct a system solver. + /// \param[in] param ParameterGroup controlling the behaviour of the linear solvers + /// \param[in] parallelInformation In the case of a parallel run + /// with dune-istl the information about the parallelization. + ISTLSolverEbos(const ParameterGroup& param, + const boost::any& parallelInformation_arg=boost::any()) + : iterations_( 0 ), + parallelInformation_(parallelInformation_arg), + isIORank_(isIORank(parallelInformation_arg)), + parameters_( param ) + { + } + + // dummy method that is not implemented for this class + SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual&) const + { + OPM_THROW(std::logic_error,"This method is not implemented"); + return SolutionVector(); + } + + /// Solve the system of linear equations Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] residual residual object containing A and b. + /// \return the solution x + + /// \copydoc NewtonIterationBlackoilInterface::iterations + int iterations () const { return iterations_; } + + /// \copydoc NewtonIterationBlackoilInterface::parallelInformation + const boost::any& parallelInformation() const { return parallelInformation_; } + + public: + /// \brief construct the CPR preconditioner and the solver. + /// \tparam P The type of the parallel information. + /// \param parallelInformation the information about the parallelization. +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + template +#else + template +#endif + void constructPreconditionerAndSolve(LinearOperator& linearOperator, + Vector& x, Vector& istlb, + const POrComm& parallelInformation_arg, + Dune::InverseOperatorResult& result) const + { + // Construct scalar product. +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + auto sp = Dune::createScalarProduct(parallelInformation_arg, category); +#else + typedef Dune::ScalarProductChooser ScalarProductChooser; + typedef std::unique_ptr SPPointer; + SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg)); +#endif + + // Communicate if parallel. + parallelInformation_arg.copyOwnerToAll(istlb, istlb); + +#if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available + if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_) + { + typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; + typedef typename CPRSelectorType::Operator MatrixOperator; + + std::unique_ptr< MatrixOperator > opA; + + if( ! std::is_same< LinearOperator, MatrixOperator > :: value ) + { + // create new operator in case linear operator and matrix operator differ + opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) ); + } + + const double relax = parameters_.ilu_relaxation_; + if ( parameters_.use_cpr_ ) + { + using Matrix = typename MatrixOperator::matrix_type; + using CouplingMetric = Dune::Amg::Diagonal; + using CritBase = Dune::Amg::SymmetricCriterion; + using Criterion = Dune::Amg::CoarsenCriterion; + using AMG = typename ISTLUtility + ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG; + + std::unique_ptr< AMG > amg; + // Construct preconditioner. + Criterion crit(15, 2000); + constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax ); + + // Solve. + solve(linearOperator, x, istlb, *sp, *amg, result); + } + else + { + typedef typename CPRSelectorType::AMG AMG; + std::unique_ptr< AMG > amg; + + // Construct preconditioner. + constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax ); + + // Solve. + solve(linearOperator, x, istlb, *sp, *amg, result); + } + } + else +#endif + { + // Construct preconditioner. + auto precond = constructPrecond(linearOperator, parallelInformation_arg); + + // Solve. + solve(linearOperator, x, istlb, *sp, *precond, result); + } + } + + + // 3x3 matrix block inversion was unstable at least 2.3 until and including + // 2.5.0. There may still be some issue with the 4x4 matrix block inversion + // we therefore still use the block inversion in OPM + typedef ParallelOverlappingILU0 >, + Vector, Vector> SeqPreconditioner; + + + template + std::unique_ptr constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const + { + const double relax = parameters_.ilu_relaxation_; + const int ilu_fillin = parameters_.ilu_fillin_level_; + std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax)); + return precond; + } + +#if HAVE_MPI + typedef Dune::OwnerOverlapCopyCommunication Comm; +#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1) + // 3x3 matrix block inversion was unstable from at least 2.3 until and + // including 2.5.0 + typedef ParallelOverlappingILU0 ParPreconditioner; +#else + typedef ParallelOverlappingILU0 >, + Vector, Vector, Comm> ParPreconditioner; +#endif + template + std::unique_ptr + constructPrecond(Operator& opA, const Comm& comm) const + { + typedef std::unique_ptr Pointer; + const double relax = parameters_.ilu_relaxation_; + return Pointer(new ParPreconditioner(opA.getmat(), comm, relax)); + } +#endif + + template + void + constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const + { + ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg ); + } + + + template + void + constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const + { + ISTLUtility::template createAMGPreconditionerPointer( opA, relax, comm, amg ); + } + + template + void + constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const + { + ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg, parameters_ ); + } + + + template + void + constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const + { + ISTLUtility::template createAMGPreconditionerPointer( opA, relax, comm, amg, parameters_ ); + } + /// \brief Solve the system using the given preconditioner and scalar product. + template + void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const + { + // TODO: Revise when linear solvers interface opm-core is done + // Construct linear solver. + // GMRes solver + int verbosity = ( isIORank_ ) ? parameters_.linear_solver_verbosity_ : 0; + + if ( parameters_.newton_use_gmres_ ) { + Dune::RestartedGMResSolver linsolve(opA, sp, precond, + parameters_.linear_solver_reduction_, + parameters_.linear_solver_restart_, + parameters_.linear_solver_maxiter_, + verbosity); + // Solve system. + linsolve.apply(x, istlb, result); + } + else { // BiCGstab solver + Dune::BiCGSTABSolver linsolve(opA, sp, precond, + parameters_.linear_solver_reduction_, + parameters_.linear_solver_maxiter_, + verbosity); + // Solve system. + linsolve.apply(x, istlb, result); + } + } + + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] A matrix A + /// \param[inout] x solution to be computed x + /// \param[in] b right hand side b + void solve(Matrix& A, Vector& x, Vector& b ) const + { + // Parallel version is deactivated until we figure out how to do it properly. +#if HAVE_MPI + if (parallelInformation_.type() == typeid(ParallelISTLInformation)) + { + typedef Dune::OwnerOverlapCopyCommunication Comm; + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation_); + Comm istlComm(info.communicator()); + + // Construct operator, scalar product and vectors needed. + typedef Dune::OverlappingSchwarzOperator Operator; + Operator opA(A, istlComm); + solve( opA, x, b, istlComm ); + } + else +#endif + { + // Construct operator, scalar product and vectors needed. + Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A ); + solve( opA, x, b ); + } + } + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] A matrix A + /// \param[inout] x solution to be computed x + /// \param[in] b right hand side b + template + void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const + { + Dune::InverseOperatorResult result; + // Parallel version is deactivated until we figure out how to do it properly. +#if HAVE_MPI + if (parallelInformation_.type() == typeid(ParallelISTLInformation)) + { + const size_t size = opA.getmat().N(); + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation_); + + // As we use a dune-istl with block size np the number of components + // per parallel is only one. + info.copyValuesTo(comm.indexSet(), comm.remoteIndices(), + size, 1); + // Construct operator, scalar product and vectors needed. + constructPreconditionerAndSolve(opA, x, b, comm, result); + } + else +#endif + { + OPM_THROW(std::logic_error,"this method if for parallel solve only"); + } + + checkConvergence( result ); + } + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] A matrix A + /// \param[inout] x solution to be computed x + /// \param[in] b right hand side b + template + void solve(Operator& opA, Vector& x, Vector& b ) const + { + Dune::InverseOperatorResult result; + // Construct operator, scalar product and vectors needed. + Dune::Amg::SequentialInformation info; + constructPreconditionerAndSolve(opA, x, b, info, result); + checkConvergence( result ); + } + + void checkConvergence( const Dune::InverseOperatorResult& result ) const + { + // store number of iterations + iterations_ = result.iterations; + + // Check for failure of linear solver. + if (!parameters_.ignoreConvergenceFailure_ && !result.converged) { + const std::string msg("Convergence failure for linear solver."); + OPM_THROW_NOLOG(LinearSolverProblem, msg); + } + } + protected: + mutable int iterations_; + boost::any parallelInformation_; + bool isIORank_; + + NewtonIterationBlackoilInterleavedParameters parameters_; + }; // end ISTLSolver + +} // namespace Opm +#endif