TestHelpers.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
7
8 #include <cstddef>
9 #include <optional>
10 #include <pup.h>
11
14 #include "Options/Options.hpp"
15 #include "Parallel/PupStlCpp17.hpp"
16 #include "Utilities/Gsl.hpp"
17
18 namespace TestHelpers::LinearSolver {
19
20 struct ApplyMatrix {
21  DenseMatrix<double> matrix;
22  void operator()(const gsl::not_null<DenseVector<double>*> result,
23  const DenseVector<double>& operand) const noexcept {
24  *result = matrix * operand;
25  }
26 };
27
28 // Use the exact inverse of the matrix as preconditioner. This should solve
29 // problems in 1 iteration.
31  void solve(const gsl::not_null<DenseVector<double>*> solution,
32  const ApplyMatrix& linear_operator,
33  const DenseVector<double>& source) const noexcept {
34  if (not inv_matrix_.has_value()) {
35  inv_matrix_ = blaze::inv(linear_operator.matrix);
36  }
37  *solution = *inv_matrix_ * source;
38  }
39
40  void reset() noexcept { inv_matrix_.reset(); }
41
42  void pup(PUP::er& p) noexcept { p | inv_matrix_; } // NOLINT
43
44  // Make option-creatable for factory tests
45  using options = tmpl::list<>;
46  static constexpr Options::String help{"halp"};
47
48  private:
49  mutable std::optional<DenseMatrix<double>> inv_matrix_{};
50 };
51
52 // Use the inverse of the diagonal as preconditioner.
54  void solve(const gsl::not_null<DenseVector<double>*> solution,
55  const ApplyMatrix& linear_operator,
56  const DenseVector<double>& source) const noexcept {
57  if (not inv_diagonal_.has_value()) {
58  inv_diagonal_ = DenseVector<double>(source.size(), 1.);
59  for (size_t i = 0; i < source.size(); ++i) {
60  (*inv_diagonal_)[i] /= linear_operator.matrix(i, i);
61  }
62  }
63  *solution = source;
64  for (size_t i = 0; i < solution->size(); ++i) {
65  (*solution)[i] *= (*inv_diagonal_)[i];
66  }
67  }
68
69  void reset() noexcept { inv_diagonal_.reset(); }
70
71  void pup(PUP::er& p) noexcept { p | inv_diagonal_; } // NOLINT
72
73  private:
74  mutable std::optional<DenseVector<double>> inv_diagonal_{};
75 };
76
77 // Run a few Richardson iterations as preconditioner.
79  RichardsonPreconditioner() = default;
80  RichardsonPreconditioner(const double relaxation_parameter,
81  const size_t num_iterations)
82  : relaxation_parameter_(relaxation_parameter),
83  num_iterations_(num_iterations) {}
84
85  void solve(
86  const gsl::not_null<DenseVector<double>*> initial_guess_in_solution_out,
87  const ApplyMatrix& linear_operator,
88  const DenseVector<double>& source) const noexcept {
89  for (size_t i = 0; i < num_iterations_; ++i) {
90  linear_operator(make_not_null(&correction_buffer_),
91  *initial_guess_in_solution_out);
92  correction_buffer_ *= -1.;
93  correction_buffer_ += source;
94  *initial_guess_in_solution_out +=
95  relaxation_parameter_ * correction_buffer_;
96  }
97  }
98
99  static void reset() noexcept {}
100
101  void pup(PUP::er& p) noexcept { // NOLINT
102  p | relaxation_parameter_;
103  p | num_iterations_;
104  }
105
106  private:
107  double relaxation_parameter_{};
108  size_t num_iterations_{};
109  mutable DenseVector<double> correction_buffer_{};
110 };
111
112 } // namespace TestHelpers::LinearSolver
TestHelpers::LinearSolver::RichardsonPreconditioner
Definition: TestHelpers.hpp:78
Options.hpp
TestHelpers::LinearSolver::JacobiPreconditioner
Definition: TestHelpers.hpp:53
TestingFramework.hpp
DenseVector.hpp
DenseMatrix.hpp
TestHelpers::LinearSolver::ExactInversePreconditioner
Definition: TestHelpers.hpp:30
TestHelpers::LinearSolver::ApplyMatrix
Definition: TestHelpers.hpp:20
DenseVector< double >
DenseMatrix< double >
cstddef
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
optional
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
gsl::not_null
Require a pointer to not be a nullptr