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"
12 #include "Informer/Tags.hpp"
13 #include "Informer/Verbosity.hpp"
14 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
15 #include "NumericalAlgorithms/Convergence/Tags.hpp"
16 #include "Options/Options.hpp"
17 #include "Parallel/GlobalCache.hpp"
18 #include "Parallel/Invoke.hpp"
19 #include "Parallel/Printf.hpp"
20 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/Tags/InboxTags.hpp"
21 #include "ParallelAlgorithms/LinearSolver/Observe.hpp"
23 #include "Utilities/EqualWithinRoundoff.hpp"
24 #include "Utilities/Functional.hpp"
25 #include "Utilities/Gsl.hpp"
26 #include "Utilities/Requires.hpp"
27 
28 /// \cond
29 namespace tuples {
30 template <typename...>
31 class TaggedTuple;
32 } // namespace tuples
33 /// \endcond
34 
35 namespace LinearSolver::cg::detail {
36 
37 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
38 struct InitializeResidual {
39  private:
40  using fields_tag = FieldsTag;
41  using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
43  using initial_residual_magnitude_tag =
46 
47  public:
48  template <typename ParallelComponent, typename DbTagsList,
49  typename Metavariables, typename ArrayIndex,
50  typename DataBox = db::DataBox<DbTagsList>,
52  nullptr>
53  static void apply(db::DataBox<DbTagsList>& box,
55  const ArrayIndex& /*array_index*/,
56  const double residual_square) noexcept {
57  constexpr size_t iteration_id = 0;
58  const double residual_magnitude = sqrt(residual_square);
59 
60  db::mutate<residual_square_tag, initial_residual_magnitude_tag>(
61  make_not_null(&box),
62  [residual_square, residual_magnitude](
63  const gsl::not_null<double*> local_residual_square,
64  const gsl::not_null<double*> initial_residual_magnitude) noexcept {
65  *local_residual_square = residual_square;
66  *initial_residual_magnitude = residual_magnitude;
67  });
68 
69  LinearSolver::observe_detail::contribute_to_reduction_observer<
70  OptionsGroup>(iteration_id, residual_magnitude, cache);
71 
72  // Determine whether the linear solver has converged
73  Convergence::HasConverged has_converged{
74  get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
75  residual_magnitude, residual_magnitude};
76 
77  // Do some logging
79  ::Verbosity::Verbose)) {
80  Parallel::printf("Linear solver '" +
81  Options::name<OptionsGroup>() +
82  "' initialized with residual: %e\n",
83  residual_magnitude);
84  }
86  cache) >= ::Verbosity::Quiet)) {
87  Parallel::printf("The linear solver '" +
88  Options::name<OptionsGroup>() +
89  "' has converged without any iterations: %s\n",
90  has_converged);
91  }
92 
93  Parallel::receive_data<Tags::InitialHasConverged<OptionsGroup>>(
94  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
95  // NOLINTNEXTLINE(performance-move-const-arg)
96  std::move(has_converged));
97  }
98 };
99 
100 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
101 struct ComputeAlpha {
102  private:
103  using fields_tag = FieldsTag;
104  using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
106 
107  public:
108  template <typename ParallelComponent, typename DbTagsList,
109  typename Metavariables, typename ArrayIndex,
110  typename DataBox = db::DataBox<DbTagsList>,
112  nullptr>
113  static void apply(db::DataBox<DbTagsList>& box,
115  const ArrayIndex& /*array_index*/,
116  const size_t iteration_id,
117  const double conj_grad_inner_product) noexcept {
118  Parallel::receive_data<Tags::Alpha<OptionsGroup>>(
119  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
120  get<residual_square_tag>(box) / conj_grad_inner_product);
121  }
122 };
123 
124 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
125 struct UpdateResidual {
126  private:
127  using fields_tag = FieldsTag;
128  using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
130  using initial_residual_magnitude_tag =
133 
134  public:
135  template <typename ParallelComponent, typename DbTagsList,
136  typename Metavariables, typename ArrayIndex,
137  typename DataBox = db::DataBox<DbTagsList>,
139  nullptr>
140  static void apply(db::DataBox<DbTagsList>& box,
142  const ArrayIndex& /*array_index*/,
143  const size_t iteration_id,
144  const double residual_square) noexcept {
145  // Compute the residual ratio before mutating the DataBox
146  const double res_ratio = residual_square / get<residual_square_tag>(box);
147 
148  db::mutate<residual_square_tag>(
149  make_not_null(&box),
150  [residual_square](
151  const gsl::not_null<double*> local_residual_square) noexcept {
152  *local_residual_square = residual_square;
153  });
154 
155  // At this point, the iteration is complete. We proceed with observing,
156  // logging and checking convergence before broadcasting back to the
157  // elements.
158 
159  const size_t completed_iterations = iteration_id + 1;
160  const double residual_magnitude = sqrt(residual_square);
161  LinearSolver::observe_detail::contribute_to_reduction_observer<
162  OptionsGroup>(completed_iterations, residual_magnitude, cache);
163 
164  // Determine whether the linear solver has converged
165  Convergence::HasConverged has_converged{
166  get<Convergence::Tags::Criteria<OptionsGroup>>(box),
167  completed_iterations, residual_magnitude,
168  get<initial_residual_magnitude_tag>(box)};
169 
170  // Do some logging
172  ::Verbosity::Verbose)) {
173  Parallel::printf("Linear solver '" +
174  Options::name<OptionsGroup>() +
175  "' iteration %zu done. Remaining residual: %e\n",
176  completed_iterations, residual_magnitude);
177  }
178  if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
179  cache) >= ::Verbosity::Quiet)) {
180  Parallel::printf("The linear solver '" +
181  Options::name<OptionsGroup>() +
182  "' has converged in %zu iterations: %s\n",
183  completed_iterations, has_converged);
184  }
185 
186  Parallel::receive_data<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(
187  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
188  // NOLINTNEXTLINE(performance-move-const-arg)
189  std::make_tuple(res_ratio, std::move(has_converged)));
190  }
191 };
192 
193 } // namespace LinearSolver::cg::detail
std::apply
T apply(T... args)
LinearSolver::Tags::MagnitudeSquare
The magnitude square w.r.t. the LinearSolver::inner_product
Definition: Tags.hpp:100
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
Options.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
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:68
Tags.hpp
std::sqrt
T sqrt(T... args)
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