Merge pull request #5557 from nilsvu/subdomain_matrix_output

Add option to output subdomain matrices in elliptic solver

src/DataStructures/DynamicMatrix.hpp

@@ -75,3 +75,20 @@ struct Options::create_from_yaml<blaze::DynamicMatrix<Type, SO, Alloc, Tag>> {
     return result;
   }
 };
+
+/// Write a `blaze::DynamicMatrix` to a CSV file.
+template <typename Type, bool SO, typename Alloc, typename Tag>
+std::ostream& write_csv(
+    std::ostream& os, const blaze::DynamicMatrix<Type, SO, Alloc, Tag>& matrix,
+    const std::string& delimiter = ",") {
+  for (size_t i = 0; i < matrix.rows(); ++i) {
+    for (size_t j = 0; j < matrix.columns(); ++j) {
+      os << matrix(i, j);
+      if (j + 1 != matrix.columns()) {
+        os << delimiter;
+      }
+    }
+    os << '\n';
+  }
+  return os;
+}

src/NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp

@@ -5,19 +5,26 @@
 
 #include <algorithm>
 #include <cstddef>
+#include <fstream>
+#include <string>
 #include <tuple>
 #include <vector>
 
+#include "DataStructures/DataBox/DataBox.hpp"
 #include "DataStructures/DynamicMatrix.hpp"
 #include "DataStructures/DynamicVector.hpp"
 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
 #include "NumericalAlgorithms/LinearSolver/BuildMatrix.hpp"
 #include "NumericalAlgorithms/LinearSolver/LinearSolver.hpp"
+#include "Options/Auto.hpp"
 #include "Options/String.hpp"
+#include "Parallel/Tags/ArrayIndex.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeWithValue.hpp"
+#include "Utilities/NoSuchType.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
+#include "Utilities/Serialization/PupStlCpp17.hpp"
 #include "Utilities/TMPL.hpp"
 
 namespace LinearSolver::Serial {
@@ -73,19 +80,32 @@ class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
   using Base = LinearSolver<LinearSolverRegistrars>;
 
  public:
-  using options = tmpl::list<>;
+  struct WriteMatrixToFile {
+    using type = Options::Auto<std::string, Options::AutoLabel::None>;
+    static constexpr Options::String help =
+        "Write the matrix representation of the linear operator to a "
+        "space-delimited CSV file with this name. A '.txt' extension will be "
+        "added. Also a suffix with the element ID will be added if this linear "
+        "solver runs on an array element, so one file per element will be "
+        "written.";
+  };
+
+  using options = tmpl::list<WriteMatrixToFile>;
   static constexpr Options::String help =
       "Build a matrix representation of the linear operator and invert it "
       "directly. This means that the first solve has a large initialization "
       "cost, but all subsequent solves converge immediately.";
 
-  ExplicitInverse() = default;
   ExplicitInverse(const ExplicitInverse& /*rhs*/) = default;
   ExplicitInverse& operator=(const ExplicitInverse& /*rhs*/) = default;
   ExplicitInverse(ExplicitInverse&& /*rhs*/) = default;
   ExplicitInverse& operator=(ExplicitInverse&& /*rhs*/) = default;
   ~ExplicitInverse() = default;
 
+  explicit ExplicitInverse(
+      std::optional<std::string> matrix_filename = std::nullopt)
+      : matrix_filename_(std::move(matrix_filename)) {}
+
   /// \cond
   explicit ExplicitInverse(CkMigrateMessage* m) : Base(m) {}
   using PUP::able::register_constructor;
@@ -130,6 +150,7 @@ class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
 
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p) override {
+    p | matrix_filename_;
     p | size_;
     p | inverse_;
     if (p.isUnpacking() and size_ != std::numeric_limits<size_t>::max()) {
@@ -143,6 +164,7 @@ class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
   }
 
  private:
+  std::optional<std::string> matrix_filename_{};
   // Caches for successive solves of the same operator
   // NOLINTNEXTLINE(spectre-mutable)
   mutable size_t size_ = std::numeric_limits<size_t>::max();
@@ -177,6 +199,30 @@ Convergence::HasConverged ExplicitInverse<LinearSolverRegistrars>::solve(
     auto result_buffer = make_with_value<SourceType>(used_for_size, 0.);
     build_matrix(make_not_null(&inverse_), make_not_null(&operand_buffer),
                  make_not_null(&result_buffer), linear_operator, operator_args);
+    // Write to file before inverting
+    if (UNLIKELY(matrix_filename_.has_value())) {
+      const auto filename_suffix =
+          [&operator_args]() -> std::optional<std::string> {
+        using DataBoxType =
+            std::decay_t<tmpl::front<tmpl::list<OperatorArgs..., NoSuchType>>>;
+        if constexpr (tt::is_a_v<db::DataBox, DataBoxType>) {
+          if constexpr (db::tag_is_retrievable_v<Parallel::Tags::ArrayIndex,
+                                                 DataBoxType>) {
+            const auto& box = std::get<0>(operator_args);
+            return "_" + get_output(db::get<Parallel::Tags::ArrayIndex>(box));
+          } else {
+            (void)operator_args;
+            return std::nullopt;
+          }
+        } else {
+          (void)operator_args;
+          return std::nullopt;
+        }
+      }();
+      std::ofstream matrix_file(matrix_filename_.value() +
+                                filename_suffix.value_or("") + ".txt");
+      write_csv(matrix_file, inverse_, " ");
+    }
     // Directly invert the matrix
     try {
       blaze::invert(inverse_);

tests/InputFiles/Elasticity/BentBeam.yaml

@@ -75,7 +75,9 @@ LinearSolver:
     Iterations: 3
     MaxOverlap: 2
     Verbosity: Quiet
-    SubdomainSolver: ExplicitInverse
+    SubdomainSolver:
+      ExplicitInverse:
+        WriteMatrixToFile: None
     ObservePerCoreReductions: False
 
 EventsAndTriggers:

tests/InputFiles/Elasticity/HalfSpaceMirror.yaml

@@ -101,7 +101,9 @@ LinearSolver:
         Restart: None
         Preconditioner:
           MinusLaplacian:
-            Solver: ExplicitInverse
+            Solver:
+              ExplicitInverse:
+                WriteMatrixToFile: None
             BoundaryConditions: Auto
     ObservePerCoreReductions: False
 

tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml

@@ -8,6 +8,8 @@ Testing:
 ExpectedOutput:
   - PoissonProductOfSinusoids1DReductions.h5
   - PoissonProductOfSinusoids1DVolume0.h5
+  - SubdomainMatrix_[B0,(L1I0)].txt
+  - SubdomainMatrix_[B0,(L1I1)].txt
 OutputFileChecks:
   - Label: Discretization error
     Subfile: /ErrorNorms.dat
@@ -75,7 +77,9 @@ LinearSolver:
     Iterations: 3
     MaxOverlap: 2
     Verbosity: Quiet
-    SubdomainSolver: ExplicitInverse
+    SubdomainSolver:
+      ExplicitInverse:
+        WriteMatrixToFile: "SubdomainMatrix"
     ObservePerCoreReductions: False
 
 RadiallyCompressedCoordinates: None

tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml

@@ -69,7 +69,9 @@ LinearSolver:
     Iterations: 3
     MaxOverlap: 2
     Verbosity: Quiet
-    SubdomainSolver: ExplicitInverse
+    SubdomainSolver:
+      ExplicitInverse:
+        WriteMatrixToFile: None
     ObservePerCoreReductions: False
 
 RadiallyCompressedCoordinates: None

tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml

@@ -76,7 +76,9 @@ LinearSolver:
     Iterations: 3
     MaxOverlap: 2
     Verbosity: Quiet
-    SubdomainSolver: ExplicitInverse
+    SubdomainSolver:
+      ExplicitInverse:
+        WriteMatrixToFile: None
     ObservePerCoreReductions: False
 
 RadiallyCompressedCoordinates: None

tests/InputFiles/Punctures/MultiplePunctures.yaml

@@ -92,7 +92,9 @@ LinearSolver:
         Restart: None
         Preconditioner:
           MinusLaplacian:
-            Solver: ExplicitInverse
+            Solver:
+              ExplicitInverse:
+                WriteMatrixToFile: None
             BoundaryConditions: Auto
     SkipResets: True
     ObservePerCoreReductions: False

tests/InputFiles/Xcts/BinaryBlackHole.yaml

@@ -156,7 +156,9 @@ LinearSolver:
         Restart: None
         Preconditioner:
           MinusLaplacian:
-            Solver: ExplicitInverse
+            Solver:
+              ExplicitInverse:
+                WriteMatrixToFile: None
             BoundaryConditions: Auto
     SkipResets: True
     ObservePerCoreReductions: False

tests/InputFiles/Xcts/HeadOnBns.yaml

@@ -150,7 +150,9 @@ LinearSolver:
         Restart: None
         Preconditioner:
           MinusLaplacian:
-            Solver: ExplicitInverse
+            Solver:
+              ExplicitInverse:
+                WriteMatrixToFile: None
             BoundaryConditions: Auto
     SkipResets: True
     ObservePerCoreReductions: False

tests/InputFiles/Xcts/KerrSchild.yaml

@@ -107,7 +107,9 @@ LinearSolver:
         Restart: None
         Preconditioner:
           MinusLaplacian:
-            Solver: ExplicitInverse
+            Solver:
+              ExplicitInverse:
+                WriteMatrixToFile: None
             BoundaryConditions: Auto
     SkipResets: True
     ObservePerCoreReductions: False

tests/InputFiles/Xcts/TovStar.yaml

@@ -103,7 +103,9 @@ LinearSolver:
         Restart: None
         Preconditioner:
           MinusLaplacian:
-            Solver: ExplicitInverse
+            Solver:
+              ExplicitInverse:
+                WriteMatrixToFile: None
             BoundaryConditions: Auto
     SkipResets: True
     ObservePerCoreReductions: False

tests/Unit/DataStructures/Test_BlazeInteroperability.cpp

@@ -100,6 +100,12 @@ SPECTRE_TEST_CASE("Unit.DataStructures.BlazeInteroperability",
               "[[0., 1., 2.], [3., 0., 4.]]") ==
           blaze::DynamicMatrix<double, blaze::rowMajor>{{0., 1., 2.},
                                                         {3., 0., 4.}});
+    // Test `write_csv`
+    std::stringstream ss{};
+    blaze::DynamicMatrix<double, blaze::columnMajor> matrix{{0., 1., 2.},
+                                                            {3., 0., 4.}};
+    write_csv(ss, matrix);
+    CHECK(ss.str() == "0,1,2\n3,0,4\n");
   }
   {
     INFO("CompressedMatrix");

tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp

@@ -422,7 +422,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.SubdomainPreconditioners.MinusLaplacian",
     const auto created =
         TestHelpers::test_creation<std::unique_ptr<LinearSolverType>>(
             "MinusLaplacian:\n"
-            "  Solver: ExplicitInverse\n"
+            "  Solver:\n"
+            "    ExplicitInverse:\n"
+            "      WriteMatrixToFile: None\n"
             "  BoundaryConditions: Auto");
     const auto serialized = serialize_and_deserialize(created);
     const auto cloned = serialized->get_clone();

tests/Unit/NumericalAlgorithms/LinearSolver/Test_ExplicitInverse.cpp

@@ -42,12 +42,16 @@ SPECTRE_TEST_CASE("Unit.LinearSolver.Serial.ExplicitInverse",
     const blaze::DynamicVector<double> source{1., 2.};
     const blaze::DynamicVector<double> expected_solution{-1., 5.};
     blaze::DynamicVector<double> solution(2);
-    const ExplicitInverse<> solver{};
+    const ExplicitInverse<> solver{"Matrix"};
     const auto has_converged =
         solver.solve(make_not_null(&solution), linear_operator, source);
     REQUIRE(has_converged);
     CHECK_MATRIX_APPROX(solver.matrix_representation(), blaze::inv(matrix));
     CHECK_ITERABLE_APPROX(solution, expected_solution);
+    std::ifstream matrix_file("Matrix.txt");
+    std::string matrix_csv((std::istreambuf_iterator<char>(matrix_file)),
+                           std::istreambuf_iterator<char>());
+    CHECK(matrix_csv == "4 1\n3 1\n");
     {
       INFO("Resetting");
       ExplicitInverse<> resetting_solver{};

tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml

@@ -68,7 +68,8 @@ SchwarzSmoother:
       Preconditioner:
         # Preconditioning with the explicitly-built inverse matrix, so all
         # subdomain solves should converge immediately
-        ExplicitInverse
+        ExplicitInverse:
+          WriteMatrixToFile: None
   ObservePerCoreReductions: False
 
 ConvergenceReason: NumIterations