Merge pull request #854 from nilsleiffischer/linsolv_namespace

Add linear solver namespace and IterationId

cmake/SpectreParseTests.py

@@ -25,6 +25,7 @@
                 "IO",
                 "LinearAlgebra",
                 "LinearOperators",
+                "LinearSolver",
                 "NumericalAlgorithms",
                 "Options",
                 "Parallel",

docs/GroupDefs.hpp

@@ -454,6 +454,37 @@
  * \brief Tags used for options parsed from the input file
  */
 
+/*!
+ * \defgroup LinearSolverGroup  Linear Solver
+ * \brief Algorithms to solve linear systems of equations
+ *
+ * \details In a way, the linear solver is for elliptic systems what time
+ * stepping is for the evolution code. This is because the DG scheme for an
+ * elliptic system reduces to a linear system of equations of the type
+ * \f$Ax=b\f$, where \f$A\f$ is a global matrix representing the DG
+ * discretization of the problem. Since this is one equation for each node in
+ * the computational domain it becomes unfeasible to numerically invert the
+ * global matrix \f$A\f$. Instead, we solve the problem iteratively so that we
+ * never need to construct \f$A\f$ globally but only need \f$Ax\f$ that can be
+ * evaluated locally by virtue of the DG formulation. This action of the
+ * operator we have to supply in each step of the iterative algorithms
+ * implemented here. It is where most of the computational cost goes and usually
+ * involves computing a volume contribution for each element and communicating
+ * fluxes with neighboring elements.
+ *
+ * In the iterative algorithms we usually don't work with the physical field
+ * \f$x\f$ directly. Instead we need to apply the operator to an internal
+ * variable defined by the respective algorithm. This variable is exposed as
+ * `LinearSolver::Tags::Operand`, and the algorithm expects that the computed
+ * operator action is written into
+ * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
+ * LinearSolver::Tags::Operand>` in each step.
+ *
+ * Since the iterative algorithms typically scale badly with increasing grid
+ * size, a preconditioner \f$P\f$ is needed such that \f$P^{-1}A\f$ is easier to
+ * invert.
+ */
+
 /// \defgroup MathFunctionsGroup Math Functions
 /// \brief Useful analytic functions
 

src/NumericalAlgorithms/CMakeLists.txt

@@ -5,5 +5,6 @@ add_subdirectory(Interpolation)
 add_subdirectory(DiscontinuousGalerkin)
 add_subdirectory(LinearAlgebra)
 add_subdirectory(LinearOperators)
+add_subdirectory(LinearSolver)
 add_subdirectory(RootFinding)
 add_subdirectory(Spectral)

src/NumericalAlgorithms/LinearSolver/CMakeLists.txt

@@ -0,0 +1,17 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(LIBRARY LinearSolver)
+
+set(LIBRARY_SOURCES
+    IterationId.cpp
+    )
+
+add_library(${LIBRARY} ${LIBRARY_SOURCES})
+
+target_link_libraries(
+  ${LIBRARY}
+  INTERFACE DataStructures
+  INTERFACE ErrorHandling
+  INTERFACE Utilities
+  )

src/NumericalAlgorithms/LinearSolver/IterationId.cpp

@@ -0,0 +1,53 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "NumericalAlgorithms/LinearSolver/IterationId.hpp"
+
+#include <boost/functional/hash.hpp>  // IWYU pragma: keep
+#include <cstddef>
+#include <ostream>  // IWYU pragma: keep
+#include <pup.h>  // IWYU pragma: keep
+
+namespace LinearSolver {
+
+void IterationId::pup(PUP::er& p) noexcept { p | step_number; }
+
+bool operator==(const IterationId& a, const IterationId& b) noexcept {
+  return a.step_number == b.step_number;
+}
+bool operator!=(const IterationId& a, const IterationId& b) noexcept {
+  return not(a == b);
+}
+
+bool operator<(const IterationId& a, const IterationId& b) noexcept {
+  return a.step_number < b.step_number;
+}
+bool operator<=(const IterationId& a, const IterationId& b) noexcept {
+  return not(b < a);
+}
+bool operator>(const IterationId& a, const IterationId& b) noexcept {
+  return b < a;
+}
+bool operator>=(const IterationId& a, const IterationId& b) noexcept {
+  return not(a < b);
+}
+
+std::ostream& operator<<(std::ostream& s, const IterationId& id) noexcept {
+  return s << id.step_number;
+}
+
+size_t hash_value(const IterationId& id) noexcept {
+  size_t h = 0;
+  boost::hash_combine(h, id.step_number);
+  return h;
+}
+
+}  // namespace LinearSolver
+
+// clang-tidy: do not modify std namespace (okay for hash)
+namespace std {  // NOLINT
+size_t hash<LinearSolver::IterationId>::operator()(
+    const LinearSolver::IterationId& id) const noexcept {
+  return boost::hash<LinearSolver::IterationId>{}(id);
+}
+}  // namespace std

src/NumericalAlgorithms/LinearSolver/IterationId.hpp

@@ -0,0 +1,52 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+/// \file
+/// Defines class IterationId.
+
+#pragma once
+
+#include <cstddef>
+#include <functional>
+#include <iosfwd>
+
+/// \cond
+namespace PUP {
+class er;
+}  // namespace PUP
+/// \endcond
+
+namespace LinearSolver {
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief Identifies a step in the linear solver algorithm
+ */
+struct IterationId {
+  size_t step_number{0};
+
+  IterationId() = default;
+
+  // clang-tidy: google-runtime-references
+  void pup(PUP::er& p) noexcept;  // NOLINT
+};
+
+bool operator==(const IterationId& a, const IterationId& b) noexcept;
+bool operator!=(const IterationId& a, const IterationId& b) noexcept;
+bool operator<(const IterationId& a, const IterationId& b) noexcept;
+bool operator<=(const IterationId& a, const IterationId& b) noexcept;
+bool operator>(const IterationId& a, const IterationId& b) noexcept;
+bool operator>=(const IterationId& a, const IterationId& b) noexcept;
+
+std::ostream& operator<<(std::ostream& s, const IterationId& id) noexcept;
+
+size_t hash_value(const IterationId& id) noexcept;
+
+}  // namespace LinearSolver
+
+namespace std {
+template <>
+struct hash<LinearSolver::IterationId> {
+  size_t operator()(const LinearSolver::IterationId& id) const noexcept;
+};
+}  // namespace std

src/NumericalAlgorithms/LinearSolver/Tags.hpp

@@ -0,0 +1,89 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+/// \file
+/// Defines DataBox tags for the linear solver
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/DataBox/DataBoxTag.hpp"
+#include "NumericalAlgorithms/LinearSolver/IterationId.hpp"
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief Functionality for solving linear systems of equations
+ */
+namespace LinearSolver {
+namespace Tags {
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief Holds an `IterationId` that identifies a step in the linear solver
+ * algorithm
+ */
+struct IterationId : db::SimpleTag {
+  static std::string name() noexcept { return "IterationId"; }
+  using type = LinearSolver::IterationId;
+};
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief The operand that the local linear operator \f$A\f$ is applied to
+ *
+ * \details The result of the operation should be wrapped in
+ * `LinearSolver::Tags::OperatorAppliedTo`.
+ */
+template <typename Tag>
+struct Operand : db::PrefixTag, db::SimpleTag {
+  static std::string name() noexcept {
+    return "LinearOperand(" + Tag::name() + ")";
+  }
+  using type = typename Tag::type;
+  using tag = Tag;
+  static constexpr bool should_be_sliced_to_boundary = true;
+};
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief The linear operator \f$A\f$ applied to the data in `Tag`
+ */
+template <typename Tag>
+struct OperatorAppliedTo : db::PrefixTag, db::SimpleTag {
+  static std::string name() noexcept {
+    return "LinearOperatorAppliedTo(" + Tag::name() + ")";
+  }
+  using type = typename Tag::type;
+  using tag = Tag;
+};
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief The residual \f$r=b - Ax\f$
+ */
+template <typename Tag>
+struct Residual : db::PrefixTag, db::SimpleTag {
+  static std::string name() noexcept {
+    return "LinearResidual(" + Tag::name() + ")";
+  }
+  using type = typename Tag::type;
+  using tag = Tag;
+};
+
+/*!
+ * \ingroup LinearSolverGroup
+ * \brief The magnitude square of the residual \f$\langle r,r\rangle_A\f$ w.r.t.
+ * the `LinearSolver::inner_product`
+ */
+template <typename Tag>
+struct ResidualMagnitudeSquare : db::PrefixTag, db::SimpleTag {
+  static std::string name() noexcept {
+    return "LinearResidualMagnitudeSquare(" + Tag::name() + ")";
+  }
+  using type = double;
+  using tag = Tag;
+};
+
+}  // namespace Tags
+}  // namespace LinearSolver

tests/Unit/NumericalAlgorithms/CMakeLists.txt

@@ -5,5 +5,6 @@ add_subdirectory(DiscontinuousGalerkin)
 add_subdirectory(Interpolation)
 add_subdirectory(LinearAlgebra)
 add_subdirectory(LinearOperators)
+add_subdirectory(LinearSolver)
 add_subdirectory(RootFinding)
 add_subdirectory(Spectral)

tests/Unit/NumericalAlgorithms/LinearSolver/CMakeLists.txt

@@ -0,0 +1,16 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(LIBRARY "Test_LinearSolver")
+
+set(LIBRARY_SOURCES
+  Test_IterationId.cpp
+  Test_Tags.cpp
+  )
+
+add_test_library(
+  ${LIBRARY}
+  "NumericalAlgorithms/LinearSolver/"
+  "${LIBRARY_SOURCES}"
+  "LinearSolver"
+  )

tests/Unit/NumericalAlgorithms/LinearSolver/Test_IterationId.cpp

@@ -0,0 +1,33 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "tests/Unit/TestingFramework.hpp"
+
+#include <functional>
+#include <string>
+
+#include "NumericalAlgorithms/LinearSolver/IterationId.hpp"
+#include "NumericalAlgorithms/LinearSolver/Tags.hpp"
+#include "Utilities/GetOutput.hpp"
+#include "tests/Unit/TestHelpers.hpp"
+
+SPECTRE_TEST_CASE("Unit.Numerical.LinearSolver.IterationId",
+                  "[Unit][NumericalAlgorithms][LinearSolver]") {
+  using Hash = std::hash<LinearSolver::IterationId>;
+
+  const LinearSolver::IterationId id{2};
+
+  CHECK(id.step_number == 2);
+  CHECK(id == LinearSolver::IterationId{2});
+  CHECK(id != LinearSolver::IterationId{3});
+  CHECK(Hash{}(id) == Hash{}(LinearSolver::IterationId{2}));
+  CHECK(Hash{}(id) != Hash{}(LinearSolver::IterationId{3}));
+  check_cmp(id, LinearSolver::IterationId{3});
+  CHECK(get_output(id) == "2");
+  test_serialization(id);
+  test_copy_semantics(id);
+
+  SECTION("Tag") {
+    CHECK(LinearSolver::Tags::IterationId::name() == "IterationId");
+  }
+}

tests/Unit/NumericalAlgorithms/LinearSolver/Test_Tags.cpp

@@ -0,0 +1,26 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "tests/Unit/TestingFramework.hpp"
+
+#include <string>
+
+#include "DataStructures/DataBox/DataBoxTag.hpp"
+#include "NumericalAlgorithms/LinearSolver/Tags.hpp"
+
+namespace {
+struct Tag : db::SimpleTag {
+  static std::string name() noexcept { return "Tag"; }
+  using type = int;
+};
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.Numerical.LinearSolver.Tags",
+                  "[Unit][NumericalAlgorithms][LinearSolver]") {
+  CHECK(LinearSolver::Tags::Operand<Tag>::name() == "LinearOperand(Tag)");
+  CHECK(LinearSolver::Tags::OperatorAppliedTo<Tag>::name() ==
+        "LinearOperatorAppliedTo(Tag)");
+  CHECK(LinearSolver::Tags::Residual<Tag>::name() == "LinearResidual(Tag)");
+  CHECK(LinearSolver::Tags::ResidualMagnitudeSquare<Tag>::name() ==
+        "LinearResidualMagnitudeSquare(Tag)");
+}