From 694efffd2718d8edcb81d91c995c58fec7a13809 Mon Sep 17 00:00:00 2001
From: IvBu <ivan.buntic@iws.uni-stuttgart.de>
Date: Thu, 27 Feb 2025 15:42:08 +0100
Subject: [PATCH] [doc] Add wiki page on printing the system matrix.
---
doc/doxygen/DoxygenDumuxLayout.xml | 1 +
doc/doxygen/pages/printing-system-matrix.md | 56 +++++++++++++++++++++
2 files changed, 57 insertions(+)
create mode 100644 doc/doxygen/pages/printing-system-matrix.md
diff --git a/doc/doxygen/DoxygenDumuxLayout.xml b/doc/doxygen/DoxygenDumuxLayout.xml
index 506067984e..2243a0458d 100644
--- a/doc/doxygen/DoxygenDumuxLayout.xml
+++ b/doc/doxygen/DoxygenDumuxLayout.xml
@@ -36,6 +36,7 @@ SPDX-License-Identifier: GPL-3.0-or-later
<tab type="user" url="@ref changing-property-name" visible="yes" title="Changing property name" intro=""/>
<tab type="user" url="@ref custom-input-data" visible="yes" title="Custom input data" intro=""/>
<tab type="user" url="@ref gmsh-with-alugrid" visible="yes" title="Gmsh with ALUGrid" intro=""/>
+ <tab type="user" url="@ref printing-system-matrix" visible="yes" title="Printing system matrix" intro=""/>
</tab>
<tab type="user" url="@ref citelist" visible="yes" title="Bibliography" intro=""/>
<tab type="namespaces" visible="yes" title="">
diff --git a/doc/doxygen/pages/printing-system-matrix.md b/doc/doxygen/pages/printing-system-matrix.md
new file mode 100644
index 0000000000..3c98c94539
--- /dev/null
+++ b/doc/doxygen/pages/printing-system-matrix.md
@@ -0,0 +1,56 @@
+# Printing system matrix
+
+Debugging may require inspecting the system (Jacobian) matrix from time to time.
+
+The following example shows a possible way to achieve this for a staggered-grid free-flow problem:
+
+We add the output within the `solveLinearSystem_` overload for `Dune::MultiTypeBlockMatrix` Jacobians in `newtonsolver.hh`:
+
+```c++
+ // MultiTypeBlockMatrix is not supported for printing, needs to be converted first
+
+ // create the bcrs matrix the IterativeSolver backend can handle
+ const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
+
+ // get the new matrix sizes
+ const std::size_t numRows = M.N();
+ assert(numRows == M.M());
+
+ // create the vector the IterativeSolver backend can handle
+ const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
+ assert(bTmp.size() == numRows);
+
+ // create a blockvector to which the linear solver writes the solution
+ using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
+ using BlockVector = typename Dune::BlockVector<VectorBlock>;
+ BlockVector y(numRows);
+
+ /////// BEGIN PRINTING /////////////////////////////////////////////////////
+ static const bool printmatrix = getParam<bool>("Problem.PrintMatrix", false);
+
+ if (printmatrix)
+ {
+ static int counter = 0;
+ const auto rank = Dune::MPIHelper::getCommunication().rank();
+
+ Dune::storeMatrixMarket(M, "matrix_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log");
+ Dune::storeMatrixMarket(bTmp, "rhs_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log");
+
+ ++counter;
+ }
+ /////// END PRINTING /////////////////////////////////////////////////////
+
+// ....
+```
+
+This will write out files for the matrix and the residual vector for each Newton iteration and each process. `A` and `b` can also be accessed via `assembler.jacobian()` and `assembler.residual()` respectively.
+
+You may then compare files using, e.g., `kompare` or `meld`.
+
+Keep in mind that the `MatrixMarket`format uses indices starting with `1` instead of `0`.
+
+For rather small matrices, you may also directly print the output to the terminal, giving you a graphical impression of the structure:
+
+```c++
+Dune::printmatrix(std::cout, M, "", "", 10/*width*/, 2/*precision*/);
+```
\ No newline at end of file
--
GitLab