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 #include <utility>
9 
10 #include "DataStructures/DataBox/DataBox.hpp"
11 #include "DataStructures/DataBox/PrefixHelpers.hpp"
12 #include "NumericalAlgorithms/Convergence/Tags.hpp"
13 #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
14 #include "Parallel/GlobalCache.hpp"
15 #include "Parallel/Invoke.hpp"
16 #include "Parallel/Reduction.hpp"
17 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitorActions.hpp"
18 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/Tags/InboxTags.hpp"
19 #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
20 #include "Utilities/Functional.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/Requires.hpp"
23 #include "Utilities/TMPL.hpp"
24 
25 /// \cond
26 namespace Convergence {
27 struct HasConverged;
28 } // namespace Convergence
29 namespace tuples {
30 template <typename...>
31 class TaggedTuple;
32 } // namespace tuples
33 namespace LinearSolver::cg::detail {
34 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
35 struct ResidualMonitor;
36 template <typename FieldsTag, typename OptionsGroup, typename Label>
37 struct UpdateOperand;
38 } // namespace LinearSolver::cg::detail
39 /// \endcond
40 
41 namespace LinearSolver::cg::detail {
42 
43 template <typename FieldsTag, typename OptionsGroup, typename Label,
44  typename SourceTag>
45 struct PrepareSolve {
46  private:
47  using fields_tag = FieldsTag;
48  using source_tag = SourceTag;
49  using operator_applied_to_fields_tag =
50  db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, fields_tag>;
51  using operand_tag =
52  db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
53  using residual_tag =
54  db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
55 
56  public:
57  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
58  typename ArrayIndex, typename ActionList,
59  typename ParallelComponent>
60  static std::tuple<db::DataBox<DbTagsList>&&> apply(
61  db::DataBox<DbTagsList>& box,
62  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
63  Parallel::GlobalCache<Metavariables>& cache,
64  const ArrayIndex& array_index, const ActionList /*meta*/,
65  const ParallelComponent* const /*meta*/) noexcept {
66  db::mutate<Convergence::Tags::IterationId<OptionsGroup>, operand_tag,
67  residual_tag>(
68  make_not_null(&box),
69  [](const gsl::not_null<size_t*> iteration_id, const auto operand,
70  const auto residual, const auto& source,
71  const auto& operator_applied_to_fields) noexcept {
72  *iteration_id = 0;
73  *operand = source - operator_applied_to_fields;
74  *residual = *operand;
75  },
76  get<source_tag>(box), get<operator_applied_to_fields_tag>(box));
77 
78  // Perform global reduction to compute initial residual magnitude square for
79  // residual monitor
80  const auto& residual = get<residual_tag>(box);
81  Parallel::contribute_to_reduction<
82  InitializeResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
83  Parallel::ReductionData<
84  Parallel::ReductionDatum<double, funcl::Plus<>>>{
85  inner_product(residual, residual)},
86  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
87  Parallel::get_parallel_component<
88  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
89 
90  return {std::move(box)};
91  }
92 };
93 
94 template <typename FieldsTag, typename OptionsGroup, typename Label>
95 struct InitializeHasConverged {
96  using inbox_tags = tmpl::list<Tags::InitialHasConverged<OptionsGroup>>;
97 
98  template <typename DbTags, typename... InboxTags, typename Metavariables,
99  typename ArrayIndex>
100  static bool is_ready(const db::DataBox<DbTags>& box,
101  const tuples::TaggedTuple<InboxTags...>& inboxes,
102  const Parallel::GlobalCache<Metavariables>& /*cache*/,
103  const ArrayIndex& /*array_index*/) noexcept {
104  const auto& inbox = get<Tags::InitialHasConverged<OptionsGroup>>(inboxes);
105  return inbox.find(db::get<Convergence::Tags::IterationId<OptionsGroup>>(
106  box)) != inbox.end();
107  }
108 
109  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
110  typename ArrayIndex, typename ActionList,
111  typename ParallelComponent>
112  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
113  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
114  const Parallel::GlobalCache<Metavariables>& /*cache*/,
115  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
116  const ParallelComponent* const /*meta*/) noexcept {
117  auto has_converged = std::move(
118  tuples::get<Tags::InitialHasConverged<OptionsGroup>>(inboxes)
119  .extract(db::get<Convergence::Tags::IterationId<OptionsGroup>>(box))
120  .mapped());
121 
122  db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
123  make_not_null(&box),
124  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
125  local_has_converged) noexcept {
126  *local_has_converged = std::move(has_converged);
127  });
128 
129  // Skip steps entirely if the solve has already converged
130  constexpr size_t step_end_index =
131  tmpl::index_of<ActionList,
132  UpdateOperand<FieldsTag, OptionsGroup, Label>>::value;
133  constexpr size_t this_action_index =
134  tmpl::index_of<ActionList, InitializeHasConverged>::value;
135  return {std::move(box), false,
136  has_converged ? (step_end_index + 1) : (this_action_index + 1)};
137  }
138 };
139 
140 template <typename FieldsTag, typename OptionsGroup, typename Label>
141 struct PerformStep {
142  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
143  typename ArrayIndex, typename ActionList,
144  typename ParallelComponent>
145  static std::tuple<db::DataBox<DbTagsList>&&> apply(
146  db::DataBox<DbTagsList>& box,
147  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
148  Parallel::GlobalCache<Metavariables>& cache,
149  const ArrayIndex& array_index, const ActionList /*meta*/,
150  const ParallelComponent* const /*meta*/) noexcept {
151  using fields_tag = FieldsTag;
152  using operand_tag =
153  db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
154  using operator_tag =
155  db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>;
156 
157  // At this point Ap must have been computed in a previous action
158  // We compute the inner product <p,p> w.r.t A. This requires a global
159  // reduction.
160  const double local_conj_grad_inner_product =
161  inner_product(get<operand_tag>(box), get<operator_tag>(box));
162 
163  Parallel::contribute_to_reduction<
164  ComputeAlpha<FieldsTag, OptionsGroup, ParallelComponent>>(
165  Parallel::ReductionData<
166  Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
167  Parallel::ReductionDatum<double, funcl::Plus<>>>{
168  get<Convergence::Tags::IterationId<OptionsGroup>>(box),
169  local_conj_grad_inner_product},
170  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
171  Parallel::get_parallel_component<
172  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
173 
174  return {std::move(box)};
175  }
176 };
177 
178 template <typename FieldsTag, typename OptionsGroup, typename Label>
179 struct UpdateFieldValues {
180  private:
181  using fields_tag = FieldsTag;
182  using operand_tag =
183  db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
184  using operator_tag =
185  db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>;
186  using residual_tag =
187  db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
188 
189  public:
190  using inbox_tags = tmpl::list<Tags::Alpha<OptionsGroup>>;
191 
192  template <typename DbTags, typename... InboxTags, typename Metavariables,
193  typename ArrayIndex>
194  static bool is_ready(const db::DataBox<DbTags>& box,
195  const tuples::TaggedTuple<InboxTags...>& inboxes,
196  const Parallel::GlobalCache<Metavariables>& /*cache*/,
197  const ArrayIndex& /*array_index*/) noexcept {
198  const auto& inbox = get<Tags::Alpha<OptionsGroup>>(inboxes);
199  return inbox.find(db::get<Convergence::Tags::IterationId<OptionsGroup>>(
200  box)) != inbox.end();
201  }
202 
203  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
204  typename ArrayIndex, typename ActionList,
205  typename ParallelComponent>
206  static std::tuple<db::DataBox<DbTagsList>&&> apply(
207  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
208  Parallel::GlobalCache<Metavariables>& cache,
209  const ArrayIndex& array_index, const ActionList /*meta*/,
210  const ParallelComponent* const /*meta*/) noexcept {
211  const double alpha = std::move(
212  tuples::get<Tags::Alpha<OptionsGroup>>(inboxes)
213  .extract(db::get<Convergence::Tags::IterationId<OptionsGroup>>(box))
214  .mapped());
215 
216  db::mutate<residual_tag, fields_tag>(
217  make_not_null(&box),
218  [alpha](const auto residual, const auto fields, const auto& operand,
219  const auto& operator_applied_to_operand) noexcept {
220  *fields += alpha * operand;
221  *residual -= alpha * operator_applied_to_operand;
222  },
223  get<operand_tag>(box), get<operator_tag>(box));
224 
225  // Compute new residual norm in a second global reduction
226  const auto& residual = get<residual_tag>(box);
227  const double local_residual_magnitude_square =
228  inner_product(residual, residual);
229 
230  Parallel::contribute_to_reduction<
231  UpdateResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
232  Parallel::ReductionData<
233  Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
234  Parallel::ReductionDatum<double, funcl::Plus<>>>{
235  get<Convergence::Tags::IterationId<OptionsGroup>>(box),
236  local_residual_magnitude_square},
237  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
238  Parallel::get_parallel_component<
239  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
240 
241  return {std::move(box)};
242  }
243 };
244 
245 template <typename FieldsTag, typename OptionsGroup, typename Label>
246 struct UpdateOperand {
247  private:
248  using fields_tag = FieldsTag;
249  using operand_tag =
250  db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
251  using residual_tag =
252  db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
253 
254  public:
255  using inbox_tags =
256  tmpl::list<Tags::ResidualRatioAndHasConverged<OptionsGroup>>;
257 
258  template <typename DbTags, typename... InboxTags, typename Metavariables,
259  typename ArrayIndex>
260  static bool is_ready(const db::DataBox<DbTags>& box,
261  const tuples::TaggedTuple<InboxTags...>& inboxes,
262  const Parallel::GlobalCache<Metavariables>& /*cache*/,
263  const ArrayIndex& /*array_index*/) noexcept {
264  const auto& inbox =
265  get<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(inboxes);
266  return inbox.find(db::get<Convergence::Tags::IterationId<OptionsGroup>>(
267  box)) != inbox.end();
268  }
269 
270  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
271  typename ArrayIndex, typename ActionList,
272  typename ParallelComponent>
273  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
274  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
275  const Parallel::GlobalCache<Metavariables>& /*cache*/,
276  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
277  const ParallelComponent* const /*meta*/) noexcept {
278  auto received_data = std::move(
279  tuples::get<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(inboxes)
280  .extract(db::get<Convergence::Tags::IterationId<OptionsGroup>>(box))
281  .mapped());
282  const double res_ratio = get<0>(received_data);
283  auto& has_converged = get<1>(received_data);
284 
285  db::mutate<operand_tag, Convergence::Tags::IterationId<OptionsGroup>,
286  Convergence::Tags::HasConverged<OptionsGroup>>(
287  make_not_null(&box),
288  [res_ratio, &has_converged](
289  const auto operand, const gsl::not_null<size_t*> iteration_id,
290  const gsl::not_null<Convergence::HasConverged*> local_has_converged,
291  const auto& residual) noexcept {
292  *operand = residual + res_ratio * *operand;
293  ++(*iteration_id);
294  *local_has_converged = std::move(has_converged);
295  },
296  get<residual_tag>(box));
297 
298  // Repeat steps until the solve has converged
299  constexpr size_t this_action_index =
300  tmpl::index_of<ActionList, UpdateOperand>::value;
301  constexpr size_t prepare_step_index =
302  tmpl::index_of<ActionList, InitializeHasConverged<
303  FieldsTag, OptionsGroup, Label>>::value +
304  1;
305  return {std::move(box), false,
306  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
307  ? (this_action_index + 1)
308  : prepare_step_index};
309  }
310 };
311 
312 } // namespace LinearSolver::cg::detail
std::apply
T apply(T... args)
utility
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
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: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:535
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:786
Convergence
Items related to checking the convergence of numerical algorithms.
Definition: Convergence.hpp:10
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
ActionTesting::is_ready
bool is_ready(MockRuntimeSystem< Metavariables > &runner, const typename Component::array_index &array_index) noexcept
Runs the is_ready function and returns the result for the next action in the current phase on the arr...
Definition: MockRuntimeSystemFreeFunctions.hpp:286
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
Convergence::Tags::IterationId
Identifies a step in an iterative algorithm.
Definition: Tags.hpp:74
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:380
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
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: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
© Copyright 2017 - 2021 SXS Collaboration, MIT License