ResidualMonitorActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <tuple>
8 #include <utility>
9 
11 #include "DataStructures/DataBox/PrefixHelpers.hpp"
14 #include "Informer/Tags.hpp"
15 #include "Informer/Verbosity.hpp"
16 #include "NumericalAlgorithms/Convergence/Tags.hpp"
17 #include "Parallel/GlobalCache.hpp"
18 #include "Parallel/Invoke.hpp"
19 #include "Parallel/Printf.hpp"
20 #include "ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp"
21 #include "ParallelAlgorithms/LinearSolver/Observe.hpp"
23 #include "Utilities/EqualWithinRoundoff.hpp"
24 #include "Utilities/Gsl.hpp"
25 #include "Utilities/Requires.hpp"
26 
27 /// \cond
28 namespace tuples {
29 template <typename...>
30 class TaggedTuple;
31 } // namespace tuples
32 /// \endcond
33 
34 namespace LinearSolver::gmres::detail {
35 
36 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
37 struct InitializeResidualMagnitude {
38  private:
39  using fields_tag = FieldsTag;
40  using initial_residual_magnitude_tag =
43  using orthogonalization_history_tag =
45 
46  public:
47  template <typename ParallelComponent, typename DbTagsList,
48  typename Metavariables, typename ArrayIndex,
49  typename DataBox = db::DataBox<DbTagsList>,
50  Requires<db::tag_is_retrievable_v<initial_residual_magnitude_tag,
51  DataBox>> = nullptr>
52  static void apply(db::DataBox<DbTagsList>& box,
54  const ArrayIndex& /*array_index*/,
55  const double residual_magnitude) noexcept {
56  constexpr size_t iteration_id = 0;
57 
58  db::mutate<initial_residual_magnitude_tag>(
59  make_not_null(&box),
60  [residual_magnitude](
61  const gsl::not_null<double*> initial_residual_magnitude) noexcept {
62  *initial_residual_magnitude = residual_magnitude;
63  });
64 
65  LinearSolver::observe_detail::contribute_to_reduction_observer<
66  OptionsGroup, ParallelComponent>(iteration_id, residual_magnitude,
67  cache);
68 
69  // Determine whether the linear solver has already converged
70  Convergence::HasConverged has_converged{
71  get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
72  residual_magnitude, residual_magnitude};
73 
74  // Do some logging
76  ::Verbosity::Quiet)) {
77  Parallel::printf("%s initialized with residual: %e\n",
78  Options::name<OptionsGroup>(), residual_magnitude);
79  }
81  cache) >= ::Verbosity::Quiet)) {
82  Parallel::printf("%s has converged without any iterations: %s\n",
83  Options::name<OptionsGroup>(), has_converged);
84  }
85 
86  Parallel::receive_data<Tags::InitialOrthogonalization<OptionsGroup>>(
87  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
88  // NOLINTNEXTLINE(performance-move-const-arg)
89  std::make_tuple(residual_magnitude, std::move(has_converged)));
90  }
91 };
92 
93 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
94 struct StoreOrthogonalization {
95  private:
96  using fields_tag = FieldsTag;
97  using initial_residual_magnitude_tag =
100  using orthogonalization_history_tag =
102 
103  public:
104  template <typename ParallelComponent, typename DbTagsList,
105  typename Metavariables, typename ArrayIndex,
106  typename DataBox = db::DataBox<DbTagsList>,
107  Requires<db::tag_is_retrievable_v<orthogonalization_history_tag,
108  DataBox>> = nullptr>
109  static void apply(db::DataBox<DbTagsList>& box,
111  const ArrayIndex& /*array_index*/,
112  const size_t iteration_id,
113  const size_t orthogonalization_iteration_id,
114  const double orthogonalization) noexcept {
115  if (UNLIKELY(orthogonalization_iteration_id == 0)) {
116  // Append a row and a column to the orthogonalization history. Zero the
117  // entries that won't be set during the orthogonalization procedure below.
118  db::mutate<orthogonalization_history_tag>(
119  make_not_null(&box),
120  [iteration_id](const auto orthogonalization_history) noexcept {
121  orthogonalization_history->resize(iteration_id + 2,
122  iteration_id + 1);
123  for (size_t j = 0; j < orthogonalization_history->columns() - 1;
124  ++j) {
125  (*orthogonalization_history)(
126  orthogonalization_history->rows() - 1, j) = 0.;
127  }
128  });
129  }
130 
131  // While the orthogonalization procedure is not complete, store the
132  // orthogonalization, broadcast it back to all elements and return early
133  if (orthogonalization_iteration_id <= iteration_id) {
134  db::mutate<orthogonalization_history_tag>(
135  make_not_null(&box),
136  [orthogonalization, iteration_id, orthogonalization_iteration_id](
137  const auto orthogonalization_history) noexcept {
138  (*orthogonalization_history)(orthogonalization_iteration_id,
139  iteration_id) = orthogonalization;
140  });
141 
142  Parallel::receive_data<Tags::Orthogonalization<OptionsGroup>>(
143  Parallel::get_parallel_component<BroadcastTarget>(cache),
144  iteration_id, orthogonalization);
145  return;
146  }
147 
148  // At this point, the orthogonalization procedure is complete.
149  db::mutate<orthogonalization_history_tag>(
150  make_not_null(&box),
151  [orthogonalization, iteration_id, orthogonalization_iteration_id](
152  const auto orthogonalization_history) noexcept {
153  (*orthogonalization_history)(orthogonalization_iteration_id,
154  iteration_id) = sqrt(orthogonalization);
155  });
156 
157  // Perform a QR decomposition of the Hessenberg matrix that was built during
158  // the orthogonalization
159  const auto& orthogonalization_history =
160  get<orthogonalization_history_tag>(box);
161  const auto num_rows = orthogonalization_iteration_id + 1;
162  DenseMatrix<double> qr_Q;
163  DenseMatrix<double> qr_R;
164  blaze::qr(orthogonalization_history, qr_Q, qr_R);
165  // Compute the residual vector from the QR decomposition
166  DenseVector<double> beta(num_rows, 0.);
167  beta[0] = get<initial_residual_magnitude_tag>(box);
168  DenseVector<double> minres = blaze::inv(qr_R) * blaze::trans(qr_Q) * beta;
169  const double residual_magnitude =
170  blaze::length(beta - orthogonalization_history * minres);
171 
172  // At this point, the iteration is complete. We proceed with observing,
173  // logging and checking convergence before broadcasting back to the
174  // elements.
175 
176  const size_t completed_iterations = iteration_id + 1;
177  LinearSolver::observe_detail::contribute_to_reduction_observer<
178  OptionsGroup, ParallelComponent>(completed_iterations,
179  residual_magnitude, cache);
180 
181  // Determine whether the linear solver has converged
182  Convergence::HasConverged has_converged{
183  get<Convergence::Tags::Criteria<OptionsGroup>>(box),
184  completed_iterations, residual_magnitude,
185  get<initial_residual_magnitude_tag>(box)};
186 
187  // Do some logging
189  ::Verbosity::Quiet)) {
190  Parallel::printf("%s(%zu) iteration complete. Remaining residual: %e\n",
191  Options::name<OptionsGroup>(), completed_iterations,
192  residual_magnitude);
193  }
194  if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
195  cache) >= ::Verbosity::Quiet)) {
196  Parallel::printf("%s has converged in %zu iterations: %s\n",
197  Options::name<OptionsGroup>(), completed_iterations,
198  has_converged);
199  }
200 
201  Parallel::receive_data<Tags::FinalOrthogonalization<OptionsGroup>>(
202  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
203  std::make_tuple(sqrt(orthogonalization), std::move(minres),
204  // NOLINTNEXTLINE(performance-move-const-arg)
205  std::move(has_converged)));
206  }
207 };
208 
209 } // namespace LinearSolver::gmres::detail
std::apply
T apply(T... args)
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
sqrt
auto sqrt(const TensorExpression< T, typename T::type, tmpl::list<>, tmpl::list<>, tmpl::list<>> &t)
Returns the tensor expression representing the square root of a tensor expression that evaluates to a...
Definition: SquareRoot.hpp:80
Tags::Initial
Prefix indicating the initial value of a quantity.
Definition: Prefixes.hpp:85
GlobalCache.hpp
Parallel::printf
void printf(const std::string &format, Args &&... args)
Print an atomic message to stdout with C printf usage.
Definition: Printf.hpp:103
DenseVector.hpp
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
DenseMatrix.hpp
tuple
DenseVector< double >
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:68
Tags.hpp
DenseMatrix< double >
Printf.hpp
DataBox.hpp
cstddef
logging::Tags::Verbosity
Tag for putting Verbosity in a DataBox.
Definition: Tags.hpp:33
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:380
Gsl.hpp
LinearSolver::Tags::Magnitude
The magnitude w.r.t. the LinearSolver::inner_product
Definition: Tags.hpp:114
Requires.hpp
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
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
std::beta
T beta(T... args)
LinearSolver::Tags::OrthogonalizationHistory
A Hessenberg matrix built up during an orthogonalization procedure.
Definition: Tags.hpp:140