From b1ca77088901ba2acdeff634b88813f570b57540 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 10 Oct 2023 23:27:27 -0700 Subject: [PATCH 1/2] Write blaze::DynamicMatrix to CSV file --- src/DataStructures/DynamicMatrix.hpp | 17 +++++++++++++++++ .../Test_BlazeInteroperability.cpp | 6 ++++++ 2 files changed, 23 insertions(+) diff --git a/src/DataStructures/DynamicMatrix.hpp b/src/DataStructures/DynamicMatrix.hpp index aefe31a90607..cde8ca752b6b 100644 --- a/src/DataStructures/DynamicMatrix.hpp +++ b/src/DataStructures/DynamicMatrix.hpp @@ -75,3 +75,20 @@ struct Options::create_from_yaml> { return result; } }; + +/// Write a `blaze::DynamicMatrix` to a CSV file. +template +std::ostream& write_csv( + std::ostream& os, const blaze::DynamicMatrix& 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; +} diff --git a/tests/Unit/DataStructures/Test_BlazeInteroperability.cpp b/tests/Unit/DataStructures/Test_BlazeInteroperability.cpp index 34eb08b17f1e..870010a426c1 100644 --- a/tests/Unit/DataStructures/Test_BlazeInteroperability.cpp +++ b/tests/Unit/DataStructures/Test_BlazeInteroperability.cpp @@ -100,6 +100,12 @@ SPECTRE_TEST_CASE("Unit.DataStructures.BlazeInteroperability", "[[0., 1., 2.], [3., 0., 4.]]") == blaze::DynamicMatrix{{0., 1., 2.}, {3., 0., 4.}}); + // Test `write_csv` + std::stringstream ss{}; + blaze::DynamicMatrix matrix{{0., 1., 2.}, + {3., 0., 4.}}; + write_csv(ss, matrix); + CHECK(ss.str() == "0,1,2\n3,0,4\n"); } { INFO("CompressedMatrix"); From 6cf76626c3d6af552f5400ec58fb13244ec50b53 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 10 Oct 2023 23:27:59 -0700 Subject: [PATCH 2/2] Add option to output subdomain matrices in elliptic solver --- .../LinearSolver/ExplicitInverse.hpp | 50 ++++++++++++++++++- tests/InputFiles/Elasticity/BentBeam.yaml | 4 +- .../Elasticity/HalfSpaceMirror.yaml | 4 +- .../Poisson/ProductOfSinusoids1D.yaml | 6 ++- .../Poisson/ProductOfSinusoids2D.yaml | 4 +- .../Poisson/ProductOfSinusoids3D.yaml | 4 +- .../Punctures/MultiplePunctures.yaml | 4 +- tests/InputFiles/Xcts/BinaryBlackHole.yaml | 4 +- tests/InputFiles/Xcts/HeadOnBns.yaml | 4 +- tests/InputFiles/Xcts/KerrSchild.yaml | 4 +- tests/InputFiles/Xcts/TovStar.yaml | 4 +- .../Test_MinusLaplacian.cpp | 4 +- .../LinearSolver/Test_ExplicitInverse.cpp | 6 ++- .../Schwarz/Test_SchwarzAlgorithm.yaml | 3 +- 14 files changed, 90 insertions(+), 15 deletions(-) diff --git a/src/NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp b/src/NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp index d671e8e9f821..6eb031ae5a55 100644 --- a/src/NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp +++ b/src/NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp @@ -5,19 +5,26 @@ #include #include +#include +#include #include #include +#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 { using Base = LinearSolver; public: - using options = tmpl::list<>; + struct WriteMatrixToFile { + using type = Options::Auto; + 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; 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 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 { // 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::max()) { @@ -143,6 +164,7 @@ class ExplicitInverse : public LinearSolver { } private: + std::optional matrix_filename_{}; // Caches for successive solves of the same operator // NOLINTNEXTLINE(spectre-mutable) mutable size_t size_ = std::numeric_limits::max(); @@ -177,6 +199,30 @@ Convergence::HasConverged ExplicitInverse::solve( auto result_buffer = make_with_value(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 { + using DataBoxType = + std::decay_t>>; + if constexpr (tt::is_a_v) { + if constexpr (db::tag_is_retrievable_v) { + const auto& box = std::get<0>(operator_args); + return "_" + get_output(db::get(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_); diff --git a/tests/InputFiles/Elasticity/BentBeam.yaml b/tests/InputFiles/Elasticity/BentBeam.yaml index e4e3e6b0163b..15eb18367d6c 100644 --- a/tests/InputFiles/Elasticity/BentBeam.yaml +++ b/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: diff --git a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml index c6be537b001c..075662f9e791 100644 --- a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml +++ b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml @@ -101,7 +101,9 @@ LinearSolver: Restart: None Preconditioner: MinusLaplacian: - Solver: ExplicitInverse + Solver: + ExplicitInverse: + WriteMatrixToFile: None BoundaryConditions: Auto ObservePerCoreReductions: False diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml index 7b242e21a790..9cd4acfbaae8 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml +++ b/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 diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml index f7f3de788c39..36c06bb76b6b 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml +++ b/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 diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml index c615718f0e16..f9dc2ca7c269 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml +++ b/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 diff --git a/tests/InputFiles/Punctures/MultiplePunctures.yaml b/tests/InputFiles/Punctures/MultiplePunctures.yaml index 7c28c34d9542..9d6da7a81390 100644 --- a/tests/InputFiles/Punctures/MultiplePunctures.yaml +++ b/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 diff --git a/tests/InputFiles/Xcts/BinaryBlackHole.yaml b/tests/InputFiles/Xcts/BinaryBlackHole.yaml index f8eb33bedb16..2c6ee5946cb0 100644 --- a/tests/InputFiles/Xcts/BinaryBlackHole.yaml +++ b/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 diff --git a/tests/InputFiles/Xcts/HeadOnBns.yaml b/tests/InputFiles/Xcts/HeadOnBns.yaml index d882065f2e8b..b5855c217586 100644 --- a/tests/InputFiles/Xcts/HeadOnBns.yaml +++ b/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 diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml index 27bc77b77a9d..229230f9fa4e 100644 --- a/tests/InputFiles/Xcts/KerrSchild.yaml +++ b/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 diff --git a/tests/InputFiles/Xcts/TovStar.yaml b/tests/InputFiles/Xcts/TovStar.yaml index 8fd72f4ab58f..53142b22e528 100644 --- a/tests/InputFiles/Xcts/TovStar.yaml +++ b/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 diff --git a/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp b/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp index 74add3178474..490e61427989 100644 --- a/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp +++ b/tests/Unit/Elliptic/SubdomainPreconditioners/Test_MinusLaplacian.cpp @@ -422,7 +422,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.SubdomainPreconditioners.MinusLaplacian", const auto created = TestHelpers::test_creation>( "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(); diff --git a/tests/Unit/NumericalAlgorithms/LinearSolver/Test_ExplicitInverse.cpp b/tests/Unit/NumericalAlgorithms/LinearSolver/Test_ExplicitInverse.cpp index 5098a2ca6666..f85a18ec19eb 100644 --- a/tests/Unit/NumericalAlgorithms/LinearSolver/Test_ExplicitInverse.cpp +++ b/tests/Unit/NumericalAlgorithms/LinearSolver/Test_ExplicitInverse.cpp @@ -42,12 +42,16 @@ SPECTRE_TEST_CASE("Unit.LinearSolver.Serial.ExplicitInverse", const blaze::DynamicVector source{1., 2.}; const blaze::DynamicVector expected_solution{-1., 5.}; blaze::DynamicVector 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(matrix_file)), + std::istreambuf_iterator()); + CHECK(matrix_csv == "4 1\n3 1\n"); { INFO("Resetting"); ExplicitInverse<> resetting_solver{}; diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml b/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml index e0d5361523b4..3bc7415f8286 100644 --- a/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml +++ b/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