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, ParallelComponent>(iteration_id, residual_magnitude,
71  cache);
72 
73  // Determine whether the linear solver has converged
74  Convergence::HasConverged has_converged{
75  get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
76  residual_magnitude, residual_magnitude};
77 
78  // Do some logging
80  ::Verbosity::Quiet)) {
81  Parallel::printf("%s initialized with residual: %e\n",
82  Options::name<OptionsGroup>(), residual_magnitude);
83  }
85  cache) >= ::Verbosity::Quiet)) {
87  "%s has converged without any iterations: %s\n",
88  Options::name<OptionsGroup>(), has_converged);
89  }
90 
91  Parallel::receive_data<Tags::InitialHasConverged<OptionsGroup>>(
92  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
93  // NOLINTNEXTLINE(performance-move-const-arg)
94  std::move(has_converged));
95  }
96 };
97 
98 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
99 struct ComputeAlpha {
100  private:
101  using fields_tag = FieldsTag;
102  using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
104 
105  public:
106  template <typename ParallelComponent, typename DbTagsList,
107  typename Metavariables, typename ArrayIndex,
108  typename DataBox = db::DataBox<DbTagsList>,
110  nullptr>
111  static void apply(db::DataBox<DbTagsList>& box,
113  const ArrayIndex& /*array_index*/,
114  const size_t iteration_id,
115  const double conj_grad_inner_product) noexcept {
116  Parallel::receive_data<Tags::Alpha<OptionsGroup>>(
117  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
118  get<residual_square_tag>(box) / conj_grad_inner_product);
119  }
120 };
121 
122 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
123 struct UpdateResidual {
124  private:
125  using fields_tag = FieldsTag;
126  using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
128  using initial_residual_magnitude_tag =
131 
132  public:
133  template <typename ParallelComponent, typename DbTagsList,
134  typename Metavariables, typename ArrayIndex,
135  typename DataBox = db::DataBox<DbTagsList>,
137  nullptr>
138  static void apply(db::DataBox<DbTagsList>& box,
140  const ArrayIndex& /*array_index*/,
141  const size_t iteration_id,
142  const double residual_square) noexcept {
143  // Compute the residual ratio before mutating the DataBox
144  const double res_ratio = residual_square / get<residual_square_tag>(box);
145 
146  db::mutate<residual_square_tag>(
147  make_not_null(&box),
148  [residual_square](
149  const gsl::not_null<double*> local_residual_square) noexcept {
150  *local_residual_square = residual_square;
151  });
152 
153  // At this point, the iteration is complete. We proceed with observing,
154  // logging and checking convergence before broadcasting back to the
155  // elements.
156 
157  const size_t completed_iterations = iteration_id + 1;
158  const double residual_magnitude = sqrt(residual_square);
159  LinearSolver::observe_detail::contribute_to_reduction_observer<
160  OptionsGroup, ParallelComponent>(completed_iterations,
161  residual_magnitude, cache);
162 
163  // Determine whether the linear solver has converged
164  Convergence::HasConverged has_converged{
165  get<Convergence::Tags::Criteria<OptionsGroup>>(box),
166  completed_iterations, residual_magnitude,
167  get<initial_residual_magnitude_tag>(box)};
168 
169  // Do some logging
171  ::Verbosity::Quiet)) {
173  "%s(%zu) iteration complete. Remaining residual: %e\n",
174  Options::name<OptionsGroup>(), completed_iterations,
175  residual_magnitude);
176  }
177  if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
178  cache) >= ::Verbosity::Quiet)) {
180  "%s has converged in %zu iterations: %s\n",
181  Options::name<OptionsGroup>(), completed_iterations, has_converged);
182  }
183 
184  Parallel::receive_data<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(
185  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
186  // NOLINTNEXTLINE(performance-move-const-arg)
187  std::make_tuple(res_ratio, std::move(has_converged)));
188  }
189 };
190 
191 } // 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: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
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
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