From 1a348c0d29988d144a29301e95a4aabfb7c05adb Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 24 Aug 2022 09:57:49 +0200 Subject: [PATCH] changed: get rid of duplicate MatrixBlock headers/classes this has already led to some confusion. move some of the code upstream to opm-models and remove the rest of the duplicated code. the remainder of MatrixBlock.hpp is renamed to SmallDenseMatrixUtils.hpp --- CMakeLists_files.cmake | 2 +- opm/simulators/linalg/ISTLSolverEbos.hpp | 10 +- opm/simulators/linalg/MatrixBlock.hpp | 782 ------------------ .../linalg/SmallDenseMatrixUtils.hpp | 143 ++++ .../linalg/bda/opencl/ChowPatelIlu.cpp | 1 - .../wells/MultisegmentWell_impl.hpp | 4 +- opm/simulators/wells/StandardWell_impl.hpp | 14 +- tests/test_invert.cpp | 18 +- tests/test_multmatrixtransposed.cpp | 6 +- 9 files changed, 167 insertions(+), 813 deletions(-) delete mode 100644 opm/simulators/linalg/MatrixBlock.hpp create mode 100644 opm/simulators/linalg/SmallDenseMatrixUtils.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 12d58957d..92e591d7c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -299,7 +299,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/FlowLinearSolverParameters.hpp opm/simulators/linalg/GraphColoring.hpp opm/simulators/linalg/ISTLSolverEbos.hpp - opm/simulators/linalg/MatrixBlock.hpp opm/simulators/linalg/MatrixMarketSpecializations.hpp opm/simulators/linalg/OwningBlockPreconditioner.hpp opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp @@ -311,6 +310,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/PreconditionerFactory.hpp opm/simulators/linalg/PreconditionerWithUpdate.hpp opm/simulators/linalg/PropertyTree.hpp + opm/simulators/linalg/SmallDenseMatrixUtils.hpp opm/simulators/linalg/WellOperators.hpp opm/simulators/linalg/WriteSystemMatrixHelper.hpp opm/simulators/linalg/findOverlapRowsAndColumns.hpp diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 0c9d9f867..429540e6d 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -345,14 +345,6 @@ namespace Opm const std::any& parallelInformation() const { return parallelInformation_; } protected: - // 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, Dune::Amg::SequentialInformation> SeqPreconditioner; - #if HAVE_MPI typedef Dune::OwnerOverlapCopyCommunication Comm; // 3x3 matrix block inversion was unstable from at least 2.3 until and diff --git a/opm/simulators/linalg/MatrixBlock.hpp b/opm/simulators/linalg/MatrixBlock.hpp deleted file mode 100644 index c2e35bda8..000000000 --- a/opm/simulators/linalg/MatrixBlock.hpp +++ /dev/null @@ -1,782 +0,0 @@ -/* - Copyright 2016 IRIS AS - Copyright 2019 NORCE - - 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_MATRIX_BLOCK_HEADER_INCLUDED -#define OPM_MATRIX_BLOCK_HEADER_INCLUDED - -#include -#include -#include -#include -#include -#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; -} - -template -static inline K invertMatrix(const DynamicMatrix& matrix, DynamicMatrix& inverse) -{ - // this function is only for 4 X 4 matrix - assert (matrix.rows() == 4); - - 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) -{ -#if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 ) - Dune::FMatrixPrecision::set_singular_limit(1.e-20); -#endif - matrix.invert(); -} - -//! invert matrix by calling matrix.invert -template -static inline void invertMatrix(Dune::DynamicMatrix& matrix) -{ - // for 4 X 4 matrix, using the invertMatrix() function above - // it is for temporary usage, mainly to reduce the huge burden of testing - // what algorithm should be used to invert 4 X 4 matrix will be handled - // as a seperate issue - if (matrix.rows() == 4) { - Dune::DynamicMatrix A = matrix; - FMatrixHelp::invertMatrix(A, matrix); - return; - } - -#if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 ) - Dune::FMatrixPrecision::set_singular_limit(1.e-30); -#endif - 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 * B - template< class TA, class TB, class TC, class PositiveSign > - static inline void multMatrixImpl( const TA &A, // n x m - const TB &B, // n x p - TC &ret, // m x p - const PositiveSign ) - { - typedef typename TA :: size_type size_type; - typedef typename TA :: field_type K; - assert( A.N() == B.N() ); - assert( A.M() == ret.N() ); - assert( B.M() == ret.M() ); - - const size_type n = A.N(); - const size_type m = ret.N(); - const size_type p = B.M(); - for( size_type i = 0; i < m; ++i ) - { - for( size_type j = 0; j < p; ++j ) - { - K sum = 0; - for( size_type k = 0; k < n; ++k ) - { - sum += A[ i ][ k ] * B[ k ][ j ]; - } - // set value depending on given sign - ret[ i ][ j ] = PositiveSign::value ? sum : -sum; - } - } - } - - //! calculates ret = sign * (A^T * B) - //! TA, TB, and TC are not necessarily FieldMatrix, but those should - //! follow the Dune::DenseMatrix interface. - template< class TA, class TB, class TC, class PositiveSign > - static inline void multMatrixTransposedImpl ( const TA &A, // n x m - const TB &B, // n x p - TC &ret, // m x p - const PositiveSign ) - { - typedef typename TA :: size_type size_type; - typedef typename TA :: field_type K; - assert( A.N() == B.N() ); - assert( A.M() == ret.N() ); - assert( B.M() == ret.M() ); - - const size_type n = A.N(); - const size_type m = ret.N(); - const size_type p = B.M(); - for( size_type i = 0; i < m; ++i ) - { - for( size_type j = 0; j < p; ++j ) - { - K sum = 0; - for( size_type k = 0; k < n; ++k ) - { - sum += A[ k ][ i ] * B[ k ][ j ]; - } - // set value depending on given sign - ret[ i ][ j ] = PositiveSign::value ? sum : -sum; - } - } - } - - //! calculates ret = A^T * B - template - static inline void multMatrixTransposed(const DenseMatrixA& A, - const DenseMatrixB& B, - DenseMatrixC& ret) - { - multMatrixTransposedImpl( A, B, ret, std::true_type() ); - } - - //! calculates ret = -A^T * B - template - static inline void negativeMultMatrixTransposed(const DenseMatrixA& A, - const DenseMatrixB& B, - DenseMatrixC& ret) - { - multMatrixTransposedImpl( A, B, ret, std::false_type() ); - } - - //! calculates ret = A * B - template< class K> - static inline void multMatrix(const Dune::DynamicMatrix& A, - const Dune::DynamicMatrix& B, - Dune::DynamicMatrix& ret ) - { - typedef typename Dune::DynamicMatrix :: size_type size_type; - - const size_type m = A.rows(); - const size_type n = A.cols(); - - assert(n == B.rows() ); - - const size_type p = B.cols(); - - ret.resize(m, p); - - 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[ i ][ k ] * B[ k ][ j ]; - } - } - } - - //! perform out of place matrix inversion on C-style arrays - //! must have a specified block_size - template - struct Inverter - { - template - void operator()(const K *matrix [[maybe_unused]], K *inverse [[maybe_unused]]) - { - throw std::logic_error("Not implemented"); - } - }; - - //! perform out of place matrix inversion on C-style arrays - template <> - struct Inverter<4> - { - template - void operator()(const K *matrix, K *inverse) - { - // based on Dune::FMatrixHelp::invertMatrix - inverse[0] = matrix[5] * matrix[10] * matrix[15] - - matrix[5] * matrix[11] * matrix[14] - - matrix[9] * matrix[6] * matrix[15] + - matrix[9] * matrix[7] * matrix[14] + - matrix[13] * matrix[6] * matrix[11] - - matrix[13] * matrix[7] * matrix[10]; - - inverse[4] = -matrix[4] * matrix[10] * matrix[15] + - matrix[4] * matrix[11] * matrix[14] + - matrix[8] * matrix[6] * matrix[15] - - matrix[8] * matrix[7] * matrix[14] - - matrix[12] * matrix[6] * matrix[11] + - matrix[12] * matrix[7] * matrix[10]; - - inverse[8] = matrix[4] * matrix[9] * matrix[15] - - matrix[4] * matrix[11] * matrix[13] - - matrix[8] * matrix[5] * matrix[15] + - matrix[8] * matrix[7] * matrix[13] + - matrix[12] * matrix[5] * matrix[11] - - matrix[12] * matrix[7] * matrix[9]; - - inverse[12] = -matrix[4] * matrix[9] * matrix[14] + - matrix[4] * matrix[10] * matrix[13] + - matrix[8] * matrix[5] * matrix[14] - - matrix[8] * matrix[6] * matrix[13] - - matrix[12] * matrix[5] * matrix[10] + - matrix[12] * matrix[6] * matrix[9]; - - inverse[1]= -matrix[1] * matrix[10] * matrix[15] + - matrix[1] * matrix[11] * matrix[14] + - matrix[9] * matrix[2] * matrix[15] - - matrix[9] * matrix[3] * matrix[14] - - matrix[13] * matrix[2] * matrix[11] + - matrix[13] * matrix[3] * matrix[10]; - - inverse[5] = matrix[0] * matrix[10] * matrix[15] - - matrix[0] * matrix[11] * matrix[14] - - matrix[8] * matrix[2] * matrix[15] + - matrix[8] * matrix[3] * matrix[14] + - matrix[12] * matrix[2] * matrix[11] - - matrix[12] * matrix[3] * matrix[10]; - - inverse[9] = -matrix[0] * matrix[9] * matrix[15] + - matrix[0] * matrix[11] * matrix[13] + - matrix[8] * matrix[1] * matrix[15] - - matrix[8] * matrix[3] * matrix[13] - - matrix[12] * matrix[1] * matrix[11] + - matrix[12] * matrix[3] * matrix[9]; - - inverse[13] = matrix[0] * matrix[9] * matrix[14] - - matrix[0] * matrix[10] * matrix[13] - - matrix[8] * matrix[1] * matrix[14] + - matrix[8] * matrix[2] * matrix[13] + - matrix[12] * matrix[1] * matrix[10] - - matrix[12] * matrix[2] * matrix[9]; - - inverse[2] = matrix[1] * matrix[6] * matrix[15] - - matrix[1] * matrix[7] * matrix[14] - - matrix[5] * matrix[2] * matrix[15] + - matrix[5] * matrix[3] * matrix[14] + - matrix[13] * matrix[2] * matrix[7] - - matrix[13] * matrix[3] * matrix[6]; - - inverse[6] = -matrix[0] * matrix[6] * matrix[15] + - matrix[0] * matrix[7] * matrix[14] + - matrix[4] * matrix[2] * matrix[15] - - matrix[4] * matrix[3] * matrix[14] - - matrix[12] * matrix[2] * matrix[7] + - matrix[12] * matrix[3] * matrix[6]; - - inverse[10] = matrix[0] * matrix[5] * matrix[15] - - matrix[0] * matrix[7] * matrix[13] - - matrix[4] * matrix[1] * matrix[15] + - matrix[4] * matrix[3] * matrix[13] + - matrix[12] * matrix[1] * matrix[7] - - matrix[12] * matrix[3] * matrix[5]; - - inverse[14] = -matrix[0] * matrix[5] * matrix[14] + - matrix[0] * matrix[6] * matrix[13] + - matrix[4] * matrix[1] * matrix[14] - - matrix[4] * matrix[2] * matrix[13] - - matrix[12] * matrix[1] * matrix[6] + - matrix[12] * matrix[2] * matrix[5]; - - inverse[3] = -matrix[1] * matrix[6] * matrix[11] + - matrix[1] * matrix[7] * matrix[10] + - matrix[5] * matrix[2] * matrix[11] - - matrix[5] * matrix[3] * matrix[10] - - matrix[9] * matrix[2] * matrix[7] + - matrix[9] * matrix[3] * matrix[6]; - - inverse[7] = matrix[0] * matrix[6] * matrix[11] - - matrix[0] * matrix[7] * matrix[10] - - matrix[4] * matrix[2] * matrix[11] + - matrix[4] * matrix[3] * matrix[10] + - matrix[8] * matrix[2] * matrix[7] - - matrix[8] * matrix[3] * matrix[6]; - - inverse[11] = -matrix[0] * matrix[5] * matrix[11] + - matrix[0] * matrix[7] * matrix[9] + - matrix[4] * matrix[1] * matrix[11] - - matrix[4] * matrix[3] * matrix[9] - - matrix[8] * matrix[1] * matrix[7] + - matrix[8] * matrix[3] * matrix[5]; - - inverse[15] = matrix[0] * matrix[5] * matrix[10] - - matrix[0] * matrix[6] * matrix[9] - - matrix[4] * matrix[1] * matrix[10] + - matrix[4] * matrix[2] * matrix[9] + - matrix[8] * matrix[1] * matrix[6] - - matrix[8] * matrix[2] * matrix[5]; - - K det = matrix[0] * inverse[0] + matrix[1] * inverse[4] + - matrix[2] * inverse[8] + matrix[3] * inverse[12]; - - // return identity for singular or nearly singular matrices. - if (std::abs(det) < 1e-40) { - for (int i = 0; i < 4; ++i){ - inverse[4*i + i] = 1.0; - } - } - K inv_det = 1.0 / det; - - for (unsigned int i = 0; i < 4 * 4; ++i) { - inverse[i] *= inv_det; - } - } - }; - - //! perform out of place matrix inversion on C-style arrays - template <> - struct Inverter<3> - { - template - void operator()(const K *matrix, K *inverse) - { - // code generated by maple, copied from Dune::DenseMatrix - K t4 = matrix[0] * matrix[4]; - K t6 = matrix[0] * matrix[5]; - K t8 = matrix[1] * matrix[3]; - K t10 = matrix[2] * matrix[3]; - K t12 = matrix[1] * matrix[6]; - K t14 = matrix[2] * matrix[6]; - - K det = (t4 * matrix[8] - t6 * matrix[7] - t8 * matrix[8] + - t10 * matrix[7] + t12 * matrix[5] - t14 * matrix[4]); - K t17 = 1.0 / det; - - inverse[0] = (matrix[4] * matrix[8] - matrix[5] * matrix[7]) * t17; - inverse[1] = -(matrix[1] * matrix[8] - matrix[2] * matrix[7]) * t17; - inverse[2] = (matrix[1] * matrix[5] - matrix[2] * matrix[4]) * t17; - inverse[3] = -(matrix[3] * matrix[8] - matrix[5] * matrix[6]) * t17; - inverse[4] = (matrix[0] * matrix[8] - t14) * t17; - inverse[5] = -(t6 - t10) * t17; - inverse[6] = (matrix[3] * matrix[7] - matrix[4] * matrix[6]) * t17; - inverse[7] = -(matrix[0] * matrix[7] - t12) * t17; - inverse[8] = (t4 - t8) * t17; - } - }; - - //! perform out of place matrix inversion on C-style arrays - template <> - struct Inverter<2> - { - template - void operator()(const K *matrix, K *inverse) - { - // code based on Dune::DenseMatrix - K detinv = matrix[0] * matrix[3] - matrix[1] * matrix[2]; - detinv = 1 / detinv; - inverse[0] = matrix[3] * detinv; - inverse[1] = -matrix[1] * detinv; - inverse[2] = -matrix[2] * detinv; - inverse[3] = matrix[0] * detinv; - } - }; - - //! perform out of place matrix inversion on C-style arrays - template <> - struct Inverter<1> - { - template - void operator()(const K *matrix, K *inverse) - { - inverse[0] = 1.0 / matrix[0]; - } - }; - -} // namespace Detail -} // namespace Opm - -#endif diff --git a/opm/simulators/linalg/SmallDenseMatrixUtils.hpp b/opm/simulators/linalg/SmallDenseMatrixUtils.hpp new file mode 100644 index 000000000..18909445f --- /dev/null +++ b/opm/simulators/linalg/SmallDenseMatrixUtils.hpp @@ -0,0 +1,143 @@ +/* + Copyright 2016 IRIS AS + Copyright 2019 NORCE + + 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_SMALL_DENSE_MATRIX_UTILS_HEADER_INCLUDED +#define OPM_SMALL_DENSE_MATRIX_UTILS_HEADER_INCLUDED + +#include + +namespace Opm +{ +namespace detail +{ + //! calculates ret = A * B + template< class TA, class TB, class TC, class PositiveSign > + static inline void multMatrixImpl( const TA &A, // n x m + const TB &B, // n x p + TC &ret, // m x p + const PositiveSign ) + { + using size_type = typename TA :: size_type; + using K = typename TA :: field_type; + assert( A.N() == B.N() ); + assert( A.M() == ret.N() ); + assert( B.M() == ret.M() ); + + const size_type n = A.N(); + const size_type m = ret.N(); + const size_type p = B.M(); + for( size_type i = 0; i < m; ++i ) + { + for( size_type j = 0; j < p; ++j ) + { + K sum = 0; + for( size_type k = 0; k < n; ++k ) + { + sum += A[ i ][ k ] * B[ k ][ j ]; + } + // set value depending on given sign + ret[ i ][ j ] = PositiveSign::value ? sum : -sum; + } + } + } + + //! calculates ret = sign * (A^T * B) + //! TA, TB, and TC are not necessarily FieldMatrix, but those should + //! follow the Dune::DenseMatrix interface. + template< class TA, class TB, class TC, class PositiveSign > + static inline void multMatrixTransposedImpl ( const TA &A, // n x m + const TB &B, // n x p + TC &ret, // m x p + const PositiveSign ) + { + using size_type = typename TA :: size_type; + using K = typename TA :: field_type; + assert( A.N() == B.N() ); + assert( A.M() == ret.N() ); + assert( B.M() == ret.M() ); + + const size_type n = A.N(); + const size_type m = ret.N(); + const size_type p = B.M(); + for( size_type i = 0; i < m; ++i ) + { + for( size_type j = 0; j < p; ++j ) + { + K sum = 0; + for( size_type k = 0; k < n; ++k ) + { + sum += A[ k ][ i ] * B[ k ][ j ]; + } + // set value depending on given sign + ret[ i ][ j ] = PositiveSign::value ? sum : -sum; + } + } + } + + //! calculates ret = A^T * B + template + static inline void multMatrixTransposed(const DenseMatrixA& A, + const DenseMatrixB& B, + DenseMatrixC& ret) + { + multMatrixTransposedImpl( A, B, ret, std::true_type() ); + } + + //! calculates ret = -A^T * B + template + static inline void negativeMultMatrixTransposed(const DenseMatrixA& A, + const DenseMatrixB& B, + DenseMatrixC& ret) + { + multMatrixTransposedImpl( A, B, ret, std::false_type() ); + } + + //! calculates ret = A * B + template< class K> + static inline void multMatrix(const Dune::DynamicMatrix& A, + const Dune::DynamicMatrix& B, + Dune::DynamicMatrix& ret ) + { + using size_type = typename Dune::DynamicMatrix :: size_type; + + const size_type m = A.rows(); + const size_type n = A.cols(); + + assert(n == B.rows() ); + + const size_type p = B.cols(); + + ret.resize(m, p); + + 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[ i ][ k ] * B[ k ][ j ]; + } + } + } + +} // namespace detail +} // namespace Opm + +#endif diff --git a/opm/simulators/linalg/bda/opencl/ChowPatelIlu.cpp b/opm/simulators/linalg/bda/opencl/ChowPatelIlu.cpp index c4d269429..e81d2408b 100644 --- a/opm/simulators/linalg/bda/opencl/ChowPatelIlu.cpp +++ b/opm/simulators/linalg/bda/opencl/ChowPatelIlu.cpp @@ -22,7 +22,6 @@ #include #include #include -#include #include #include diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 4f9af6d54..aaf486a5e 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -741,9 +741,9 @@ namespace Opm for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) { const auto col_index = colB.index(); OffDiagMatrixBlockWellType tmp1; - Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); + detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); typename SparseMatrixAdapter::MatrixBlock tmp2; - Detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type()); + detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type()); jacobian.addToBlock(row_index, col_index, tmp2); } } diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 3194634a1..c9dc82a5e 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -21,7 +21,7 @@ #include #include -#include +#include #include #include @@ -552,7 +552,13 @@ namespace Opm // do the local inversion of D. try { this->invDuneD_ = this->duneD_; // Not strictly need if not cpr with well contributions is used - Dune::ISTLUtility::invertMatrix(this->invDuneD_[0][0]); + detail::invertMatrix(this->invDuneD_[0][0]); + } catch (NumericalProblem&) { + // for singular matrices, use identity as the inverse + this->invDuneD_[0][0] = 0.0; + for (size_t i = 0; i < this->invDuneD_[0][0].rows(); ++i) { + this->invDuneD_[0][0][i][i] = 1.0; + } } catch( ... ) { OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger); } @@ -2171,8 +2177,8 @@ namespace Opm for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB ) { - Detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp); - Detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat); + detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp); + detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat); jacobian.addToBlock( row_index, colB.index(), tmpMat ); } } diff --git a/tests/test_invert.cpp b/tests/test_invert.cpp index 5a339968b..4127961e4 100644 --- a/tests/test_invert.cpp +++ b/tests/test_invert.cpp @@ -21,7 +21,7 @@ #define BOOST_TEST_MODULE InvertSpecializationTest #include -#include +#include void checkIdentity(Dune::FieldMatrix M) { @@ -41,7 +41,7 @@ void checkIdentity(Dune::FieldMatrix M) { BOOST_AUTO_TEST_CASE(Invert4x4) { - typedef Dune::FieldMatrix BaseType; + using BaseType = Dune::FieldMatrix; BaseType matrix; BaseType eye; BaseType inverse; @@ -56,19 +56,13 @@ BOOST_AUTO_TEST_CASE(Invert4x4) matrix[3][0] = 5; matrix[0][3] = 14; - double det = Dune::FMatrixHelp::invertMatrix(matrix, inverse); + double det = Opm::detail::invertMatrix4(matrix, inverse); BOOST_CHECK_CLOSE(4, det, 1e-14); // check matrix * inverse close to identiy checkIdentity(matrix.rightmultiply(inverse)); - // check return identity matrix if singular matrix - inverse = 0.0; - double det2 = Dune::FMatrixHelp::invertMatrix(matrix_sing, inverse); - BOOST_CHECK_CLOSE(1.0, det2, 1e-14); - checkIdentity(inverse); + // check singular matrix + BOOST_CHECK_THROW(Opm::detail::invertMatrix4(matrix_sing, inverse), + Opm::NumericalProblem); } - - - - diff --git a/tests/test_multmatrixtransposed.cpp b/tests/test_multmatrixtransposed.cpp index 8564c1362..829d6cbe0 100644 --- a/tests/test_multmatrixtransposed.cpp +++ b/tests/test_multmatrixtransposed.cpp @@ -21,10 +21,12 @@ #define BOOST_TEST_MODULE MultMatrixTransposed #include -#include +#include + +#include using namespace Dune; -using namespace Opm::Detail; +using namespace Opm::detail; BOOST_AUTO_TEST_CASE(testmultmatrixtrans) {