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>(iteration_id, residual_magnitude, cache);
67 
68  // Determine whether the linear solver has already converged
69  Convergence::HasConverged has_converged{
70  get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
71  residual_magnitude, residual_magnitude};
72 
73  // Do some logging
75  ::Verbosity::Verbose)) {
76  Parallel::printf("Linear solver '" +
77  Options::name<OptionsGroup>() +
78  "' initialized with residual: %e\n",
79  residual_magnitude);
80  }
82  cache) >= ::Verbosity::Quiet)) {
83  Parallel::printf("The linear solver '" +
84  Options::name<OptionsGroup>() +
85  "' has converged without any iterations: %s\n",
86  has_converged);
87  }
88 
89  Parallel::receive_data<Tags::InitialOrthogonalization<OptionsGroup>>(
90  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
91  // NOLINTNEXTLINE(performance-move-const-arg)
92  std::make_tuple(residual_magnitude, std::move(has_converged)));
93  }
94 };
95 
96 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
97 struct StoreOrthogonalization {
98  private:
99  using fields_tag = FieldsTag;
100  using initial_residual_magnitude_tag =
103  using orthogonalization_history_tag =
105 
106  public:
107  template <typename ParallelComponent, typename DbTagsList,
108  typename Metavariables, typename ArrayIndex,
109  typename DataBox = db::DataBox<DbTagsList>,
110  Requires<db::tag_is_retrievable_v<orthogonalization_history_tag,
111  DataBox>> = nullptr>
112  static void apply(db::DataBox<DbTagsList>& box,
114  const ArrayIndex& /*array_index*/,
115  const size_t iteration_id,
116  const size_t orthogonalization_iteration_id,
117  const double orthogonalization) noexcept {
118  if (UNLIKELY(orthogonalization_iteration_id == 0)) {
119  // Append a row and a column to the orthogonalization history. Zero the
120  // entries that won't be set during the orthogonalization procedure below.
121  db::mutate<orthogonalization_history_tag>(
122  make_not_null(&box),
123  [iteration_id](const auto orthogonalization_history) noexcept {
124  orthogonalization_history->resize(iteration_id + 2,
125  iteration_id + 1);
126  for (size_t j = 0; j < orthogonalization_history->columns() - 1;
127  ++j) {
128  (*orthogonalization_history)(
129  orthogonalization_history->rows() - 1, j) = 0.;
130  }
131  });
132  }
133 
134  // While the orthogonalization procedure is not complete, store the
135  // orthogonalization, broadcast it back to all elements and return early
136  if (orthogonalization_iteration_id <= iteration_id) {
137  db::mutate<orthogonalization_history_tag>(
138  make_not_null(&box),
139  [orthogonalization, iteration_id, orthogonalization_iteration_id](
140  const auto orthogonalization_history) noexcept {
141  (*orthogonalization_history)(orthogonalization_iteration_id,
142  iteration_id) = orthogonalization;
143  });
144 
145  Parallel::receive_data<Tags::Orthogonalization<OptionsGroup>>(
146  Parallel::get_parallel_component<BroadcastTarget>(cache),
147  iteration_id, orthogonalization);
148  return;
149  }
150 
151  // At this point, the orthogonalization procedure is complete.
152  db::mutate<orthogonalization_history_tag>(
153  make_not_null(&box),
154  [orthogonalization, iteration_id, orthogonalization_iteration_id](
155  const auto orthogonalization_history) noexcept {
156  (*orthogonalization_history)(orthogonalization_iteration_id,
157  iteration_id) = sqrt(orthogonalization);
158  });
159 
160  // Perform a QR decomposition of the Hessenberg matrix that was built during
161  // the orthogonalization
162  const auto& orthogonalization_history =
163  get<orthogonalization_history_tag>(box);
164  const auto num_rows = orthogonalization_iteration_id + 1;
165  DenseMatrix<double> qr_Q;
166  DenseMatrix<double> qr_R;
167  blaze::qr(orthogonalization_history, qr_Q, qr_R);
168  // Compute the residual vector from the QR decomposition
169  DenseVector<double> beta(num_rows, 0.);
170  beta[0] = get<initial_residual_magnitude_tag>(box);
171  DenseVector<double> minres = blaze::inv(qr_R) * blaze::trans(qr_Q) * beta;
172  const double residual_magnitude =
173  blaze::length(beta - orthogonalization_history * minres);
174 
175  // At this point, the iteration is complete. We proceed with observing,
176  // logging and checking convergence before broadcasting back to the
177  // elements.
178 
179  const size_t completed_iterations = iteration_id + 1;
180  LinearSolver::observe_detail::contribute_to_reduction_observer<
181  OptionsGroup>(completed_iterations, residual_magnitude, cache);
182 
183  // Determine whether the linear solver has converged
184  Convergence::HasConverged has_converged{
185  get<Convergence::Tags::Criteria<OptionsGroup>>(box),
186  completed_iterations, residual_magnitude,
187  get<initial_residual_magnitude_tag>(box)};
188 
189  // Do some logging
191  ::Verbosity::Verbose)) {
192  Parallel::printf("Linear solver '" +
193  Options::name<OptionsGroup>() +
194  "' iteration %zu done. Remaining residual: %e\n",
195  completed_iterations, residual_magnitude);
196  }
197  if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
198  cache) >= ::Verbosity::Quiet)) {
199  Parallel::printf("The linear solver '" +
200  Options::name<OptionsGroup>() +
201  "' has converged in %zu iterations: %s\n",
202  completed_iterations, has_converged);
203  }
204 
205  Parallel::receive_data<Tags::FinalOrthogonalization<OptionsGroup>>(
206  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
207  std::make_tuple(sqrt(orthogonalization), std::move(minres),
208  // NOLINTNEXTLINE(performance-move-const-arg)
209  std::move(has_converged)));
210  }
211 };
212 
213 } // 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:638
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
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
std::sqrt
T sqrt(T... args)
DenseMatrix< double >
Printf.hpp
DataBox.hpp
cstddef
logging::Tags::Verbosity
Tag for putting Verbosity in a DataBox.
Definition: Tags.hpp:33
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