ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <limits>
8 #include <optional>
9 #include <string>
10 #include <tuple>
11 #include <utility>
12 #include <vector>
13 
15 #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 #include "IO/Logging/Tags.hpp"
17 #include "IO/Logging/Verbosity.hpp"
18 #include "IO/Observer/Actions/RegisterWithObservers.hpp"
19 #include "IO/Observer/ArrayComponentId.hpp"
20 #include "IO/Observer/GetSectionObservationKey.hpp"
21 #include "IO/Observer/Helpers.hpp"
22 #include "IO/Observer/ObservationId.hpp"
23 #include "IO/Observer/ObserverComponent.hpp"
24 #include "IO/Observer/Protocols/ReductionDataFormatter.hpp"
25 #include "IO/Observer/ReductionActions.hpp"
26 #include "IO/Observer/Tags.hpp"
27 #include "IO/Observer/TypeOfObservation.hpp"
28 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
29 #include "NumericalAlgorithms/Convergence/Tags.hpp"
30 #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
32 #include "Parallel/GlobalCache.hpp"
33 #include "Parallel/Invoke.hpp"
34 #include "Parallel/Printf.hpp"
35 #include "Parallel/Reduction.hpp"
36 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
38 #include "Utilities/Functional.hpp"
39 #include "Utilities/GetOutput.hpp"
40 #include "Utilities/Gsl.hpp"
41 #include "Utilities/PrettyType.hpp"
42 #include "Utilities/ProtocolHelpers.hpp"
43 #include "Utilities/Requires.hpp"
44 #include "Utilities/TMPL.hpp"
45 
46 /// \cond
47 namespace tuples {
48 template <typename...>
49 class TaggedTuple;
50 } // namespace tuples
52 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
53  typename Label, typename ArraySectionIdTag,
54  bool ObserveInitialResidual>
55 struct CompleteStep;
56 } // namespace LinearSolver::async_solvers
57 /// \endcond
58 
59 /// Functionality shared between parallel linear solvers that have no global
60 /// synchronization points
62 
63 using reduction_data = Parallel::ReductionData<
64  // Iteration
66  // Residual
68 
69 template <typename OptionsGroup>
71  : tt::ConformsTo<observers::protocols::ReductionDataFormatter> {
72  using reduction_data = async_solvers::reduction_data;
73  ResidualReductionFormatter() noexcept = default;
74  ResidualReductionFormatter(std::string local_section_observation_key) noexcept
75  : section_observation_key(std::move(local_section_observation_key)) {}
76  std::string operator()(const size_t iteration_id,
77  const double residual) const noexcept {
78  if (iteration_id == 0) {
79  return Options::name<OptionsGroup>() + section_observation_key +
80  " initialized with residual: " + get_output(residual);
81  } else {
82  return Options::name<OptionsGroup>() + section_observation_key + "(" +
83  get_output(iteration_id) +
84  ") iteration complete. Remaining residual: " +
85  get_output(residual);
86  }
87  }
88  // NOLINTNEXTLINE(google-runtime-references)
89  void pup(PUP::er& p) noexcept { p | section_observation_key; }
90  std::string section_observation_key{};
91 };
92 
93 template <typename OptionsGroup, typename ParallelComponent,
94  typename Metavariables, typename ArrayIndex>
95 void contribute_to_residual_observation(
96  const size_t iteration_id, const double residual_magnitude_square,
97  Parallel::GlobalCache<Metavariables>& cache, const ArrayIndex& array_index,
98  const std::string& section_observation_key) noexcept {
99  auto& local_observer =
100  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
101  cache)
102  .ckLocalBranch();
103  auto formatter =
105  ::Verbosity::Quiet)
106  ? std::make_optional(ResidualReductionFormatter<OptionsGroup>{
107  section_observation_key})
108  : std::nullopt;
109  Parallel::simple_action<observers::Actions::ContributeReductionData>(
110  local_observer,
112  iteration_id,
113  pretty_type::get_name<OptionsGroup>() + section_observation_key),
116  Parallel::ArrayIndex<ArrayIndex>(array_index)},
117  std::string{"/" + Options::name<OptionsGroup>() +
118  section_observation_key + "Residuals"},
119  std::vector<std::string>{"Iteration", "Residual"},
120  reduction_data{iteration_id, residual_magnitude_square},
121  std::move(formatter));
123  ::Verbosity::Debug)) {
124  if (iteration_id == 0) {
125  Parallel::printf("%s %s initialized with local residual: %e\n",
126  get_output(array_index),
127  Options::name<OptionsGroup>() + section_observation_key,
128  sqrt(residual_magnitude_square));
129  } else {
131  "%s %s(%zu) iteration complete. Remaining local residual: %e\n",
132  get_output(array_index),
133  Options::name<OptionsGroup>() + section_observation_key, iteration_id,
134  sqrt(residual_magnitude_square));
135  }
136  }
137 }
138 
139 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
141  private:
142  using fields_tag = FieldsTag;
143  using operator_applied_to_fields_tag =
145  using source_tag = SourceTag;
146  using residual_tag =
150 
151  public:
152  using const_global_cache_tags =
153  tmpl::list<Convergence::Tags::Iterations<OptionsGroup>>;
154 
155  using simple_tags = tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
157  operator_applied_to_fields_tag>;
158  using compute_tags =
159  tmpl::list<LinearSolver::Tags::ResidualCompute<fields_tag, source_tag>>;
160 
161  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
162  typename ArrayIndex, typename ActionList,
163  typename ParallelComponent>
164  static auto apply(db::DataBox<DbTagsList>& box,
165  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
166  const Parallel::GlobalCache<Metavariables>& /*cache*/,
167  const ArrayIndex& /*array_index*/,
168  const ActionList /*meta*/,
169  const ParallelComponent* const /*meta*/) noexcept {
170  // The `PrepareSolve` action populates these tags with initial
171  // values, except for `operator_applied_to_fields_tag` which is
172  // expected to be updated in every iteration of the algorithm
174  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>>>(
176  return std::make_tuple(std::move(box));
177  }
178 };
179 
180 template <typename OptionsGroup, typename ArraySectionIdTag = void>
182  template <typename ParallelComponent, typename DbTagsList,
183  typename ArrayIndex>
185  register_info(const db::DataBox<DbTagsList>& box,
186  const ArrayIndex& /*array_index*/) noexcept {
187  // Get the observation key, or "Unused" if the element does not belong
188  // to a section with this tag. In the latter case, no observations will
189  // ever be contributed.
190  const std::optional<std::string>& section_observation_key =
191  observers::get_section_observation_key<ArraySectionIdTag>(box);
192  ASSERT(section_observation_key != "Unused",
193  "The identifier 'Unused' is reserved to indicate that no "
194  "observations with this key will be contributed. Use a different "
195  "key, or change the identifier 'Unused' to something else.");
196  return {
198  observers::ObservationKey{pretty_type::get_name<OptionsGroup>() +
199  section_observation_key.value_or("Unused")}};
200  }
201 };
202 
203 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
204  typename ArraySectionIdTag = void>
207 
208 /*!
209  * \brief Prepare the asynchronous linear solver for a solve
210  *
211  * This action resets the asynchronous linear solver to its initial state, and
212  * optionally observes the initial residual. If the initial residual should be
213  * observed, both the `SourceTag` as well as the
214  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, FieldsTag>`
215  * must be up-to-date at the time this action is invoked.
216  *
217  * This action also provides an anchor point in the action list for looping the
218  * linear solver, in the sense that the algorithm jumps back to the action
219  * immediately following this one when a step is complete and the solver hasn't
220  * yet converged (see `LinearSolver::async_solvers::CompleteStep`).
221  *
222  * \tparam FieldsTag The data `x` in the linear equation `Ax = b` to be solved.
223  * Should hold the initial guess `x_0` at this point in the algorithm.
224  * \tparam OptionsGroup An options group identifying the linear solver
225  * \tparam SourceTag The data `b` in `Ax = b`
226  * \tparam Label An optional compile-time label for the solver to distinguish
227  * different solves with the same solver in the action list
228  * \tparam ArraySectionIdTag Observe the residual norm separately for each
229  * array section identified by this tag (see `Parallel::Section`). Set to `void`
230  * to observe the residual norm over all elements of the array (default). The
231  * section only affects observations of residuals and has no effect on the
232  * solver algorithm.
233  * \tparam ObserveInitialResidual Whether or not to observe the initial residual
234  * `b - A x_0` (default: `true`). Disable when `b` or `A x_0` are not yet
235  * available at the preparation stage.
236  */
237 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
238  typename Label = OptionsGroup, typename ArraySectionIdTag = void,
239  bool ObserveInitialResidual = true>
240 struct PrepareSolve {
241  private:
242  using fields_tag = FieldsTag;
243  using residual_tag =
245 
246  public:
247  using const_global_cache_tags =
248  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
249 
250  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
251  typename ArrayIndex, typename ActionList,
252  typename ParallelComponent>
253  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
254  db::DataBox<DbTagsList>& box,
255  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
257  const ArrayIndex& array_index, const ActionList /*meta*/,
258  const ParallelComponent* const /*meta*/) noexcept {
259  constexpr size_t iteration_id = 0;
260 
262  ::Verbosity::Debug)) {
263  Parallel::printf("%s %s: Prepare solve\n", get_output(array_index),
264  Options::name<OptionsGroup>());
265  }
266 
267  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
269  make_not_null(&box),
270  [](const gsl::not_null<size_t*> local_iteration_id,
271  const gsl::not_null<Convergence::HasConverged*> has_converged,
272  const size_t num_iterations) noexcept {
273  *local_iteration_id = iteration_id;
274  *has_converged =
275  Convergence::HasConverged{num_iterations, iteration_id};
276  },
277  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
278 
279  if constexpr (ObserveInitialResidual) {
280  // Observe the initial residual even if no steps are going to be performed
281  const std::optional<std::string> section_observation_key =
282  observers::get_section_observation_key<ArraySectionIdTag>(box);
283  if (section_observation_key.has_value()) {
284  const auto& residual = get<residual_tag>(box);
285  const double residual_magnitude_square =
286  inner_product(residual, residual);
287  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
288  iteration_id, residual_magnitude_square, cache, array_index,
289  *section_observation_key);
290  }
291  }
292 
293  // Skip steps entirely if the solve has already converged
294  constexpr size_t step_end_index = tmpl::index_of<
295  ActionList,
296  CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
297  ArraySectionIdTag, ObserveInitialResidual>>::value;
298  constexpr size_t this_action_index =
299  tmpl::index_of<ActionList, PrepareSolve>::value;
300  return {std::move(box), false,
301  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
302  ? (step_end_index + 1)
303  : (this_action_index + 1)};
304  }
305 };
306 
307 /*!
308  * \brief Complete a step of the asynchronous linear solver
309  *
310  * This action prepares the next step of the asynchronous linear solver, and
311  * observes the residual. To observe the correct residual, make sure the
312  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, FieldsTag>` is
313  * up-to-date at the time this action is invoked.
314  *
315  * This action checks if the algorithm has converged, i.e. it has completed the
316  * requested number of steps. If it hasn't, the algorithm jumps back to the
317  * action immediately following the `LinearSolver::async_solvers::PrepareSolve`
318  * to perform another iteration. Make sure both actions use the same template
319  * parameters.
320  *
321  * \tparam FieldsTag The data `x` in the linear equation `Ax = b` to be solved.
322  * \tparam OptionsGroup An options group identifying the linear solver
323  * \tparam SourceTag The data `b` in `Ax = b`
324  * \tparam Label An optional compile-time label for the solver to distinguish
325  * different solves with the same solver in the action list
326  * \tparam ArraySectionIdTag Observe the residual norm separately for each
327  * array section identified by this tag (see `Parallel::Section`). Set to `void`
328  * to observe the residual norm over all elements of the array (default). The
329  * section only affects observations of residuals and has no effect on the
330  * solver algorithm.
331  * \tparam ObserveInitialResidual Whether or not to observe the _initial_
332  * residual `b - A x_0`. This parameter should match the one passed to
333  * `PrepareSolve`.
334  */
335 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
336  typename Label = OptionsGroup, typename ArraySectionIdTag = void,
337  bool ObserveInitialResidual = true>
338 struct CompleteStep {
339  private:
340  using fields_tag = FieldsTag;
341  using residual_tag =
343 
344  public:
345  using const_global_cache_tags =
346  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
347 
348  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
349  typename ArrayIndex, typename ActionList,
350  typename ParallelComponent>
351  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
352  db::DataBox<DbTagsList>& box,
355  const ArrayIndex& array_index, const ActionList /*meta*/,
356  const ParallelComponent* const /*meta*/) noexcept {
357  // Prepare for next iteration
358  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
360  make_not_null(&box),
361  [](const gsl::not_null<size_t*> iteration_id,
362  const gsl::not_null<Convergence::HasConverged*> has_converged,
363  const size_t num_iterations) noexcept {
364  ++(*iteration_id);
365  *has_converged =
366  Convergence::HasConverged{num_iterations, *iteration_id};
367  },
368  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
369 
370  // Observe element-local residual magnitude
371  const std::optional<std::string> section_observation_key =
372  observers::get_section_observation_key<ArraySectionIdTag>(box);
373  if (section_observation_key.has_value()) {
374  const size_t completed_iterations =
375  get<Convergence::Tags::IterationId<OptionsGroup>>(box);
376  const auto& residual = get<residual_tag>(box);
377  const double residual_magnitude_square =
378  inner_product(residual, residual);
379  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
380  completed_iterations, residual_magnitude_square, cache, array_index,
381  *section_observation_key);
382  }
383 
384  // Repeat steps until the solve has converged
385  constexpr size_t step_begin_index =
386  tmpl::index_of<
387  ActionList,
388  PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label,
389  ArraySectionIdTag, ObserveInitialResidual>>::value +
390  1;
391  constexpr size_t this_action_index =
392  tmpl::index_of<ActionList, CompleteStep>::value;
393  return {std::move(box), false,
394  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
395  ? (this_action_index + 1)
396  : step_begin_index};
397  }
398 };
399 
400 } // namespace LinearSolver::async_solvers
observers::ObservationId
A unique identifier for an observation representing the type of observation and the instance (e....
Definition: ObservationId.hpp:71
LinearSolver::Tags::MagnitudeSquare
The magnitude square w.r.t. the LinearSolver::inner_product
Definition: Tags.hpp:99
std::string
LinearSolver::async_solvers::ResidualReductionFormatter
Definition: ElementActions.hpp:70
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
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:83
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:65
InnerProduct.hpp
std::pair
GlobalCache.hpp
vector
PrettyType.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
observers::ObservationKey
Used as a key in maps to keep track of how many elements have registered.
Definition: ObservationId.hpp:28
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
LinearSolver::async_solvers::RegisterObservers
Definition: ElementActions.hpp:181
LinearSolver::async_solvers::InitializeElement
Definition: ElementActions.hpp:140
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:68
Tags.hpp
get_output
std::string get_output(const T &t) noexcept
Get the streamed output of t as a std::string
Definition: GetOutput.hpp:14
Initialization::mutate_assign
constexpr void mutate_assign(const gsl::not_null< db::DataBox< BoxTags > * > box, Args &&... args) noexcept
Perform a mutation to the DataBox box, assigning the args to the tags in MutateTagList in order.
Definition: MutateAssign.hpp:39
funcl::Sqrt
Functional for computing sqrt on an object.
Definition: Functional.hpp:275
std::add_pointer_t
Printf.hpp
LinearSolver::inner_product
double inner_product(const Lhs &lhs, const Rhs &rhs) noexcept
The local part of the Euclidean inner product on the vector space w.r.t. which the addition and scala...
Definition: InnerProduct.hpp:71
DataBox.hpp
cstddef
logging::Tags::Verbosity
Tag for putting Verbosity in a DataBox.
Definition: Tags.hpp:33
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
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:382
PhaseControl::reduction_data
Parallel::ReductionData< Parallel::ReductionDatum< tuples::tagged_tuple_from_typelist< TagsAndCombinesPresent >, TaggedTupleCombine > > reduction_data
A Parallel::ReductionData with a single Parallel::ReductionDatum for a given tagged tuple type determ...
Definition: PhaseControlReductionHelpers.hpp:83
limits
Gsl.hpp
LinearSolver::async_solvers
Functionality shared between parallel linear solvers that have no global synchronization points.
Definition: ElementActions.hpp:61
Requires.hpp
observers::TypeOfObservation::Reduction
@ Reduction
The sender will only perform reduction observations.
optional
observers::ArrayComponentId
An ID type that identifies both the parallel component and the index in the parallel component.
Definition: ArrayComponentId.hpp:27
LinearSolver::async_solvers::PrepareSolve
Prepare the asynchronous linear solver for a solve.
Definition: ElementActions.hpp:240
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
Parallel::ArrayIndex
The array index used for indexing Chare Arrays, mostly an implementation detail.
Definition: ArrayIndex.hpp:28
std::numeric_limits::max
T max(T... args)
Convergence::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the iterative algorithm has converged,...
Definition: Tags.hpp:86
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
observers::Actions::RegisterWithObservers
Register an observation ID with the observers.
Definition: RegisterWithObservers.hpp:47
LinearSolver::async_solvers::CompleteStep
Complete a step of the asynchronous linear solver.
Definition: ElementActions.hpp:338
string