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 <tuple>
9 #include <utility>
10 
12 #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 #include "NumericalAlgorithms/Convergence/Tags.hpp"
15 #include "Parallel/AlgorithmMetafunctions.hpp"
16 #include "Parallel/GlobalCache.hpp"
17 #include "Parallel/Invoke.hpp"
18 #include "Parallel/Reduction.hpp"
19 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitorActions.hpp"
20 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/Tags/InboxTags.hpp"
22 #include "Utilities/Functional.hpp"
23 #include "Utilities/Gsl.hpp"
24 #include "Utilities/Requires.hpp"
25 #include "Utilities/TMPL.hpp"
26 
27 /// \cond
28 namespace Convergence {
29 struct HasConverged;
30 } // namespace Convergence
31 namespace tuples {
32 template <typename...>
33 class TaggedTuple;
34 } // namespace tuples
35 namespace LinearSolver::cg::detail {
36 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
37 struct ResidualMonitor;
38 template <typename FieldsTag, typename OptionsGroup, typename Label>
39 struct UpdateOperand;
40 } // namespace LinearSolver::cg::detail
41 /// \endcond
42 
43 namespace LinearSolver::cg::detail {
44 
45 template <typename FieldsTag, typename OptionsGroup, typename Label,
46  typename SourceTag>
47 struct PrepareSolve {
48  private:
49  using fields_tag = FieldsTag;
50  using source_tag = SourceTag;
51  using operator_applied_to_fields_tag =
53  using operand_tag =
55  using residual_tag =
57 
58  public:
59  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
60  typename ArrayIndex, typename ActionList,
61  typename ParallelComponent>
63  db::DataBox<DbTagsList>& box,
64  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
66  const ArrayIndex& array_index, const ActionList /*meta*/,
67  const ParallelComponent* const /*meta*/) noexcept {
68  db::mutate<Convergence::Tags::IterationId<OptionsGroup>, operand_tag,
69  residual_tag>(
70  make_not_null(&box),
71  [](const gsl::not_null<size_t*> iteration_id, const auto operand,
72  const auto residual, const auto& source,
73  const auto& operator_applied_to_fields) noexcept {
74  *iteration_id = 0;
75  *operand = source - operator_applied_to_fields;
76  *residual = *operand;
77  },
78  get<source_tag>(box), get<operator_applied_to_fields_tag>(box));
79 
80  // Perform global reduction to compute initial residual magnitude square for
81  // residual monitor
82  const auto& residual = get<residual_tag>(box);
84  InitializeResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
85  Parallel::ReductionData<
87  inner_product(residual, residual)},
88  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
90  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
91 
92  return {std::move(box)};
93  }
94 };
95 
96 template <typename FieldsTag, typename OptionsGroup, typename Label>
97 struct InitializeHasConverged {
98  using inbox_tags = tmpl::list<Tags::InitialHasConverged<OptionsGroup>>;
99 
100  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
101  typename ArrayIndex, typename ActionList,
102  typename ParallelComponent>
104  size_t>
105  apply(db::DataBox<DbTagsList>& box,
107  const Parallel::GlobalCache<Metavariables>& /*cache*/,
108  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
109  const ParallelComponent* const /*meta*/) noexcept {
110  auto& inbox = get<Tags::InitialHasConverged<OptionsGroup>>(inboxes);
111  const auto& iteration_id =
112  db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
113  if (inbox.find(iteration_id) == inbox.end()) {
114  return {std::move(box), Parallel::AlgorithmExecution::Retry,
116  }
117 
118  auto has_converged = std::move(inbox.extract(iteration_id).mapped());
119 
120  db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
121  make_not_null(&box),
122  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
123  local_has_converged) noexcept {
124  *local_has_converged = std::move(has_converged);
125  });
126 
127  // Skip steps entirely if the solve has already converged
128  constexpr size_t step_end_index =
129  tmpl::index_of<ActionList,
130  UpdateOperand<FieldsTag, OptionsGroup, Label>>::value;
131  constexpr size_t this_action_index =
132  tmpl::index_of<ActionList, InitializeHasConverged>::value;
133  return {std::move(box), Parallel::AlgorithmExecution::Continue,
134  has_converged ? (step_end_index + 1) : (this_action_index + 1)};
135  }
136 };
137 
138 template <typename FieldsTag, typename OptionsGroup, typename Label>
139 struct PerformStep {
140  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
141  typename ArrayIndex, typename ActionList,
142  typename ParallelComponent>
144  db::DataBox<DbTagsList>& box,
145  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
147  const ArrayIndex& array_index, const ActionList /*meta*/,
148  const ParallelComponent* const /*meta*/) noexcept {
149  using fields_tag = FieldsTag;
150  using operand_tag =
152  using operator_tag =
154 
155  // At this point Ap must have been computed in a previous action
156  // We compute the inner product <p,p> w.r.t A. This requires a global
157  // reduction.
158  const double local_conj_grad_inner_product =
159  inner_product(get<operand_tag>(box), get<operator_tag>(box));
160 
162  ComputeAlpha<FieldsTag, OptionsGroup, ParallelComponent>>(
163  Parallel::ReductionData<
166  get<Convergence::Tags::IterationId<OptionsGroup>>(box),
167  local_conj_grad_inner_product},
168  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
170  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
171 
172  return {std::move(box)};
173  }
174 };
175 
176 template <typename FieldsTag, typename OptionsGroup, typename Label>
177 struct UpdateFieldValues {
178  private:
179  using fields_tag = FieldsTag;
180  using operand_tag =
182  using operator_tag =
184  using residual_tag =
186 
187  public:
188  using inbox_tags = tmpl::list<Tags::Alpha<OptionsGroup>>;
189 
190  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
191  typename ArrayIndex, typename ActionList,
192  typename ParallelComponent>
194  apply(db::DataBox<DbTagsList>& box,
197  const ArrayIndex& array_index, const ActionList /*meta*/,
198  const ParallelComponent* const /*meta*/) noexcept {
199  auto& inbox = get<Tags::Alpha<OptionsGroup>>(inboxes);
200  const auto& iteration_id =
201  db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
202  if (inbox.find(iteration_id) == inbox.end()) {
203  return {std::move(box), Parallel::AlgorithmExecution::Retry};
204  }
205 
206  const double alpha = std::move(inbox.extract(iteration_id).mapped());
207 
208  db::mutate<residual_tag, fields_tag>(
209  make_not_null(&box),
210  [alpha](const auto residual, const auto fields, const auto& operand,
211  const auto& operator_applied_to_operand) noexcept {
212  *fields += alpha * operand;
213  *residual -= alpha * operator_applied_to_operand;
214  },
215  get<operand_tag>(box), get<operator_tag>(box));
216 
217  // Compute new residual norm in a second global reduction
218  const auto& residual = get<residual_tag>(box);
219  const double local_residual_magnitude_square =
220  inner_product(residual, residual);
221 
223  UpdateResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
224  Parallel::ReductionData<
227  get<Convergence::Tags::IterationId<OptionsGroup>>(box),
228  local_residual_magnitude_square},
229  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
231  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
232 
233  return {std::move(box), Parallel::AlgorithmExecution::Continue};
234  }
235 };
236 
237 template <typename FieldsTag, typename OptionsGroup, typename Label>
238 struct UpdateOperand {
239  private:
240  using fields_tag = FieldsTag;
241  using operand_tag =
243  using residual_tag =
245 
246  public:
247  using inbox_tags =
248  tmpl::list<Tags::ResidualRatioAndHasConverged<OptionsGroup>>;
249 
250  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
251  typename ArrayIndex, typename ActionList,
252  typename ParallelComponent>
254  size_t>
255  apply(db::DataBox<DbTagsList>& box,
257  const Parallel::GlobalCache<Metavariables>& /*cache*/,
258  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
259  const ParallelComponent* const /*meta*/) noexcept {
260  auto& inbox =
261  get<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(inboxes);
262  const auto& iteration_id =
263  db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
264  if (inbox.find(iteration_id) == inbox.end()) {
265  return {std::move(box), Parallel::AlgorithmExecution::Retry,
267  }
268 
269  auto received_data = std::move(inbox.extract(iteration_id).mapped());
270  const double res_ratio = get<0>(received_data);
271  auto& has_converged = get<1>(received_data);
272 
273  db::mutate<operand_tag, Convergence::Tags::IterationId<OptionsGroup>,
275  make_not_null(&box),
276  [res_ratio, &has_converged](
277  const auto operand, const gsl::not_null<size_t*> local_iteration_id,
278  const gsl::not_null<Convergence::HasConverged*> local_has_converged,
279  const auto& residual) noexcept {
280  *operand = residual + res_ratio * *operand;
281  ++(*local_iteration_id);
282  *local_has_converged = std::move(has_converged);
283  },
284  get<residual_tag>(box));
285 
286  // Repeat steps until the solve has converged
287  constexpr size_t this_action_index =
288  tmpl::index_of<ActionList, UpdateOperand>::value;
289  constexpr size_t prepare_step_index =
290  tmpl::index_of<ActionList, InitializeHasConverged<
291  FieldsTag, OptionsGroup, Label>>::value +
292  1;
293  return {std::move(box), Parallel::AlgorithmExecution::Continue,
294  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
295  ? (this_action_index + 1)
296  : prepare_step_index};
297  }
298 };
299 
300 } // namespace LinearSolver::cg::detail
std::apply
T apply(T... args)
utility
Parallel::AlgorithmExecution::Retry
@ Retry
Temporarily stop executing iterable actions, but try the same action again after receiving data from ...
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:65
InnerProduct.hpp
GlobalCache.hpp
Parallel::get_parallel_component
auto get_parallel_component(GlobalCache< Metavariables > &cache) noexcept -> Parallel::proxy_from_parallel_component< GlobalCache_detail::get_component_if_mocked< typename Metavariables::component_list, ParallelComponentTag >> &
Access the Charm++ proxy associated with a ParallelComponent.
Definition: GlobalCache.hpp:535
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
Convergence
Items related to checking the convergence of numerical algorithms.
Definition: Convergence.hpp:10
Tags.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
Parallel::AlgorithmExecution::Continue
@ Continue
Leave the algorithm termination flag in its current state.
DataBox.hpp
cstddef
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
Parallel::AlgorithmExecution
AlgorithmExecution
The possible options for altering the current execution of the algorithm, used in the return type of ...
Definition: AlgorithmMetafunctions.hpp:31
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
limits
Gsl.hpp
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
std::numeric_limits::max
T max(T... args)
Parallel::contribute_to_reduction
void contribute_to_reduction(ReductionData< Ts... > reduction_data, const SenderProxy &sender_component, const TargetProxy &target_component, [[maybe_unused]] const gsl::not_null< SectionType * > section=&no_section()) noexcept
Perform a reduction from the sender_component (typically your own parallel component) to the target_c...
Definition: Reduction.hpp:272
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: ReadSpecPiecewisePolynomial.hpp:13