From 329036612b60333b34333a519dc7267210c761b2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 8 May 2020 10:27:11 +0200
Subject: [PATCH] Enable matrix dumping also without MPI.

Also remove the need for the reinterpret_cast by adding some extra
implementation details for more general block matrices.
---
 .../linalg/ISTLSolverEbosFlexible.hpp         | 13 +----
 .../linalg/WriteSystemMatrixHelper.hpp        | 56 +++++++++++++++++--
 2 files changed, 54 insertions(+), 15 deletions(-)

diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
index d20dd450c..9a9044fd5 100644
--- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
+++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
@@ -67,6 +67,8 @@ class ISTLSolverEbosFlexible
     using MatrixType = typename SparseMatrixAdapter::IstlMatrix;
 #if HAVE_MPI
     using Communication = Dune::OwnerOverlapCopyCommunication<int, int>;
+#else
+    using Communication = int; // Dummy type.
 #endif
     using SolverType = Dune::FlexibleSolver<MatrixType, VectorType>;
 
@@ -281,21 +283,14 @@ protected:
 
     void writeMatrix()
     {
-#ifdef HAVE_MPI
         const int verbosity = prm_.get<int>("verbosity");
         const bool write_matrix = verbosity > 10;
         if (write_matrix) {
-            using block_type = typename MatrixType::block_type;
-            using value_type = typename block_type::value_type;
-            using BaseBlockType = Dune::FieldMatrix<value_type,block_type::rows,block_type::cols>;
-            using BaseMatrixType = Dune::BCRSMatrix<BaseBlockType>;
-            const BaseMatrixType* matrix =  reinterpret_cast<BaseMatrixType*>(this->matrix_);
             Opm::Helper::writeSystem(this->simulator_, //simulator is only used to get names
-                                     *matrix,
+                                     *(this->matrix_),
                                      this->rhs_,
                                      comm_.get());
         }
-#endif
     }
 
     const Simulator& simulator_;
@@ -306,9 +301,7 @@ protected:
     VectorType rhs_;
     Dune::InverseOperatorResult res_;
     std::any parallelInformation_;
-#if HAVE_MPI
     std::unique_ptr<Communication> comm_;
-#endif
     std::vector<int> overlapRows_;
     std::vector<int> interiorRows_;
 }; // end ISTLSolverEbosFlexible
diff --git a/opm/simulators/linalg/WriteSystemMatrixHelper.hpp b/opm/simulators/linalg/WriteSystemMatrixHelper.hpp
index f76985176..13c07489e 100644
--- a/opm/simulators/linalg/WriteSystemMatrixHelper.hpp
+++ b/opm/simulators/linalg/WriteSystemMatrixHelper.hpp
@@ -22,6 +22,46 @@
 
 #include <dune/istl/matrixmarket.hh>
 
+namespace Dune
+{
+namespace MatrixMarketImpl
+{
+
+    template <typename Block, typename A>
+    struct mm_header_printer<BCRSMatrix<Block, A>>
+    {
+        static void print(std::ostream& os)
+        {
+            using Value = typename Block::value_type;
+            os << "%%MatrixMarket matrix coordinate ";
+            os << mm_numeric_type<Value>::str() << " general" << std::endl;
+        }
+    };
+
+    template <typename Block, typename A>
+    struct mm_block_structure_header<BCRSMatrix<Block, A>>
+    {
+        using M = BCRSMatrix<Block, A>;
+        static void print(std::ostream& os, const M&)
+        {
+            os << "% ISTL_STRUCT blocked ";
+            os << Block::rows << " " << Block::cols << std::endl;
+        }
+    };
+
+} // namespace MatrixMarketImpl
+
+template <typename Block, typename A>
+struct mm_multipliers<BCRSMatrix<Block, A>>
+{
+    enum {
+        rows = Block::rows,
+        cols = Block::cols
+    };
+};
+
+} // namespace Dune
+
 namespace Opm
 {
 namespace Helper
@@ -30,7 +70,7 @@ namespace Helper
     void writeSystem(const SimulatorType& simulator,
                      const MatrixType& matrix,
                      const VectorType& rhs,
-                     const Communicator* comm)
+                     [[maybe_unused]] const Communicator* comm)
     {
         std::string dir = simulator.problem().outputDir();
         if (dir == ".") {
@@ -56,18 +96,24 @@ namespace Helper
         std::string prefix = full_path.string();
         {
             std::string filename = prefix + "matrix_istl";
+#if HAVE_MPI
             if (comm != nullptr) { // comm is not set in serial runs
                 Dune::storeMatrixMarket(matrix, filename, *comm, true);
-            } else {
-                Dune::storeMatrixMarket(matrix, filename);
+            } else
+#endif
+            {
+                Dune::storeMatrixMarket(matrix, filename + ".mm");
             }
         }
         {
             std::string filename = prefix + "rhs_istl";
+#if HAVE_MPI
             if (comm != nullptr) { // comm is not set in serial runs
                 Dune::storeMatrixMarket(rhs, filename, *comm, true);
-            } else {
-                Dune::storeMatrixMarket(rhs, filename);
+            } else
+#endif
+            {
+                Dune::storeMatrixMarket(rhs, filename + ".mm");
             }
         }
     }