ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <tuple>
7 
10 #include "NumericalAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitorActions.hpp"
15 #include "Parallel/Info.hpp"
16 #include "Parallel/Invoke.hpp"
17 #include "Parallel/Reduction.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>
30 struct ResidualMonitor;
31 } // namespace cg_detail
32 } // namespace LinearSolver
33 /// \endcond
34 
35 namespace LinearSolver {
36 namespace cg_detail {
37 
38 struct InitializeHasConverged {
39  template <typename... DbTags, typename... InboxTags, typename Metavariables,
40  typename ArrayIndex, typename ActionList,
41  typename ParallelComponent,
42  Requires<sizeof...(DbTags) != 0> = nullptr>
43  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
44  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
46  const ArrayIndex& /*array_index*/,
47  const ActionList /*meta*/,
48  const ParallelComponent* const /*meta*/,
50  has_converged) noexcept {
51  db::mutate<LinearSolver::Tags::HasConverged>(
52  make_not_null(&box), [&has_converged](
55  local_has_converged) noexcept {
56  *local_has_converged = has_converged;
57  });
58  }
59 };
60 
61 struct PerformStep {
62  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
63  typename ArrayIndex, typename ActionList,
64  typename ParallelComponent>
65  static auto apply(db::DataBox<DbTagsList>& box,
66  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
68  const ArrayIndex& array_index, const ActionList /*meta*/,
69  const ParallelComponent* const /*meta*/) noexcept {
70  using fields_tag = typename Metavariables::system::fields_tag;
71  using operand_tag =
73  using operator_tag =
75 
76  // At this point Ap must have been computed in a previous action
77  // We compute the inner product <p,p> w.r.t A. This requires a global
78  // reduction.
79  const double local_conj_grad_inner_product =
80  inner_product(get<operand_tag>(box), get<operator_tag>(box));
81 
82  Parallel::contribute_to_reduction<ComputeAlpha<ParallelComponent>>(
83  Parallel::ReductionData<
85  local_conj_grad_inner_product},
86  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
87  Parallel::get_parallel_component<ResidualMonitor<Metavariables>>(
88  cache));
89 
90  // Terminate algorithm for now. The reduction will be broadcast to the
91  // next action which is responsible for restarting the algorithm.
92  return std::tuple<db::DataBox<DbTagsList>&&, bool>(std::move(box), true);
93  }
94 };
95 
96 struct UpdateFieldValues {
97  template <typename... DbTags, typename... InboxTags, typename Metavariables,
98  typename ArrayIndex, typename ActionList,
99  typename ParallelComponent,
100  Requires<sizeof...(DbTags) != 0> = nullptr>
101  static auto apply(db::DataBox<tmpl::list<DbTags...>>& box,
102  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
104  const ArrayIndex& array_index, const ActionList /*meta*/,
105  const ParallelComponent* const /*meta*/,
106  const double alpha) noexcept {
107  using fields_tag = typename Metavariables::system::fields_tag;
108  using operand_tag =
110  using operator_tag =
112  using residual_tag =
114 
115  // Received global reduction result, proceed with conjugate gradient.
116  db::mutate<residual_tag, fields_tag>(
117  make_not_null(&box),
121  const db::item_type<operator_tag>& Ap) noexcept {
122  *x += alpha * p;
123  *r -= alpha * Ap;
124  },
125  get<operand_tag>(box), get<operator_tag>(box));
126 
127  // Compute new residual norm in a second global reduction
128  const auto& r = get<residual_tag>(box);
129  const double local_residual_magnitude_square = inner_product(r, r);
130 
131  Parallel::contribute_to_reduction<UpdateResidual<ParallelComponent>>(
132  Parallel::ReductionData<
134  local_residual_magnitude_square},
135  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
136  Parallel::get_parallel_component<ResidualMonitor<Metavariables>>(
137  cache));
138  }
139 };
140 
141 struct UpdateOperand {
142  template <typename... DbTags, typename... InboxTags, typename Metavariables,
143  typename ArrayIndex, typename ActionList,
144  typename ParallelComponent,
145  Requires<sizeof...(DbTags) != 0> = nullptr>
146  static auto apply(db::DataBox<tmpl::list<DbTags...>>& box,
147  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
149  const ArrayIndex& array_index, const ActionList /*meta*/,
150  const ParallelComponent* const /*meta*/,
151  const double res_ratio,
153  has_converged) noexcept {
154  using fields_tag = typename Metavariables::system::fields_tag;
155  using operand_tag =
157  using residual_tag =
159 
160  // Prepare conjugate gradient for next iteration
164  make_not_null(&box),
165  [
166  res_ratio, &has_converged
169  local_has_converged,
170  const gsl::not_null<IterationId*> iteration_id,
171  const gsl::not_null<IterationId*> next_iteration_id,
172  const db::item_type<residual_tag>& r) noexcept {
173  *p = r + res_ratio * *p;
174  *local_has_converged = has_converged;
175  iteration_id->step_number++;
176  next_iteration_id->step_number = iteration_id->step_number + 1;
177  },
178  get<residual_tag>(box));
179 
180  // Proceed with algorithm
181  // We use `ckLocal()` here since this is essentially retrieving "self",
182  // which is guaranteed to be on the local processor. This ensures the calls
183  // are evaluated in order.
184  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
185  .ckLocal()
186  ->set_terminate(false);
187  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
188  .perform_algorithm();
189  }
190 };
191 
192 } // namespace cg_detail
193 } // 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: DataBoxTag.hpp:533
Definition: Variables.hpp:46
Defines functions for interfacing with the parallelization framework.
Defines DataBox tags for the linear solver.
void mutate(const gsl::not_null< DataBox< TagList > *> box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:1099
Functionality for solving linear systems of equations.
Definition: TerminateIfConverged.hpp:22
Definition: TaggedTuple.hpp:25
Define prefixes for DataBox tags.
Holds a LinearSolver::HasConverged flag that signals the linear solver has converged, along with the reason for convergence.
Definition: Tags.hpp:204
Defines the type alias Requires.
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:66
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
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:76
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:163
Prefix indicating the value a quantity will take on the next iteration of the algorithm.
Definition: Prefixes.hpp:98
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:410
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:863
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: SolvePoissonProblem.hpp:38
Defines class template ConstGlobalCache.
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: ConservativeFromPrimitive.hpp:12
Defines class IterationId.