Add print dilu levelsets option

CMakeLists_files.cmake

@@ -506,6 +506,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/bda/WellContributions.hpp
   opm/simulators/linalg/amgcpr.hh
   opm/simulators/linalg/DILU.hpp
+  opm/simulators/linalg/DILUUtil.hpp
   opm/simulators/linalg/twolevelmethodcpr.hh
   opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp
   opm/simulators/linalg/ExtraSmoothers.hpp

opm/simulators/linalg/DILU.hpp

@@ -19,6 +19,7 @@
 
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
+#include <opm/simulators/linalg/DILUUtil.hpp>
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 
 #include <dune/common/fmatrix.hh>
@@ -36,10 +37,15 @@
 #include <omp.h>
 #endif
 
+static std::once_flag print_dilu_parallelism;
+
 // TODO: rewrite factory and constructor to keep track of a number of threads variable
 namespace Dune
 {
 
+// Global variable in this file to ensure that the parallelism strucutre of the matrix
+// is only printed once.
+
 /*! \brief The OpenMP thread parallelized DILU preconditioner.
  *  \details Safe to run serially without OpenMP. When run in parallel
              the matrix is assumed to be symmetric.
@@ -71,10 +77,10 @@ class MultithreadDILU : public PreconditionerWithUpdate<X, Y>
         OPM_TIMEBLOCK(prec_construct);
         // TODO: rewrite so this value is set by an argument to the constructor
 #if HAVE_OPENMP
-        use_multithreading = omp_get_max_threads() > 1;
+        use_multithreading_ = omp_get_max_threads() > 1;
 #endif
 
-        if (use_multithreading) {
+        if (use_multithreading_) {
             A_reordered_.emplace(A_.N(), A_.N(), A_.nonzeroes(), M::row_wise);
 
             //! Assuming symmetric matrices using a lower triangular coloring to construct
@@ -103,14 +109,20 @@ class MultithreadDILU : public PreconditionerWithUpdate<X, Y>
         update();
     }
 
+    MultithreadDILU(const M& A, int verbosity) : MultithreadDILU(A){
+        if (verbosity > 0 && use_multithreading_){
+            std::call_once(print_dilu_parallelism, Opm::DILUUtils::writeSparseTableRowSizesToFile<size_t>, &level_sets_);
+        }
+    }
+
     /*!
        \brief Update the preconditioner.
        \copydoc Preconditioner::update()
     */
     void update() override
     {
         OPM_TIMEBLOCK(prec_update);
-        if (use_multithreading) {
+        if (use_multithreading_) {
             parallelUpdate();
         } else {
             serialUpdate();
@@ -135,7 +147,7 @@ class MultithreadDILU : public PreconditionerWithUpdate<X, Y>
     void apply(X& v, const Y& d) override
     {
         OPM_TIMEBLOCK(prec_apply);
-        if (use_multithreading) {
+        if (use_multithreading_) {
             parallelApply(v, d);
         } else {
             serialApply(v, d);
@@ -177,7 +189,7 @@ class MultithreadDILU : public PreconditionerWithUpdate<X, Y>
     //! \brief converts from index in natural ordered structure to index reordered strucutre
     std::vector<std::size_t> natural_to_reorder_;
     //! \brief Boolean value describing whether or not to use multithreaded version of functions
-    bool use_multithreading{false};
+    bool use_multithreading_{false};
 
     void serialUpdate()
     {

opm/simulators/linalg/DILUUtil.hpp

@@ -0,0 +1,58 @@
+
+/*
+  Copyright 2022-2023 SINTEF 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_DILUUTIL_HEADER_INCLUDED
+#define OPM_DILUUTIL_HEADER_INCLUDED
+
+#include <config.h>
+#include <fstream>
+#include <mpi.h>
+#include <opm/grid/utility/SparseTable.hpp>
+
+namespace Opm::DILUUtils{
+
+// TODO: make proper doxygen
+/// @brief This function writes to a file a sparse table and the sizes of each row
+/// The intended use is to be called from a DILU preconditioner to explore how
+/// parallelizable the upper and lower solves are. For now the function creats a file
+/// in the directory from which flow was called, only intended to be done when verbosity>0
+/// @tparam T Any type, it does not matter what objects are stored in the table
+/// @param sparse_table A pointer to a SparseTable, this function does not change the object
+template <class T>
+void writeSparseTableRowSizesToFile(const Opm::SparseTable<T> *sparse_table){
+
+        int rank = 0;
+        int size = sparse_table->size();
+#ifdef HAVE_MPI
+        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+#endif
+        std::string filename = "DILU_parallelism_distribution_on_rank_" + std::to_string(rank) + ".txt";
+        std::ofstream file(filename);
+
+        // print brief exeplenation of how to interpret the data in the file
+        file << "First is the number of levels, then the size of each level set" << std::endl;
+        file << size << std::endl;
+
+        for (int i = 0; i < size; i++){
+            file << sparse_table->rowSize(i) << std::endl;
+        }
+}
+
+template void writeSparseTableRowSizesToFile(const Opm::SparseTable<size_t>*);
+
+} // namespace Opm::DILUUtils
+
+#endif

opm/simulators/linalg/PreconditionerFactory_impl.hpp

@@ -177,8 +177,8 @@ struct StandardPreconditioners
           return createParILU(op, prm, comm, prm.get<int>("ilulevel", 0));
         });
         F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
-          DUNE_UNUSED_PARAMETER(prm);
-          return wrapBlockPreconditioner<MultithreadDILU<M, V, V>>(comm, op.getmat());
+          const int preconditioner_verbosity = prm.get<int>("verbosity", 0);
+          return wrapBlockPreconditioner<MultithreadDILU<M, V, V>>(comm, op.getmat(), preconditioner_verbosity);
         });
         F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&,
                      std::size_t, const C& comm) {
@@ -286,9 +286,10 @@ struct StandardPreconditioners
         });
 
         F::addCreator("CUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
+            const int preconditioner_verbosity = prm.get<int>("verbosity", 0);
             using field_type = typename V::field_type;
             using CuDILU = typename Opm::cuistl::CuDILU<M, Opm::cuistl::CuVector<field_type>, Opm::cuistl::CuVector<field_type>>;
-            auto cuDILU = std::make_shared<CuDILU>(op.getmat());
+            auto cuDILU = std::make_shared<CuDILU>(op.getmat(), preconditioner_verbosity);
 
             auto adapted = std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CuDILU>>(cuDILU);
             auto wrapped = std::make_shared<Opm::cuistl::CuBlockPreconditioner<V, V, Comm>>(adapted, comm);
@@ -373,8 +374,8 @@ struct StandardPreconditioners<Operator,Dune::Amg::SequentialInformation>
                 op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
         });
         F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
-            DUNE_UNUSED_PARAMETER(prm);
-            return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat());
+            const int preconditioner_verbosity = prm.get<int>("verbosity", 0);
+            return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat(), preconditioner_verbosity);
         });
         F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
             const int n = prm.get<int>("repeats", 1);
@@ -516,9 +517,10 @@ struct StandardPreconditioners<Operator,Dune::Amg::SequentialInformation>
         });
 
         F::addCreator("CUDILU", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {
+            const int preconditioner_verbosity = prm.get<int>("verbosity", 0);
             using field_type = typename V::field_type;
             using CUDILU = typename Opm::cuistl::CuDILU<M, Opm::cuistl::CuVector<field_type>, Opm::cuistl::CuVector<field_type>>;
-            return std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CUDILU>>(std::make_shared<CUDILU>(op.getmat()));
+            return std::make_shared<Opm::cuistl::PreconditionerAdapter<V, V, CUDILU>>(std::make_shared<CUDILU>(op.getmat(), preconditioner_verbosity));
         });
 
         F::addCreator("CUDILUFloat", [](const O& op, [[maybe_unused]] const P& prm, const std::function<V()>&, std::size_t) {

opm/simulators/linalg/cuistl/CuDILU.cpp

@@ -28,8 +28,14 @@
 #include <opm/simulators/linalg/cuistl/detail/cusparse_matrix_operations.hpp>
 #include <opm/simulators/linalg/cuistl/detail/safe_conversion.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <string>
 #include <vector>
+#include <opm/simulators/linalg/DILUUtil.hpp>
+#include <mutex>
 
+// Global variable in this file to ensure that the parallelism strucutre of the matrix
+// is only printed once.
+std::once_flag printDiluParallelism;
 namespace
 {
 std::vector<int>
@@ -125,6 +131,13 @@ CuDILU<M, X, Y, l>::CuDILU(const M& A)
     update();
 }
 
+template <class M, class X, class Y, int l>
+CuDILU<M, X, Y, l>::CuDILU(const M& A, int verbosity) : CuDILU(A){
+    if (verbosity > 0){
+        std::call_once(printDiluParallelism, Opm::DILUUtils::writeSparseTableRowSizesToFile<size_t>, &m_levelSets);
+    }
+}
+
 template <class M, class X, class Y, int l>
 void
 CuDILU<M, X, Y, l>::pre([[maybe_unused]] X& x, [[maybe_unused]] Y& b)

opm/simulators/linalg/cuistl/CuDILU.hpp

@@ -60,10 +60,17 @@ class CuDILU : public Dune::PreconditionerWithUpdate<X, Y>
     //!
     //!  Constructor gets all parameters to operate the prec.
     //! \param A The matrix to operate on.
-    //! \param w The relaxation factor.
     //!
     explicit CuDILU(const M& A);
 
+    //! \brief Constructor accounting for verbosity.
+    //!
+    //!  Constructor gets all parameters to operate the prec.
+    //! \param A The matrix to operate on.
+    //! \param verbosity The verbosity level.
+    //!
+    explicit CuDILU(const M& A, int verbosity);
+
     //! \brief Prepare the preconditioner.
     //! \note Does nothing at the time being.
     void pre(X& x, Y& b) override;