Merge pull request #2594 from atgeirr/writematrix-revised

Add a writeSystem() helper and use it to write system matrices.

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -24,6 +24,7 @@
 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/setupPropertyTree.hpp>
+#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
 
 #include <opm/common/ErrorMacros.hpp>
 
@@ -199,9 +200,11 @@ public:
         if (recreate_solver || !solver_) {
             if (isParallel()) {
 #if HAVE_MPI
+                matrix_ = &mat.istlMatrix();
                 solver_.reset(new SolverType(prm_, mat.istlMatrix(), weightsCalculator, *comm_));
 #endif
             } else {
+                matrix_ = &mat.istlMatrix();
                 solver_.reset(new SolverType(prm_, mat.istlMatrix(), weightsCalculator));
             }
             rhs_ = b;
@@ -214,6 +217,7 @@ public:
     bool solve(VectorType& x)
     {
         solver_->apply(x, rhs_, res_);
+        this->writeMatrix();
         return res_.converged;
     }
 
@@ -275,8 +279,27 @@ protected:
         return weights;
     }
 
-    const Simulator& simulator_;
+    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->rhs_,
+                                     comm_.get());
+        }
+#endif
+    }
 
+    const Simulator& simulator_;
+    MatrixType* matrix_;
     std::unique_ptr<SolverType> solver_;
     FlowLinearSolverParameters parameters_;
     boost::property_tree::ptree prm_;

opm/simulators/linalg/WriteSystemMatrixHelper.hpp

@@ -0,0 +1,78 @@
+/*
+  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
+
+  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_WRITESYSTEMMATRIXHELPER_HEADER_INCLUDED
+#define OPM_WRITESYSTEMMATRIXHELPER_HEADER_INCLUDED
+
+#include <dune/istl/matrixmarket.hh>
+
+namespace Opm
+{
+namespace Helper
+{
+    template <class SimulatorType, class MatrixType, class VectorType, class Communicator>
+    void writeSystem(const SimulatorType& simulator,
+                     const MatrixType& matrix,
+                     const VectorType& rhs,
+                     const Communicator* comm)
+    {
+        std::string dir = simulator.problem().outputDir();
+        if (dir == ".") {
+            dir = "";
+        } else if (!dir.empty() && dir.back() != '/') {
+            dir += "/";
+        }
+        namespace fs = Opm::filesystem;
+        fs::path output_dir(dir);
+        fs::path subdir("reports");
+        output_dir = output_dir / subdir;
+        if (!(fs::exists(output_dir))) {
+            fs::create_directory(output_dir);
+        }
+        // Combine and return.
+        std::ostringstream oss;
+        oss << "prob_" << simulator.episodeIndex() << "_time_";
+        oss << std::setprecision(15) << std::setw(12) << std::setfill('0') << simulator.time() << "_";
+        int nit = simulator.model().newtonMethod().numIterations();
+        oss << "_nit_" << nit << "_";
+        std::string output_file(oss.str());
+        fs::path full_path = output_dir / output_file;
+        std::string prefix = full_path.string();
+        {
+            std::string filename = prefix + "matrix_istl";
+            if (comm != nullptr) { // comm is not set in serial runs
+                Dune::storeMatrixMarket(matrix, filename, *comm, true);
+            } else {
+                Dune::storeMatrixMarket(matrix, filename);
+            }
+        }
+        {
+            std::string filename = prefix + "rhs_istl";
+            if (comm != nullptr) { // comm is not set in serial runs
+                Dune::storeMatrixMarket(rhs, filename, *comm, true);
+            } else {
+                Dune::storeMatrixMarket(rhs, filename);
+            }
+        }
+    }
+
+
+} // namespace Helper
+} // namespace Opm
+#endif