ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <tuple>
7 
12 #include "Parallel/Info.hpp"
13 #include "Parallel/Invoke.hpp"
14 #include "Parallel/Reduction.hpp"
15 #include "ParallelAlgorithms/LinearSolver/Gmres/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 gmres_detail {
29 template <typename Metavariables, typename FieldsTag>
30 struct ResidualMonitor;
31 } // namespace gmres_detail
32 } // namespace LinearSolver
33 /// \endcond
34 
35 namespace LinearSolver {
36 namespace gmres_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  using orthogonalization_iteration_id_tag =
51 
52  db::mutate<LinearSolver::Tags::IterationId,
53  orthogonalization_iteration_id_tag>(
54  make_not_null(&box),
56  iteration_id,
57  const gsl::not_null<
59  orthogonalization_iteration_id,
61  LinearSolver::Tags::IterationId>>& next_iteration_id) noexcept {
62  *iteration_id = next_iteration_id;
63  *orthogonalization_iteration_id = 0;
64  },
65  get<::Tags::Next<LinearSolver::Tags::IterationId>>(box));
66  return {std::move(box)};
67  }
68 };
69 
70 template <typename FieldsTag>
71 struct PerformStep {
72  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
73  typename ArrayIndex, typename ActionList,
74  typename ParallelComponent>
77  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
79  const ArrayIndex& array_index,
80  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
81  const ActionList /*meta*/,
82  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
83  const ParallelComponent* const /*meta*/) noexcept {
84  using fields_tag = FieldsTag;
85  using operand_tag =
87  using operator_tag =
89  using basis_history_tag =
91 
92  ASSERT(get<LinearSolver::Tags::IterationId>(box) !=
94  "Linear solve iteration ID is at initial state. Did you forget to "
95  "invoke 'PrepareStep'?");
96 
97  db::mutate<operand_tag>(
98  make_not_null(&box),
100  const db::const_item_type<operator_tag>& operator_action) noexcept {
101  *operand = db::item_type<operand_tag>(operator_action);
102  },
103  get<operator_tag>(box));
104 
106  StoreOrthogonalization<FieldsTag, ParallelComponent>>(
107  Parallel::ReductionData<
109  get<basis_history_tag>(box)[0], get<operand_tag>(box))},
110  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
112  ResidualMonitor<Metavariables, FieldsTag>>(cache));
113 
114  // Terminate algorithm for now. The `ResidualMonitor` will receive the
115  // reduction that is performed above and then broadcast to the following
116  // action, which is responsible for restarting the algorithm.
117  return {std::move(box), true};
118  }
119 };
120 
121 template <typename FieldsTag>
122 struct OrthogonalizeOperand {
123  private:
124  using fields_tag = FieldsTag;
125  using operand_tag =
127  using orthogonalization_iteration_id_tag =
130  using basis_history_tag = LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>;
131 
132  public:
133  template <typename ParallelComponent, typename DbTagsList,
134  typename Metavariables, typename ArrayIndex,
135  typename DataBox = db::DataBox<DbTagsList>,
137  db::tag_is_retrievable_v<operand_tag, DataBox> and
139  orthogonalization_iteration_id_tag, DataBox> and
140  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
141  db::tag_is_retrievable_v<LinearSolver::Tags::IterationId,
142  DataBox>> = nullptr>
143  static void apply(db::DataBox<DbTagsList>& box,
145  const ArrayIndex& array_index,
146  const double orthogonalization) noexcept {
147  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
148  make_not_null(&box),
149  [orthogonalization](
151  const gsl::not_null<
153  orthogonalization_iteration_id,
155  basis_history) noexcept {
156  *operand -= orthogonalization *
157  gsl::at(basis_history, *orthogonalization_iteration_id);
158  (*orthogonalization_iteration_id)++;
159  },
160  get<basis_history_tag>(box));
161 
162  const auto& next_orthogonalization_iteration_id =
163  get<orthogonalization_iteration_id_tag>(box);
164  const auto& iteration_id = get<LinearSolver::Tags::IterationId>(box);
165 
166  if (next_orthogonalization_iteration_id <= iteration_id) {
168  StoreOrthogonalization<FieldsTag, ParallelComponent>>(
169  Parallel::ReductionData<
171  inner_product(gsl::at(get<basis_history_tag>(box),
172  next_orthogonalization_iteration_id),
173  get<operand_tag>(box))},
174  Parallel::get_parallel_component<ParallelComponent>(
175  cache)[array_index],
177  ResidualMonitor<Metavariables, FieldsTag>>(cache));
178  } else {
180  StoreFinalOrthogonalization<FieldsTag, ParallelComponent>>(
181  Parallel::ReductionData<
183  inner_product(get<operand_tag>(box), get<operand_tag>(box))},
184  Parallel::get_parallel_component<ParallelComponent>(
185  cache)[array_index],
187  ResidualMonitor<Metavariables, FieldsTag>>(cache));
188  }
189  }
190 };
191 
192 template <typename FieldsTag>
193 struct NormalizeOperandAndUpdateField {
194  private:
195  using fields_tag = FieldsTag;
196  using initial_fields_tag =
198  using operand_tag =
200  using basis_history_tag = LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>;
201 
202  public:
203  template <typename ParallelComponent, typename DbTagsList,
204  typename Metavariables, typename ArrayIndex,
205  typename DataBox = db::DataBox<DbTagsList>,
207  db::tag_is_retrievable_v<initial_fields_tag, DataBox> and
208  db::tag_is_retrievable_v<operand_tag, DataBox> and
209  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
211  DataBox>> = nullptr>
212  static void apply(db::DataBox<DbTagsList>& box,
214  const ArrayIndex& array_index, const double normalization,
215  const DenseVector<double>& minres,
217  has_converged) noexcept {
218  db::mutate<operand_tag, basis_history_tag, fields_tag,
219  LinearSolver::Tags::HasConverged>(
220  make_not_null(&box),
221  [
222  normalization, &minres, &has_converged
227  local_has_converged,
229  initial_field) noexcept {
230  *operand /= normalization;
231  basis_history->push_back(*operand);
232  *field = initial_field;
233  for (size_t i = 0; i < minres.size(); i++) {
234  *field += minres[i] * gsl::at(*basis_history, i);
235  }
236  *local_has_converged = has_converged;
237  },
238  get<initial_fields_tag>(box));
239 
240  // Proceed with algorithm
241  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
242  .perform_algorithm(true);
243  }
244 };
245 
246 } // namespace gmres_detail
247 } // 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:630
Defines functions for interfacing with the parallelization framework.
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:1097
Functionality for solving linear systems of equations.
Definition: TerminateIfConverged.hpp:22
Definition: TaggedTuple.hpp:27
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:206
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.
A set of vectors that form a basis of the -th Krylov subspace .
Definition: Tags.hpp:191
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:76
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:273
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:72
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:135
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:1628
The prefix for tags related to an orthogonalization procedurce.
Definition: Tags.hpp:155
Defines class DenseVector.
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:222
Prefix indicating the value a quantity will take on the next iteration of the algorithm.
Definition: Prefixes.hpp:150
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:475
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
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
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid...
Definition: Gsl.hpp:124