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/ReductionActions.hpp"
21 #include "IO/Observer/TypeOfObservation.hpp"
22 #include "Informer/Tags.hpp"
23 #include "Informer/Verbosity.hpp"
24 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
25 #include "NumericalAlgorithms/Convergence/Tags.hpp"
26 #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
28 #include "Parallel/GlobalCache.hpp"
29 #include "Parallel/Invoke.hpp"
30 #include "Parallel/Printf.hpp"
31 #include "Parallel/Reduction.hpp"
32 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
34 #include "Utilities/Functional.hpp"
35 #include "Utilities/GetOutput.hpp"
36 #include "Utilities/Gsl.hpp"
37 #include "Utilities/PrettyType.hpp"
38 #include "Utilities/Requires.hpp"
39 #include "Utilities/TMPL.hpp"
40 
41 /// \cond
42 namespace tuples {
43 template <typename...>
44 class TaggedTuple;
45 } // namespace tuples
47 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
48  typename Label>
49 struct CompleteStep;
50 } // namespace LinearSolver::async_solvers
51 /// \endcond
52 
53 /// Functionality shared between parallel linear solvers that have no global
54 /// synchronization points
56 
57 using reduction_data = Parallel::ReductionData<
58  // Iteration
60  // Residual
62 
63 template <typename OptionsGroup, typename ParallelComponent,
64  typename Metavariables, typename ArrayIndex>
65 void contribute_to_residual_observation(
66  const size_t iteration_id, const double residual_magnitude_square,
68  const ArrayIndex& array_index) noexcept {
69  auto& local_observer =
70  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
71  cache)
72  .ckLocalBranch();
73  Parallel::simple_action<observers::Actions::ContributeReductionData>(
74  local_observer,
75  observers::ObservationId(iteration_id,
76  pretty_type::get_name<OptionsGroup>()),
80  std::string{"/" + Options::name<OptionsGroup>() + "Residuals"},
81  std::vector<std::string>{"Iteration", "Residual"},
82  reduction_data{iteration_id, residual_magnitude_square});
84  ::Verbosity::Verbose)) {
85  if (iteration_id == 0) {
87  "Linear solver '" + Options::name<OptionsGroup>() +
88  "' initialized on element %s. Remaining local residual: %e\n",
89  get_output(array_index), sqrt(residual_magnitude_square));
90  } else {
91  Parallel::printf("Linear solver '" +
92  Options::name<OptionsGroup>() +
93  "' iteration %zu done on element %s. Remaining "
94  "local residual: %e\n",
95  iteration_id, get_output(array_index),
96  sqrt(residual_magnitude_square));
97  }
98  }
99 }
100 
101 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
103  private:
104  using fields_tag = FieldsTag;
105  using operator_applied_to_fields_tag =
107  using source_tag = SourceTag;
108  using residual_tag =
112 
113  public:
114  using const_global_cache_tags =
115  tmpl::list<Convergence::Tags::Iterations<OptionsGroup>>;
116 
117  using simple_tags = tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
119  operator_applied_to_fields_tag>;
120  using compute_tags =
121  tmpl::list<LinearSolver::Tags::ResidualCompute<fields_tag, source_tag>>;
122 
123  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
124  typename ArrayIndex, typename ActionList,
125  typename ParallelComponent>
126  static auto apply(db::DataBox<DbTagsList>& box,
127  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
128  const Parallel::GlobalCache<Metavariables>& /*cache*/,
129  const ArrayIndex& /*array_index*/,
130  const ActionList /*meta*/,
131  const ParallelComponent* const /*meta*/) noexcept {
132  // The `PrepareSolve` action populates these tags with initial
133  // values, except for `operator_applied_to_fields_tag` which is
134  // expected to be updated in every iteration of the algorithm
136  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>>>(
138  return std::make_tuple(std::move(box));
139  }
140 };
141 
142 template <typename OptionsGroup>
144  template <typename ParallelComponent, typename DbTagsList,
145  typename ArrayIndex>
147  register_info(const db::DataBox<DbTagsList>& /*box*/,
148  const ArrayIndex& /*array_index*/) noexcept {
150  observers::ObservationKey{pretty_type::get_name<OptionsGroup>()}};
151  }
152 };
153 
154 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
155 using RegisterElement =
157 
158 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
159  typename Label>
160 struct PrepareSolve {
161  private:
162  using fields_tag = FieldsTag;
163  using residual_tag =
165 
166  public:
167  using const_global_cache_tags =
168  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
169 
170  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
171  typename ArrayIndex, typename ActionList,
172  typename ParallelComponent>
173  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
174  db::DataBox<DbTagsList>& box,
175  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
177  const ArrayIndex& array_index, const ActionList /*meta*/,
178  const ParallelComponent* const /*meta*/) noexcept {
179  constexpr size_t iteration_id = 0;
180 
181  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
183  make_not_null(&box),
184  [](const gsl::not_null<size_t*> local_iteration_id,
185  const gsl::not_null<Convergence::HasConverged*> has_converged,
186  const size_t num_iterations) noexcept {
187  *local_iteration_id = iteration_id;
188  *has_converged =
189  Convergence::HasConverged{num_iterations, iteration_id};
190  },
191  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
192 
193  // Observe the initial residual even if no steps are going to be performed
194  const auto& residual = get<residual_tag>(box);
195  const double residual_magnitude_square = inner_product(residual, residual);
196  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
197  iteration_id, residual_magnitude_square, cache, array_index);
198 
199  // Skip steps entirely if the solve has already converged
200  constexpr size_t step_end_index =
201  tmpl::index_of<ActionList, CompleteStep<FieldsTag, OptionsGroup,
202  SourceTag, Label>>::value;
203  constexpr size_t this_action_index =
204  tmpl::index_of<ActionList, PrepareSolve>::value;
205  return {std::move(box), false,
206  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
207  ? (step_end_index + 1)
208  : (this_action_index + 1)};
209  }
210 };
211 
212 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
213  typename Label>
214 struct CompleteStep {
215  private:
216  using fields_tag = FieldsTag;
217  using residual_tag =
219 
220  public:
221  using const_global_cache_tags =
222  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
223 
224  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
225  typename ArrayIndex, typename ActionList,
226  typename ParallelComponent>
227  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
228  db::DataBox<DbTagsList>& box,
231  const ArrayIndex& array_index, const ActionList /*meta*/,
232  const ParallelComponent* const /*meta*/) noexcept {
233  // Prepare for next iteration
234  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
236  make_not_null(&box),
237  [](const gsl::not_null<size_t*> iteration_id,
238  const gsl::not_null<Convergence::HasConverged*> has_converged,
239  const size_t num_iterations) noexcept {
240  ++(*iteration_id);
241  *has_converged =
242  Convergence::HasConverged{num_iterations, *iteration_id};
243  },
244  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
245 
246  // Observe element-local residual magnitude
247  const size_t completed_iterations =
248  get<Convergence::Tags::IterationId<OptionsGroup>>(box);
249  const auto& residual = get<residual_tag>(box);
250  const double residual_magnitude_square = inner_product(residual, residual);
251  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
252  completed_iterations, residual_magnitude_square, cache, array_index);
253 
254  // Repeat steps until the solve has converged
255  constexpr size_t step_begin_index =
256  tmpl::index_of<ActionList, PrepareSolve<FieldsTag, OptionsGroup,
257  SourceTag, Label>>::value +
258  1;
259  constexpr size_t this_action_index =
260  tmpl::index_of<ActionList, CompleteStep>::value;
261  return {std::move(box), false,
262  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
263  ? (this_action_index + 1)
264  : step_begin_index};
265  }
266 };
267 
268 } // 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
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
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:143
LinearSolver::async_solvers::InitializeElement
Definition: ElementActions.hpp:102
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:273
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
limits
Gsl.hpp
LinearSolver::async_solvers
Functionality shared between parallel linear solvers that have no global synchronization points.
Definition: ElementActions.hpp:55
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:160
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
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:39
LinearSolver::async_solvers::CompleteStep
Definition: ElementActions.hpp:214
string