ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <tuple>
7 
9 #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 #include "Parallel/Info.hpp"
14 #include "Parallel/Invoke.hpp"
15 #include "Parallel/Reduction.hpp"
16 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitorActions.hpp"
18 #include "Utilities/Functional.hpp"
19 #include "Utilities/Gsl.hpp"
20 #include "Utilities/Requires.hpp"
21 
22 /// \cond
23 namespace tuples {
24 template <typename...>
25 class TaggedTuple;
26 } // namespace tuples
27 namespace LinearSolver {
28 namespace cg_detail {
29 template <typename Metavariables, typename FieldsTag>
30 struct ResidualMonitor;
31 } // namespace cg_detail
32 } // namespace LinearSolver
33 /// \endcond
34 
35 namespace LinearSolver {
36 namespace cg_detail {
37 
38 struct PrepareStep {
39  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
40  typename ArrayIndex, typename ActionList,
41  typename ParallelComponent>
44  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
46  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
47  const ParallelComponent* const /*meta*/) noexcept {
48  db::mutate<LinearSolver::Tags::IterationId>(
49  make_not_null(&box),
51  iteration_id,
53  LinearSolver::Tags::IterationId>>& next_iteration_id) noexcept {
54  *iteration_id = next_iteration_id;
55  },
56  get<::Tags::Next<LinearSolver::Tags::IterationId>>(box));
57  return {std::move(box)};
58  }
59 };
60 
61 template <typename FieldsTag>
62 struct PerformStep {
63  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
64  typename ArrayIndex, typename ActionList,
65  typename ParallelComponent>
68  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
70  const ArrayIndex& array_index,
71  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
72  const ActionList /*meta*/,
73  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
74  const ParallelComponent* const /*meta*/) noexcept {
75  using fields_tag = FieldsTag;
76  using operand_tag =
78  using operator_tag =
80 
81  ASSERT(get<LinearSolver::Tags::IterationId>(box) !=
83  "Linear solve iteration ID is at initial state. Did you forget to "
84  "invoke 'PrepareStep'?");
85 
86  // At this point Ap must have been computed in a previous action
87  // We compute the inner product <p,p> w.r.t A. This requires a global
88  // reduction.
89  const double local_conj_grad_inner_product =
90  inner_product(get<operand_tag>(box), get<operator_tag>(box));
91 
93  ComputeAlpha<FieldsTag, ParallelComponent>>(
94  Parallel::ReductionData<
96  local_conj_grad_inner_product},
97  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
99  ResidualMonitor<Metavariables, FieldsTag>>(cache));
100 
101  // Terminate algorithm for now. The `ResidualMonitor` will receive the
102  // reduction that is performed above and then broadcast to the following
103  // action, which is responsible for restarting the algorithm.
104  return {std::move(box), true};
105  }
106 };
107 
108 template <typename FieldsTag>
109 struct UpdateFieldValues {
110  private:
111  using fields_tag = FieldsTag;
112  using operand_tag =
114  using operator_tag =
116  using residual_tag =
118 
119  public:
120  template <typename ParallelComponent, typename DbTagsList,
121  typename Metavariables, typename ArrayIndex,
122  typename DataBox = db::DataBox<DbTagsList>,
124  db::tag_is_retrievable_v<fields_tag, DataBox> and
125  db::tag_is_retrievable_v<operand_tag, DataBox> and
126  db::tag_is_retrievable_v<operator_tag, DataBox>> = nullptr>
127  static auto apply(db::DataBox<DbTagsList>& box,
129  const ArrayIndex& array_index,
130  const double alpha) noexcept {
131  // Received global reduction result, proceed with conjugate gradient.
132  db::mutate<residual_tag, fields_tag>(
133  make_not_null(&box),
137  const db::const_item_type<operator_tag>& Ap) noexcept {
138  *x += alpha * p;
139  *r -= alpha * Ap;
140  },
141  get<operand_tag>(box), get<operator_tag>(box));
142 
143  // Compute new residual norm in a second global reduction
144  const auto& r = get<residual_tag>(box);
145  const double local_residual_magnitude_square = inner_product(r, r);
146 
148  UpdateResidual<FieldsTag, ParallelComponent>>(
149  Parallel::ReductionData<
151  local_residual_magnitude_square},
152  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
154  ResidualMonitor<Metavariables, FieldsTag>>(cache));
155  }
156 };
157 
158 template <typename FieldsTag>
159 struct UpdateOperand {
160  private:
161  using fields_tag = FieldsTag;
162  using operand_tag =
164  using residual_tag =
166 
167  public:
168  template <typename ParallelComponent, typename DbTagsList,
169  typename Metavariables, typename ArrayIndex,
170  typename DataBox = db::DataBox<DbTagsList>,
173  DataBox> and
174  db::tag_is_retrievable_v<residual_tag, DataBox>> = nullptr>
175  static auto apply(db::DataBox<DbTagsList>& box,
177  const ArrayIndex& array_index, const double res_ratio,
179  has_converged) noexcept {
180  db::mutate<operand_tag, LinearSolver::Tags::HasConverged>(
181  make_not_null(&box),
182  [
183  res_ratio, &has_converged
186  local_has_converged,
187  const db::const_item_type<residual_tag>& r) noexcept {
188  *p = r + res_ratio * *p;
189  *local_has_converged = has_converged;
190  },
191  get<residual_tag>(box));
192 
193  // Proceed with algorithm
194  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
195  .perform_algorithm(true);
196  }
197 };
198 
199 } // namespace cg_detail
200 } // namespace LinearSolver
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Wrap Tag in Prefix<_, Args...>, also wrapping variables tags if Tag is a Tags::Variables.
Definition: PrefixHelpers.hpp:145
Defines functions for interfacing with the parallelization framework.
Functionality for solving linear systems of equations.
Definition: Gmres.cpp:16
Definition: TaggedTuple.hpp:26
Define prefixes for DataBox tags.
Holds a Convergence::HasConverged flag that signals the linear solver has converged, along with the reason for convergence.
Definition: Tags.hpp:209
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:243
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
Defines the type alias Requires.
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:79
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
constexpr bool tag_is_retrievable_v
Equal to true if Tag can be retrieved from a DataBox of type DataBoxType.
Definition: DataBox.hpp:77
Definition: InterpolationTargetWedgeSectionTorus.hpp:25
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
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:136
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:1444
Defines an inner product for the linear solver.
auto get_parallel_component(ConstGlobalCache< Metavariables > &cache) noexcept -> Parallel::proxy_from_parallel_component< ConstGlobalCache_detail::get_component_if_mocked< typename Metavariables::component_list, ParallelComponentTag >> &
Access the Charm++ proxy associated with a ParallelComponent.
Definition: ConstGlobalCache.hpp:223
Prefix indicating the value a quantity will take on the next iteration of the algorithm.
Definition: Prefixes.hpp:152
Defines functions and classes from the GSL.
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
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
Definition: Test_ActionTesting.cpp:365
Defines class template ConstGlobalCache.
Defines DataBox tags for the linear solver.
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182