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 <string>
9 #include <tuple>
10 #include <utility>
11 #include <vector>
12 
14 #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 #include "IO/Observer/Actions/RegisterWithObservers.hpp"
16 #include "IO/Observer/ArrayComponentId.hpp"
17 #include "IO/Observer/Helpers.hpp"
18 #include "IO/Observer/ObservationId.hpp"
19 #include "IO/Observer/ObserverComponent.hpp"
20 #include "IO/Observer/Protocols/ReductionDataFormatter.hpp"
21 #include "IO/Observer/ReductionActions.hpp"
22 #include "IO/Observer/TypeOfObservation.hpp"
23 #include "Informer/Tags.hpp"
24 #include "Informer/Verbosity.hpp"
25 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
26 #include "NumericalAlgorithms/Convergence/Tags.hpp"
27 #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
29 #include "Parallel/GlobalCache.hpp"
30 #include "Parallel/Invoke.hpp"
31 #include "Parallel/Printf.hpp"
32 #include "Parallel/Reduction.hpp"
33 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
35 #include "Utilities/Functional.hpp"
36 #include "Utilities/GetOutput.hpp"
37 #include "Utilities/Gsl.hpp"
38 #include "Utilities/PrettyType.hpp"
39 #include "Utilities/ProtocolHelpers.hpp"
40 #include "Utilities/Requires.hpp"
41 #include "Utilities/TMPL.hpp"
42 
43 /// \cond
44 namespace tuples {
45 template <typename...>
46 class TaggedTuple;
47 } // namespace tuples
49 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
50  typename Label>
51 struct CompleteStep;
52 } // namespace LinearSolver::async_solvers
53 /// \endcond
54 
55 /// Functionality shared between parallel linear solvers that have no global
56 /// synchronization points
58 
59 using reduction_data = Parallel::ReductionData<
60  // Iteration
62  // Residual
64 
65 template <typename OptionsGroup>
67  : tt::ConformsTo<observers::protocols::ReductionDataFormatter> {
68  using reduction_data = async_solvers::reduction_data;
69  std::string operator()(const size_t iteration_id, const double residual) const
70  noexcept {
71  if (iteration_id == 0) {
72  return Options::name<OptionsGroup>() +
73  " initialized with residual: " + get_output(residual);
74  } else {
75  return Options::name<OptionsGroup>() + "(" + get_output(iteration_id) +
76  ") iteration complete. Remaining residual: " +
77  get_output(residual);
78  }
79  }
80  // NOLINTNEXTLINE(google-runtime-references)
81  void pup(PUP::er& /*p*/) noexcept {}
82 };
83 
84 template <typename OptionsGroup, typename ParallelComponent,
85  typename Metavariables, typename ArrayIndex>
86 void contribute_to_residual_observation(
87  const size_t iteration_id, const double residual_magnitude_square,
89  const ArrayIndex& array_index) noexcept {
90  auto& local_observer =
91  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
92  cache)
93  .ckLocalBranch();
94  auto formatter =
96  ::Verbosity::Quiet)
97  ? std::make_optional(ResidualReductionFormatter<OptionsGroup>{})
98  : std::nullopt;
99  Parallel::simple_action<observers::Actions::ContributeReductionData>(
100  local_observer,
101  observers::ObservationId(iteration_id,
102  pretty_type::get_name<OptionsGroup>()),
105  Parallel::ArrayIndex<ArrayIndex>(array_index)},
106  std::string{"/" + Options::name<OptionsGroup>() + "Residuals"},
107  std::vector<std::string>{"Iteration", "Residual"},
108  reduction_data{iteration_id, residual_magnitude_square},
109  std::move(formatter));
111  ::Verbosity::Debug)) {
112  if (iteration_id == 0) {
113  Parallel::printf("%s %s initialized with local residual: %e\n",
114  get_output(array_index), Options::name<OptionsGroup>(),
115  sqrt(residual_magnitude_square));
116  } else {
118  "%s %s(%zu) iteration complete. Remaining local residual: %e\n",
119  get_output(array_index), Options::name<OptionsGroup>(), iteration_id,
120  sqrt(residual_magnitude_square));
121  }
122  }
123 }
124 
125 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
127  private:
128  using fields_tag = FieldsTag;
129  using operator_applied_to_fields_tag =
131  using source_tag = SourceTag;
132  using residual_tag =
136 
137  public:
138  using const_global_cache_tags =
139  tmpl::list<Convergence::Tags::Iterations<OptionsGroup>>;
140 
141  using simple_tags = tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
143  operator_applied_to_fields_tag>;
144  using compute_tags =
145  tmpl::list<LinearSolver::Tags::ResidualCompute<fields_tag, source_tag>>;
146 
147  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
148  typename ArrayIndex, typename ActionList,
149  typename ParallelComponent>
150  static auto apply(db::DataBox<DbTagsList>& box,
151  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
152  const Parallel::GlobalCache<Metavariables>& /*cache*/,
153  const ArrayIndex& /*array_index*/,
154  const ActionList /*meta*/,
155  const ParallelComponent* const /*meta*/) noexcept {
156  // The `PrepareSolve` action populates these tags with initial
157  // values, except for `operator_applied_to_fields_tag` which is
158  // expected to be updated in every iteration of the algorithm
160  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>>>(
162  return std::make_tuple(std::move(box));
163  }
164 };
165 
166 template <typename OptionsGroup>
168  template <typename ParallelComponent, typename DbTagsList,
169  typename ArrayIndex>
171  register_info(const db::DataBox<DbTagsList>& /*box*/,
172  const ArrayIndex& /*array_index*/) noexcept {
174  observers::ObservationKey{pretty_type::get_name<OptionsGroup>()}};
175  }
176 };
177 
178 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
179 using RegisterElement =
181 
182 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
183  typename Label>
184 struct PrepareSolve {
185  private:
186  using fields_tag = FieldsTag;
187  using residual_tag =
189 
190  public:
191  using const_global_cache_tags =
192  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
193 
194  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
195  typename ArrayIndex, typename ActionList,
196  typename ParallelComponent>
197  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
198  db::DataBox<DbTagsList>& box,
199  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
201  const ArrayIndex& array_index, const ActionList /*meta*/,
202  const ParallelComponent* const /*meta*/) noexcept {
203  constexpr size_t iteration_id = 0;
204 
205  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
207  make_not_null(&box),
208  [](const gsl::not_null<size_t*> local_iteration_id,
209  const gsl::not_null<Convergence::HasConverged*> has_converged,
210  const size_t num_iterations) noexcept {
211  *local_iteration_id = iteration_id;
212  *has_converged =
213  Convergence::HasConverged{num_iterations, iteration_id};
214  },
215  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
216 
217  // Observe the initial residual even if no steps are going to be performed
218  const auto& residual = get<residual_tag>(box);
219  const double residual_magnitude_square = inner_product(residual, residual);
220  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
221  iteration_id, residual_magnitude_square, cache, array_index);
222 
223  // Skip steps entirely if the solve has already converged
224  constexpr size_t step_end_index =
225  tmpl::index_of<ActionList, CompleteStep<FieldsTag, OptionsGroup,
226  SourceTag, Label>>::value;
227  constexpr size_t this_action_index =
228  tmpl::index_of<ActionList, PrepareSolve>::value;
229  return {std::move(box), false,
230  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
231  ? (step_end_index + 1)
232  : (this_action_index + 1)};
233  }
234 };
235 
236 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
237  typename Label>
238 struct CompleteStep {
239  private:
240  using fields_tag = FieldsTag;
241  using residual_tag =
243 
244  public:
245  using const_global_cache_tags =
246  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
247 
248  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
249  typename ArrayIndex, typename ActionList,
250  typename ParallelComponent>
251  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
252  db::DataBox<DbTagsList>& box,
255  const ArrayIndex& array_index, const ActionList /*meta*/,
256  const ParallelComponent* const /*meta*/) noexcept {
257  // Prepare for next iteration
258  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
260  make_not_null(&box),
261  [](const gsl::not_null<size_t*> iteration_id,
262  const gsl::not_null<Convergence::HasConverged*> has_converged,
263  const size_t num_iterations) noexcept {
264  ++(*iteration_id);
265  *has_converged =
266  Convergence::HasConverged{num_iterations, *iteration_id};
267  },
268  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
269 
270  // Observe element-local residual magnitude
271  const size_t completed_iterations =
272  get<Convergence::Tags::IterationId<OptionsGroup>>(box);
273  const auto& residual = get<residual_tag>(box);
274  const double residual_magnitude_square = inner_product(residual, residual);
275  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
276  completed_iterations, residual_magnitude_square, cache, array_index);
277 
278  // Repeat steps until the solve has converged
279  constexpr size_t step_begin_index =
280  tmpl::index_of<ActionList, PrepareSolve<FieldsTag, OptionsGroup,
281  SourceTag, Label>>::value +
282  1;
283  constexpr size_t this_action_index =
284  tmpl::index_of<ActionList, CompleteStep>::value;
285  return {std::move(box), false,
286  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
287  ? (this_action_index + 1)
288  : step_begin_index};
289  }
290 };
291 
292 } // 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:100
std::string
LinearSolver::async_solvers::ResidualReductionFormatter
Definition: ElementActions.hpp:66
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
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:63
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
LinearSolver::async_solvers::RegisterObservers
Definition: ElementActions.hpp:167
LinearSolver::async_solvers::InitializeElement
Definition: ElementActions.hpp:126
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:40
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
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
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: PhaseControlTags.hpp:92
limits
Gsl.hpp
LinearSolver::async_solvers
Functionality shared between parallel linear solvers that have no global synchronization points.
Definition: ElementActions.hpp:57
Requires.hpp
observers::TypeOfObservation::Reduction
@ Reduction
The sender will only perform reduction observations.
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
Definition: ElementActions.hpp:184
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: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
observers::Actions::RegisterWithObservers
Register an observation ID with the observers.
Definition: RegisterWithObservers.hpp:47
LinearSolver::async_solvers::CompleteStep
Definition: ElementActions.hpp:238
string