Add support for experimental alternative of Dune::BCRSMatrix

CMakeLists.txt

@@ -35,6 +35,7 @@ option(USE_HYPRE "Use the Hypre library for linear solvers?" OFF)
 set(OPM_COMPILE_COMPONENTS "2;3;4;5;6;7;8;9;10" CACHE STRING "The components to compile support for")
 option(USE_OPENCL "Enable OpenCL support?" ON)
 option(USE_DEV_SIMULATOR_IN_TESTS "Use the development simulator binaries in tests?" OFF)
+option(USE_EXPERIMENTAL_ISTL_MATRIX "Use experimental IBCRSMatrix from dune-istl?" OFF)
 
 # Wrapper for opm_add_target_options that also adds
 # compile definitions and target links from the library
@@ -126,6 +127,9 @@ macro(opm-simulators_prereqs_hook)
     target_link_libraries(opmsimulators PUBLIC HDF5::HDF5)
     target_compile_definitions(opmsimulators PUBLIC HAVE_HDF5=1)
   endif()
+  if(USE_EXPERIMENTAL_ISTL_MATRIX)
+    target_compile_definitions(opmsimulators PUBLIC OPM_USE_EXPERIMENTAL_IBCRSMATRIX=1)
+  endif()
 
   if (OPM_ENABLE_PYTHON)
     # Be backwards compatible.

CMakeLists_files.cmake

@@ -1095,6 +1095,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/amgcpr.hh
   opm/simulators/linalg/bicgstabsolver.hh
   opm/simulators/linalg/blacklist.hh
+  opm/simulators/linalg/BlockSparseMatrix.hpp
   opm/simulators/linalg/combinedcriterion.hh
   opm/simulators/linalg/convergencecriterion.hh
   opm/simulators/linalg/DILU.hpp

flowexperimental/comp/wells/CompWellEquations.hpp

@@ -20,9 +20,10 @@
 #ifndef OPM_COMP_WELL_EQUATIONS_HPP
 #define OPM_COMP_WELL_EQUATIONS_HPP
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 namespace Opm {
@@ -48,11 +49,11 @@ class CompWellEquations
 
     // the matrix type for the diagonal matrix D
     using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq>;
-    using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
+    using DiagMatWell = BlockSparseMatrix<DiagMatrixBlockWellType>;
 
     // the matrix type for the non-diagonal matrix B and C^T
     using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, numEq>;
-    using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
+    using OffDiagMatWell = BlockSparseMatrix<OffDiagMatrixBlockWellType>;
 
     void init(const int num_conn, const std::vector<std::size_t>& cells);
 

opm/simulators/flow/GenericTemperatureModel.hpp

@@ -28,7 +28,6 @@
 #ifndef OPM_GENERIC_TEMPERATURE_MODEL_HPP
 #define OPM_GENERIC_TEMPERATURE_MODEL_HPP
 
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/grid/common/CartesianIndexMapper.hpp>
 
@@ -37,6 +36,7 @@
 #include <opm/models/blackoil/blackoilmodel.hh>
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 
 #include <cstddef>
@@ -54,7 +54,7 @@ class GenericTemperatureModel
 public:
     // the jacobian matrix
     using MatrixBlockTemp = MatrixBlock<Scalar, 1, 1>;
-    using EnergyMatrix = Dune::BCRSMatrix<MatrixBlockTemp>;
+    using EnergyMatrix = BlockSparseMatrix<MatrixBlockTemp>;
     using EnergyVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     static constexpr int dimWorld = Grid::dimensionworld;

opm/simulators/flow/GenericTracerModel.hpp

@@ -28,14 +28,13 @@
 #ifndef OPM_GENERIC_TRACER_MODEL_HPP
 #define OPM_GENERIC_TRACER_MODEL_HPP
 
-#include <dune/istl/bcrsmatrix.hh>
-
 #include <opm/grid/common/CartesianIndexMapper.hpp>
 
 #include <opm/input/eclipse/EclipseState/Phase.hpp>
 
 #include <opm/models/blackoil/blackoilmodel.hh>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/wells/WellTracerRate.hpp>
 
@@ -56,7 +55,7 @@ template<class Grid, class GridView, class DofMapper, class Stencil, class Fluid
 class GenericTracerModel {
 public:
     using TracerVectorSingle = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
-    using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 2, 2>>;
+    using TracerMatrix = BlockSparseMatrix<Opm::MatrixBlock<Scalar, 2, 2>>;
     using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar, 2>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     static constexpr int dimWorld = Grid::dimensionworld;

opm/simulators/flow/NonlinearSolver.hpp

@@ -22,7 +22,6 @@
 #define OPM_NON_LINEAR_SOLVER_HPP
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/Exceptions.hpp>

opm/simulators/linalg/AmgxPreconditioner.hpp

@@ -28,7 +28,6 @@
 #include <opm/simulators/linalg/gpuistl/AmgxInterface.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <amgx_c.h>
 

opm/simulators/linalg/BlockSparseMatrix.hpp

@@ -0,0 +1,50 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ */
+#ifndef OPM_BLOCKSPARSEMATRIX_HEADER_INCLUDED
+#define OPM_BLOCKSPARSEMATRIX_HEADER_INCLUDED
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+#include <dune/istl/ibcrsmatrix.hh>
+#else
+#include <dune/istl/bcrsmatrix.hh>
+#endif
+
+namespace Opm {
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+// Use the exprimental IBCRSMatrix, which has a better performance under distributed workloads.
+template<class B, class A = std::allocator<B>>
+using BlockSparseMatrix = Dune::IBCRSMatrix<B, uint_least32_t, A>;
+#else
+// Use the classic BCRSMatrix.
+template<class B, class A = std::allocator<B>>
+using BlockSparseMatrix = Dune::BCRSMatrix<B, A>;
+#endif
+
+} // namespace Opm
+
+#endif

opm/simulators/linalg/DILU.hpp

@@ -24,7 +24,6 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/unused.hh>
 #include <dune/common/version.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <opm/simulators/linalg/GraphColoring.hpp>
 

opm/simulators/linalg/FlexibleSolver_impl.hpp

@@ -27,6 +27,7 @@
 #include <opm/simulators/linalg/ilufirstelement.hh>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/Preconditioner2InverseOperator.hpp>
 #include <opm/simulators/linalg/WellOperators.hpp>
@@ -355,7 +356,7 @@ namespace Dune
 template<class Scalar, int N>
 using BV = Dune::BlockVector<Dune::FieldVector<Scalar, N>>;
 template<class Scalar, int N>
-using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, N, N>>;
+using OBM = Opm::BlockSparseMatrix<Opm::MatrixBlock<Scalar, N, N>>;
 
 // Sequential operators.
 template<class Scalar, int N>

opm/simulators/linalg/HyprePreconditioner.hpp

@@ -34,7 +34,6 @@
 #endif
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 
 #include <HYPRE.h>
 #include <HYPRE_krylov.h>

opm/simulators/linalg/ISTLSolver.cpp

@@ -29,6 +29,7 @@
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 #include <fmt/format.h>
@@ -164,7 +165,7 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
 }
 
 template<class Scalar, int Dim>
-using BM = Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>;
+using BM = BlockSparseMatrix<MatrixBlock<Scalar,Dim,Dim>>;
 template<class Scalar, int Dim>
 using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
 

opm/simulators/linalg/ISTLSolverGpuBridge.cpp

@@ -29,6 +29,7 @@
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 #include <fmt/format.h>
@@ -237,7 +238,7 @@ copyMatToBlockJac(const Matrix& mat, Matrix& blockJac)
 }
 
 template<class Scalar, int Dim>
-using BM = Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>;
+using BM = BlockSparseMatrix<MatrixBlock<Scalar,Dim,Dim>>;
 template<class Scalar, int Dim>
 using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
 

opm/simulators/linalg/MILU.cpp

@@ -23,13 +23,13 @@
 
 #include <dune/common/version.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/ilu.hh>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 
 #include <array>
 #include <algorithm>
@@ -267,8 +267,8 @@ void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
     INSTANTIATE_ILUN(__VA_ARGS__)
 
 #define INSTANTIATE_DIM(T,Dim)                                         \
-    INSTANTIATE_FULL(T,Dune::BCRSMatrix<MatrixBlock<T,Dim,Dim>>)       \
-    INSTANTIATE_FULL(T,Dune::BCRSMatrix<Dune::FieldMatrix<T,Dim,Dim>>)
+    INSTANTIATE_FULL(T,BlockSparseMatrix<MatrixBlock<T,Dim,Dim>>)       \
+    INSTANTIATE_FULL(T,BlockSparseMatrix<Dune::FieldMatrix<T,Dim,Dim>>)
 
 #define INSTANTIATE_TYPE(T)                 \
     template T identityFunctor(const T&);   \

opm/simulators/linalg/MatrixMarketSpecializations.hpp

@@ -44,6 +44,18 @@ namespace MatrixMarketImpl
         }
     };
 
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+    template <typename T, typename I, int i, int j, typename A>
+    struct mm_header_printer<IBCRSMatrix<Opm::MatrixBlock<T,i,j>, I, A>>
+    {
+        static void print(std::ostream& os)
+        {
+            os << "%%MatrixMarket matrix coordinate ";
+            os << mm_numeric_type<T>::str() << " general" << std::endl;
+        }
+    };
+#endif
+
     template <typename T, int i, int j, typename A>
     struct mm_block_structure_header<BCRSMatrix<Opm::MatrixBlock<T,i,j>, A>>
     {
@@ -54,6 +66,19 @@ namespace MatrixMarketImpl
             os << i << " " << j << std::endl;
         }
     };
+
+#if defined(OPM_USE_EXPERIMENTAL_IBCRSMATRIX) && OPM_USE_EXPERIMENTAL_IBCRSMATRIX
+    template <typename T, typename I, int i, int j, typename A>
+    struct mm_block_structure_header<IBCRSMatrix<Opm::MatrixBlock<T,i,j>, I, A>>
+    {
+        using M = IBCRSMatrix<Opm::MatrixBlock<T,i,j>, I, A>;
+        static void print(std::ostream& os, const M&)
+        {
+            os << "% ISTL_STRUCT blocked ";
+            os << i << " " << j << std::endl;
+        }
+    };
+#endif
 } // namespace MatrixMarketImpl
 
 namespace MatrixMarketImpl

opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp

@@ -27,7 +27,6 @@
 #include <opm/common/ErrorMacros.hpp>
 
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/paamg/amg.hh>
 
 #include <fstream>

opm/simulators/linalg/ParallelOverlappingILU0.cpp

@@ -22,17 +22,19 @@
 
 #include <opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
+
 #include <dune/istl/owneroverlapcopy.hh>
 
 namespace Opm
 {
 
 #define INSTANTIATE_PAR(T, Dim, ...)                                                     \
-  template class ParallelOverlappingILU0<Dune::BCRSMatrix<MatrixBlock<T,Dim,Dim>>,       \
+  template class ParallelOverlappingILU0<BlockSparseMatrix<MatrixBlock<T,Dim,Dim>>, \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          __VA_ARGS__>;                                   \
-  template class ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::FieldMatrix<T,Dim,Dim>>, \
+  template class ParallelOverlappingILU0<BlockSparseMatrix<Dune::FieldMatrix<T,Dim,Dim>>, \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          Dune::BlockVector<Dune::FieldVector<T,Dim>>,    \
                                          __VA_ARGS__>;

opm/simulators/linalg/PreconditionerFactory_impl.hpp

@@ -25,6 +25,7 @@
 
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 
+#include <opm/simulators/linalg/BlockSparseMatrix.hpp>
 #include <opm/simulators/linalg/DILU.hpp>
 #include <opm/simulators/linalg/ExtraSmoothers.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
@@ -198,21 +199,21 @@ PreconditionerFactory<Operator, Comm>::addCreator(const std::string& type, ParCr
 using CommSeq = Dune::Amg::SequentialInformation;
 
 template<class Scalar, int Dim>
-using OpFSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
+using OpFSeq = Dune::MatrixAdapter<BlockSparseMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
 template<class Scalar, int Dim>
-using OpBSeq = Dune::MatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpBSeq = Dune::MatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
 
 template<class Scalar, int Dim, bool overlap>
-using OpW = WellModelMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpW = WellModelMatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                    Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>>;
 
 template<class Scalar, int Dim, bool overlap>
-using OpWG = WellModelGhostLastMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpWG = WellModelGhostLastMatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                              overlap>;
@@ -221,24 +222,24 @@ using OpWG = WellModelGhostLastMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar
 using CommPar = Dune::OwnerOverlapCopyCommunication<int, int>;
 
 template<class Scalar, int Dim>
-using OpFPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
+using OpFPar = Dune::OverlappingSchwarzOperator<BlockSparseMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 CommPar>;
 
 template<class Scalar, int Dim>
-using OpBPar = Dune::OverlappingSchwarzOperator<Dune::BCRSMatrix<MatrixBlock<Scalar, Dim, Dim>>,
+using OpBPar = Dune::OverlappingSchwarzOperator<BlockSparseMatrix<MatrixBlock<Scalar, Dim, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>,
                                                 CommPar>;
 template<class Scalar, int Dim>
-using OpGLFPar = Opm::GhostLastMatrixAdapter<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,Dim,Dim>>,
+using OpGLFPar = Opm::GhostLastMatrixAdapter<BlockSparseMatrix<Dune::FieldMatrix<Scalar,Dim,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              CommPar>;
 
 template<class Scalar, int Dim>
-using OpGLBPar = Opm::GhostLastMatrixAdapter<Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>,
+using OpGLBPar = Opm::GhostLastMatrixAdapter<BlockSparseMatrix<MatrixBlock<Scalar,Dim,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>,
                                              CommPar>;