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/MergeIntoDataBox.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  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
118  typename ArrayIndex, typename ActionList,
119  typename ParallelComponent>
120  static auto apply(db::DataBox<DbTagsList>& box,
121  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
122  const Parallel::GlobalCache<Metavariables>& /*cache*/,
123  const ArrayIndex& /*array_index*/,
124  const ActionList /*meta*/,
125  const ParallelComponent* const /*meta*/) noexcept {
126  return std::make_tuple(
131  operator_applied_to_fields_tag>,
134  std::move(box),
135  // The `PrepareSolve` action populates these tags with initial
136  // values, except for `operator_applied_to_fields_tag` which is
137  // expected to be updated in every iteration of the algorithm
139  typename operator_applied_to_fields_tag::type{}));
140  }
141 };
142 
143 template <typename OptionsGroup>
145  template <typename ParallelComponent, typename DbTagsList,
146  typename ArrayIndex>
148  register_info(const db::DataBox<DbTagsList>& /*box*/,
149  const ArrayIndex& /*array_index*/) noexcept {
151  observers::ObservationKey{pretty_type::get_name<OptionsGroup>()}};
152  }
153 };
154 
155 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
156 using RegisterElement =
158 
159 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
160  typename Label>
161 struct PrepareSolve {
162  private:
163  using fields_tag = FieldsTag;
164  using residual_tag =
166 
167  public:
168  using const_global_cache_tags =
169  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
170 
171  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
172  typename ArrayIndex, typename ActionList,
173  typename ParallelComponent>
174  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
176  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
178  const ArrayIndex& array_index, const ActionList /*meta*/,
179  const ParallelComponent* const /*meta*/) noexcept {
180  constexpr size_t iteration_id = 0;
181 
182  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
184  make_not_null(&box),
185  [](const gsl::not_null<size_t*> local_iteration_id,
186  const gsl::not_null<Convergence::HasConverged*> has_converged,
187  const size_t num_iterations) noexcept {
188  *local_iteration_id = iteration_id;
189  *has_converged =
190  Convergence::HasConverged{num_iterations, iteration_id};
191  },
192  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
193 
194  // Observe the initial residual even if no steps are going to be performed
195  const auto& residual = get<residual_tag>(box);
196  const double residual_magnitude_square = inner_product(residual, residual);
197  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
198  iteration_id, residual_magnitude_square, cache, array_index);
199 
200  // Skip steps entirely if the solve has already converged
201  constexpr size_t step_end_index =
202  tmpl::index_of<ActionList, CompleteStep<FieldsTag, OptionsGroup,
203  SourceTag, Label>>::value;
204  constexpr size_t this_action_index =
205  tmpl::index_of<ActionList, PrepareSolve>::value;
206  return {std::move(box), false,
207  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
208  ? (step_end_index + 1)
209  : (this_action_index + 1)};
210  }
211 };
212 
213 template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
214  typename Label>
215 struct CompleteStep {
216  private:
217  using fields_tag = FieldsTag;
218  using residual_tag =
220 
221  public:
222  using const_global_cache_tags =
223  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
224 
225  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
226  typename ArrayIndex, typename ActionList,
227  typename ParallelComponent>
228  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
232  const ArrayIndex& array_index, const ActionList /*meta*/,
233  const ParallelComponent* const /*meta*/) noexcept {
234  // Prepare for next iteration
235  db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
237  make_not_null(&box),
238  [](const gsl::not_null<size_t*> iteration_id,
239  const gsl::not_null<Convergence::HasConverged*> has_converged,
240  const size_t num_iterations) noexcept {
241  ++(*iteration_id);
242  *has_converged =
243  Convergence::HasConverged{num_iterations, *iteration_id};
244  },
245  get<Convergence::Tags::Iterations<OptionsGroup>>(box));
246 
247  // Observe element-local residual magnitude
248  const size_t completed_iterations =
249  get<Convergence::Tags::IterationId<OptionsGroup>>(box);
250  const auto& residual = get<residual_tag>(box);
251  const double residual_magnitude_square = inner_product(residual, residual);
252  contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
253  completed_iterations, residual_magnitude_square, cache, array_index);
254 
255  // Repeat steps until the solve has converged
256  constexpr size_t step_begin_index =
257  tmpl::index_of<ActionList, PrepareSolve<FieldsTag, OptionsGroup,
258  SourceTag, Label>>::value +
259  1;
260  constexpr size_t this_action_index =
261  tmpl::index_of<ActionList, CompleteStep>::value;
262  return {std::move(box), false,
263  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
264  ? (this_action_index + 1)
265  : step_begin_index};
266  }
267 };
268 
269 } // 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:101
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:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
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:52
observers::ObservationKey
Used as a key in maps to keep track of how many elements have registered.
Definition: ObservationId.hpp:28
db::AddComputeTags
tmpl::flatten< tmpl::list< Tags... > > AddComputeTags
List of Compute Item Tags to add to the DataBox.
Definition: DataBox.hpp:1153
tuple
LinearSolver::async_solvers::RegisterObservers
Definition: ElementActions.hpp:144
LinearSolver::async_solvers::InitializeElement
Definition: ElementActions.hpp:102
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:59
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
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
db::AddSimpleTags
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1132
Convergence::Tags::IterationId
Identifies a step in an iterative algorithm.
Definition: Tags.hpp:74
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.
Initialization::merge_into_databox
auto merge_into_databox(db::DataBox< DbTagsList > &&box, Args &&... args) noexcept
Add tags that are not yet in the DataBox.
Definition: MergeIntoDataBox.hpp:133
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:161
LinearSolver::Tags::ResidualCompute
Compute the residual from the SourceTag and the db::add_tag_prefix<LinearSolver::Tags::OperatorAppl...
Definition: Tags.hpp:81
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:27
Convergence::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the iterative algorithm has converged,...
Definition: Tags.hpp:86
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
std::numeric_limits
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
observers::Actions::RegisterWithObservers
Register an observation ID with the observers.
Definition: RegisterWithObservers.hpp:39
LinearSolver::async_solvers::CompleteStep
Definition: ElementActions.hpp:215
string