ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <tuple>
8 
10 #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 #include "Parallel/GlobalCache.hpp"
14 #include "Parallel/Info.hpp"
15 #include "Parallel/Invoke.hpp"
16 #include "Parallel/Reduction.hpp"
17 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitorActions.hpp"
19 #include "Utilities/Functional.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/Requires.hpp"
22 
23 /// \cond
24 namespace Convergence {
25 struct HasConverged;
26 } // namespace Convergence
27 namespace tuples {
28 template <typename...>
29 class TaggedTuple;
30 } // namespace tuples
31 namespace LinearSolver::cg::detail {
32 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
33 struct ResidualMonitor;
34 } // namespace LinearSolver::cg::detail
35 /// \endcond
36 
37 namespace LinearSolver::cg::detail {
38 
39 template <typename FieldsTag, typename OptionsGroup>
40 struct PrepareSolve {
41  private:
42  using fields_tag = FieldsTag;
44  using operator_applied_to_fields_tag =
46  using operand_tag =
48  using residual_tag =
50 
51  public:
52  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
53  typename ArrayIndex, typename ActionList,
54  typename ParallelComponent>
57  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
59  const ArrayIndex& array_index, const ActionList /*meta*/,
60  const ParallelComponent* const /*meta*/) noexcept {
61  db::mutate<LinearSolver::Tags::IterationId<OptionsGroup>, operand_tag,
62  residual_tag>(
63  make_not_null(&box),
64  [](const gsl::not_null<size_t*> iteration_id, const auto operand,
65  const auto residual, const auto& source,
66  const auto& operator_applied_to_fields) noexcept {
67  *iteration_id = 0;
68  *operand = source - operator_applied_to_fields;
69  *residual = *operand;
70  },
71  get<source_tag>(box), get<operator_applied_to_fields_tag>(box));
72 
73  // Perform global reduction to compute initial residual magnitude square for
74  // residual monitor
75  const auto& residual = get<residual_tag>(box);
77  InitializeResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
78  Parallel::ReductionData<
80  inner_product(residual, residual)},
81  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
83  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
84 
85  return {
86  std::move(box),
87  // Terminate algorithm for now. The `ResidualMonitor` will receive the
88  // reduction that is performed above and then broadcast to the following
89  // action, which is responsible for restarting the algorithm.
90  true};
91  }
92 };
93 
94 template <typename FieldsTag, typename OptionsGroup>
95 struct InitializeHasConverged {
96  template <
97  typename ParallelComponent, typename DbTagsList, typename Metavariables,
98  typename ArrayIndex, typename DataBox = db::DataBox<DbTagsList>,
99  Requires<db::tag_is_retrievable_v<
101  static void apply(db::DataBox<DbTagsList>& box,
103  const ArrayIndex& array_index,
104  const Convergence::HasConverged& has_converged) noexcept {
105  db::mutate<LinearSolver::Tags::HasConverged<OptionsGroup>>(
106  make_not_null(&box),
107  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
108  local_has_converged) noexcept {
109  *local_has_converged = has_converged;
110  });
111 
112  // Proceed with algorithm
113  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
114  .perform_algorithm(true);
115  }
116 };
117 
118 template <typename FieldsTag, typename OptionsGroup>
119 struct PrepareStep {
120  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
121  typename ArrayIndex, typename ActionList,
122  typename ParallelComponent>
125  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
126  const Parallel::GlobalCache<Metavariables>& /*cache*/,
127  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
128  const ParallelComponent* const /*meta*/) noexcept {
129  // Nothing to do before applying the linear operator to the operand
130  return {std::move(box)};
131  }
132 };
133 
134 template <typename FieldsTag, typename OptionsGroup>
135 struct PerformStep {
136  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
137  typename ArrayIndex, typename ActionList,
138  typename ParallelComponent>
141  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
143  const ArrayIndex& array_index,
144  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
145  const ActionList /*meta*/,
146  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
147  const ParallelComponent* const /*meta*/) noexcept {
148  using fields_tag = FieldsTag;
149  using operand_tag =
151  using operator_tag =
153 
154  // At this point Ap must have been computed in a previous action
155  // We compute the inner product <p,p> w.r.t A. This requires a global
156  // reduction.
157  const double local_conj_grad_inner_product =
158  inner_product(get<operand_tag>(box), get<operator_tag>(box));
159 
161  ComputeAlpha<FieldsTag, OptionsGroup, ParallelComponent>>(
162  Parallel::ReductionData<
164  local_conj_grad_inner_product},
165  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
167  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
168 
169  // Terminate algorithm for now. The `ResidualMonitor` will receive the
170  // reduction that is performed above and then broadcast to the following
171  // action, which is responsible for restarting the algorithm.
172  return {std::move(box), true};
173  }
174 };
175 
176 template <typename FieldsTag, typename OptionsGroup>
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  template <typename ParallelComponent, typename DbTagsList,
189  typename Metavariables, typename ArrayIndex,
190  typename DataBox = db::DataBox<DbTagsList>,
192  db::tag_is_retrievable_v<fields_tag, DataBox> and
193  db::tag_is_retrievable_v<operand_tag, DataBox> and
194  db::tag_is_retrievable_v<operator_tag, DataBox>> = nullptr>
195  static auto apply(db::DataBox<DbTagsList>& box,
197  const ArrayIndex& array_index,
198  const double alpha) noexcept {
199  // Received global reduction result, proceed with conjugate gradient.
200  db::mutate<residual_tag, fields_tag>(
201  make_not_null(&box),
202  [alpha](const auto residual, const auto fields, const auto& operand,
203  const auto& operator_applied_to_operand) noexcept {
204  *fields += alpha * operand;
205  *residual -= alpha * operator_applied_to_operand;
206  },
207  get<operand_tag>(box), get<operator_tag>(box));
208 
209  // Compute new residual norm in a second global reduction
210  const auto& residual = get<residual_tag>(box);
211  const double local_residual_magnitude_square =
212  inner_product(residual, residual);
213 
215  UpdateResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
216  Parallel::ReductionData<
218  local_residual_magnitude_square},
219  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
221  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
222  }
223 };
224 
225 template <typename FieldsTag, typename OptionsGroup>
226 struct UpdateOperand {
227  private:
228  using fields_tag = FieldsTag;
229  using operand_tag =
231  using residual_tag =
233 
234  public:
235  template <
236  typename ParallelComponent, typename DbTagsList, typename Metavariables,
237  typename ArrayIndex, typename DataBox = db::DataBox<DbTagsList>,
239  db::tag_is_retrievable_v<
241  db::tag_is_retrievable_v<
243  db::tag_is_retrievable_v<residual_tag, DataBox>> = nullptr>
244  static auto apply(db::DataBox<DbTagsList>& box,
246  const ArrayIndex& array_index, const double res_ratio,
247  const Convergence::HasConverged& has_converged) noexcept {
248  db::mutate<operand_tag, LinearSolver::Tags::IterationId<OptionsGroup>,
250  make_not_null(&box),
251  [res_ratio, &has_converged](
252  const auto operand, const gsl::not_null<size_t*> iteration_id,
253  const gsl::not_null<Convergence::HasConverged*> local_has_converged,
254  const auto& residual) noexcept {
255  *operand = residual + res_ratio * *operand;
256  ++(*iteration_id);
257  *local_has_converged = has_converged;
258  },
259  get<residual_tag>(box));
260 
261  // Proceed with algorithm
262  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
263  .perform_algorithm(true);
264  }
265 };
266 
267 } // namespace LinearSolver::cg::detail
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
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:223
LinearSolver::Tags::IterationId
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:84
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:52
tuple
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:59
Convergence
Items related to checking the convergence of numerical algorithms.
Definition: Convergence.hpp:10
Info.hpp
Tags.hpp
Parallel::contribute_to_reduction
void contribute_to_reduction(ReductionData< Ts... > reduction_data, const SenderProxy &sender_component, const TargetProxy &target_component) noexcept
Perform a reduction from the sender_component (typically your own parallel component) to the target_c...
Definition: Reduction.hpp:245
db::apply
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1424
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
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
LinearSolver::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the linear solver has converged,...
Definition: Tags.hpp:234
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
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
Prefixes.hpp
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183