Merge pull request #1617 from andlaus/sparse_matrix_abstraction

Sparse matrix abstraction

opm/autodiff/BlackoilAmg.hpp

@@ -19,6 +19,7 @@
 #ifndef OPM_AMG_HEADER_INCLUDED
 #define OPM_AMG_HEADER_INCLUDED
 
+#include <ewoms/linear/matrixblock.hh>
 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <dune/istl/paamg/twolevelmethod.hh>
@@ -168,6 +169,12 @@ struct ScalarType<Dune::MatrixBlock<FieldType, ROWS, COLS> >
     typedef Dune::MatrixBlock<FieldType, 1, 1> value;
 };
 
+template<typename FieldType, int ROWS, int COLS>
+struct ScalarType<Ewoms::MatrixBlock<FieldType, ROWS, COLS> >
+{
+    typedef Ewoms::MatrixBlock<FieldType, 1, 1> value;
+};
+
 template<typename BlockType, typename Allocator>
 struct ScalarType<Dune::BCRSMatrix<BlockType, Allocator> >
 {

opm/autodiff/BlackoilModelEbos.hpp

@@ -109,6 +109,7 @@ namespace Opm {
         typedef typename GET_PROP_TYPE(TypeTag, Simulator)         Simulator;
         typedef typename GET_PROP_TYPE(TypeTag, Grid)              Grid;
         typedef typename GET_PROP_TYPE(TypeTag, ElementContext)    ElementContext;
+        typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
         typedef typename GET_PROP_TYPE(TypeTag, SolutionVector)    SolutionVector ;
         typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables)  PrimaryVariables ;
         typedef typename GET_PROP_TYPE(TypeTag, FluidSystem)       FluidSystem;
@@ -126,8 +127,8 @@ namespace Opm {
         static const int temperatureIdx = Indices::temperatureIdx;
 
         typedef Dune::FieldVector<Scalar, numEq >        VectorBlockType;
-        typedef Dune::FieldMatrix<Scalar, numEq, numEq >        MatrixBlockType;
+        typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
-        typedef Dune::BCRSMatrix <MatrixBlockType>      Mat;
+        typedef typename SparseMatrixAdapter::IstlMatrix Mat;
         typedef Dune::BlockVector<VectorBlockType>      BVector;
 
         typedef ISTLSolverEbos<TypeTag> ISTLSolverType;
@@ -358,13 +359,13 @@ namespace Opm {
             ebosSimulator_.model().linearizer().linearize();
             ebosSimulator_.problem().endIteration();
 
-            auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
+            auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
             if (param_.matrix_add_well_contributions_) {
-                wellModel().addWellContributions(ebosJac);
+                wellModel().addWellContributions(ebosJac.istlMatrix());
             }
             if ( param_.preconditioner_add_well_contributions_ &&
                                   ! param_.matrix_add_well_contributions_ ) {
-                matrix_for_preconditioner_ .reset(new Mat(ebosJac));
+                matrix_for_preconditioner_ .reset(new Mat(ebosJac.istlMatrix()));
                 wellModel().addWellContributions(*matrix_for_preconditioner_);
             }
 
@@ -494,7 +495,7 @@ namespace Opm {
             // r = [r, r_well], where r is the residual and r_well the well residual.
             // r -= B^T * D^-1 r_well
 
-            auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
+            auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
             auto& ebosResid = ebosSimulator_.model().linearizer().residual();
 
             wellModel().apply(ebosResid);
@@ -502,13 +503,13 @@ namespace Opm {
             // set initial guess
             x = 0.0;
 
-            const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac;
+            const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac.istlMatrix();
             // Solve system.
             if( isParallel() )
             {	      
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
 
-                auto ebosJacIgnoreOverlap = Mat(ebosJac);
+                auto ebosJacIgnoreOverlap = Mat(ebosJac.istlMatrix());
                 //remove ghost rows in local matrix
                 makeOverlapRowsInvalid(ebosJacIgnoreOverlap);
 
@@ -522,7 +523,7 @@ namespace Opm {
             else
             {
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
-                Operator opA(ebosJac, actual_mat_for_prec, wellModel());
+                Operator opA(ebosJac.istlMatrix(), actual_mat_for_prec, wellModel());
                 istlSolver().solve( opA, x, ebosResid );
             }
         }

opm/autodiff/BlackoilWellModel.hpp

@@ -84,7 +84,7 @@ namespace Opm {
             typedef typename GET_PROP_TYPE(TypeTag, Scalar)              Scalar;
             typedef typename GET_PROP_TYPE(TypeTag, RateVector)          RateVector;
             typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector)      GlobalEqVector;
-            typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix)      JacobianMatrix;
+            typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
 
             typedef typename Ewoms::BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
 
@@ -103,7 +103,7 @@ namespace Opm {
 #else
             typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
 #endif
-            typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
+            typedef typename SparseMatrixAdapter::IstlMatrix Mat;
 
             typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
 
@@ -127,7 +127,7 @@ namespace Opm {
             void applyInitial()
             {}
 
-            void linearize(JacobianMatrix& mat , GlobalEqVector& res);
+            void linearize(SparseMatrixAdapter& mat , GlobalEqVector& res);
 
             void postSolve(GlobalEqVector& deltaX)
             {

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -94,7 +94,7 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    linearize(JacobianMatrix& mat , GlobalEqVector& res)
+    linearize(SparseMatrixAdapter& mat , GlobalEqVector& res)
     {
         if (!localWellsActive())
             return;

opm/autodiff/ISTLSolver.hpp

@@ -20,7 +20,7 @@
 #ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
 #define OPM_ISTLSOLVER_HEADER_INCLUDED
 
-#include <opm/autodiff/ISTLSolverEbos.hpp>
+#include <opm/autodiff/MPIUtilities.hpp>
 #include <opm/autodiff/BlackoilAmg.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
@@ -28,6 +28,7 @@
 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/MatrixBlock.hpp>
 
 #include <opm/common/Exceptions.hpp>
 #include <opm/core/linalg/ParallelIstlInformation.hpp>

opm/autodiff/ISTLSolverEbos.hpp

@@ -20,6 +20,7 @@
 #ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
 #define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
 
+#include <opm/autodiff/MatrixBlock.hpp>
 #include <opm/autodiff/BlackoilAmg.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
@@ -50,304 +51,13 @@ NEW_TYPE_TAG(FlowIstlSolver, INHERITS_FROM(FlowIstlSolverParams));
 
 NEW_PROP_TAG(Scalar);
 NEW_PROP_TAG(GlobalEqVector);
-NEW_PROP_TAG(JacobianMatrix);
+NEW_PROP_TAG(SparseMatrixAdapter);
 NEW_PROP_TAG(Indices);
 
 END_PROPERTIES
 
-namespace Dune
-{
-namespace FMatrixHelp {
-//! invert 4x4 Matrix without changing the original matrix
-template <typename K>
-static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &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 <typename K>
-static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
-{
-    FieldMatrix<K,1,1> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling FMatrixHelp::invert
-template <typename K>
-static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
-{
-    FieldMatrix<K,2,2> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling FMatrixHelp::invert
-template <typename K>
-static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
-{
-    FieldMatrix<K,3,3> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling FMatrixHelp::invert
-template <typename K>
-static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
-{
-    FieldMatrix<K,4,4> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling matrix.invert
-template <typename K, int n>
-static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
-{
-#if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
-    Dune::FMatrixPrecision<K>::set_singular_limit(1.e-20);
-#endif
-    matrix.invert();
-}
-
-} // end ISTLUtility
-
-template <class Scalar, int n, int m>
-class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
-{
-public:
-    typedef Dune::FieldMatrix<Scalar, n, m>  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<class K, int n, int m>
-void
-print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
-           typename FieldMatrix<K,n,m>::size_type I,
-           typename FieldMatrix<K,n,m>::size_type J,
-           typename FieldMatrix<K,n,m>::size_type therow, int width,
-           int precision)
-{
-    print_row(s, A.asBase(), I, J, therow, width, precision);
-}
-
-template<class K, int n, int m>
-K& firstmatrixelement (MatrixBlock<K,n,m>& A)
-{
-   return firstmatrixelement( A.asBase() );
-}
-
-
-
-template<typename Scalar, int n, int m>
-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<typename T, typename A, int n, int m>
-class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
-    : public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
-{
-    typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
-    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
-
-public:
-    typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
-
-    UMFPack(const RealMatrix& matrix, int verbose, bool)
-        : Base(reinterpret_cast<const Matrix&>(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<typename T, typename A, int n, int m>
-class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
-    : public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
-{
-    typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
-    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
-
-public:
-    typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
-
-    SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
-        : Base(reinterpret_cast<const Matrix&>(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
@@ -361,10 +71,12 @@ namespace Detail
     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, SparseMatrixAdapter) SparseMatrixAdapter;
         typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
         typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
 
+        typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
+
         enum { pressureIndex = Indices::pressureSwitchIdx };
 
     public:

opm/autodiff/MPIUtilities.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_MPIUTILITIES_HEADER_INCLUDED
 #define OPM_MPIUTILITIES_HEADER_INCLUDED
 
+#include <boost/any.hpp>
+
 namespace Opm
 {
     /// Return true if this is a serial run, or rank zero on an MPI run.

opm/autodiff/MatrixBlock.hpp

@@ -0,0 +1,321 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_MATRIX_BLOCK_HEADER_INCLUDED
+#define OPM_MATRIX_BLOCK_HEADER_INCLUDED
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+namespace Dune
+{
+namespace FMatrixHelp {
+//! invert 4x4 Matrix without changing the original matrix
+template <typename K>
+static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &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 <typename K>
+static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
+{
+    FieldMatrix<K,1,1> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling FMatrixHelp::invert
+template <typename K>
+static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
+{
+    FieldMatrix<K,2,2> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling FMatrixHelp::invert
+template <typename K>
+static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
+{
+    FieldMatrix<K,3,3> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling FMatrixHelp::invert
+template <typename K>
+static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
+{
+    FieldMatrix<K,4,4> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling matrix.invert
+template <typename K, int n>
+static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
+{
+#if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
+    Dune::FMatrixPrecision<K>::set_singular_limit(1.e-20);
+#endif
+    matrix.invert();
+}
+
+} // end ISTLUtility
+
+template <class Scalar, int n, int m>
+class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
+{
+public:
+    typedef Dune::FieldMatrix<Scalar, n, m>  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<class K, int n, int m>
+void
+print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
+           typename FieldMatrix<K,n,m>::size_type I,
+           typename FieldMatrix<K,n,m>::size_type J,
+           typename FieldMatrix<K,n,m>::size_type therow, int width,
+           int precision)
+{
+    print_row(s, A.asBase(), I, J, therow, width, precision);
+}
+
+template<class K, int n, int m>
+K& firstmatrixelement (MatrixBlock<K,n,m>& A)
+{
+   return firstmatrixelement( A.asBase() );
+}
+
+
+
+template<typename Scalar, int n, int m>
+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<typename T, typename A, int n, int m>
+class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
+    : public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
+{
+    typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
+    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
+
+public:
+    typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
+
+    UMFPack(const RealMatrix& matrix, int verbose, bool)
+        : Base(reinterpret_cast<const Matrix&>(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<typename T, typename A, int n, int m>
+class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
+    : public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
+{
+    typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
+    typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
+
+public:
+    typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
+
+    SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
+        : Base(reinterpret_cast<const Matrix&>(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 ];
+            }
+        }
+    }
+} // namespace Detail
+} // namespace Opm
+
+#endif

opm/autodiff/NewtonIterationBlackoilInterleaved.cpp

@@ -26,13 +26,13 @@
 // Define making clear that the simulator supports AMG
 #define FLOW_SUPPORT_AMG !defined(HAVE_UMFPACK)
 
-#include <opm/autodiff/DuneMatrix.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
 #include <opm/autodiff/NewtonIterationUtilities.hpp>
 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/DuneMatrix.hpp>
 #include <opm/common/Exceptions.hpp>
 #include <opm/core/linalg/ParallelIstlInformation.hpp>
 

opm/autodiff/WellConnectionAuxiliaryModule.hpp

@@ -36,7 +36,7 @@ class WellConnectionAuxiliaryModule
     : public Ewoms::BaseAuxiliaryModule<TypeTag>
 {
     typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
-    typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
 
 public:
 
@@ -116,7 +116,7 @@ public:
     void applyInitial()
     {}
 
-    void linearize(JacobianMatrix& , GlobalEqVector&)
+    void linearize(SparseMatrixAdapter& , GlobalEqVector&)
     {
         // Linearization is done in StandardDenseWells
     }

opm/autodiff/WellInterface.hpp

@@ -80,6 +80,7 @@ namespace Opm
         typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
         typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
         typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+        typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
         typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
 
         static const int numEq = Indices::numEq;
@@ -87,7 +88,7 @@ namespace Opm
 
         typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
         typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
-        typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
+        typedef typename SparseMatrixAdapter::IstlMatrix Mat;
         typedef Dune::BlockVector<VectorBlockType> BVector;
         typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;