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"
14 #include "Parallel/Info.hpp"
15 #include "Parallel/Invoke.hpp"
16 #include "Parallel/Reduction.hpp"
17 #include "ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp"
19 #include "Utilities/Functional.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/Requires.hpp"
22 
23 /// \cond
24 namespace tuples {
25 template <typename...>
26 class TaggedTuple;
27 } // namespace tuples
28 namespace LinearSolver {
29 namespace gmres_detail {
30 template <typename Metavariables, typename FieldsTag>
31 struct ResidualMonitor;
32 } // namespace gmres_detail
33 } // namespace LinearSolver
34 /// \endcond
35 
36 namespace LinearSolver {
37 namespace gmres_detail {
38 
39 struct PrepareStep {
40  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
41  typename ArrayIndex, typename ActionList,
42  typename ParallelComponent>
45  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
47  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
48  const ParallelComponent* const /*meta*/) noexcept {
49  using orthogonalization_iteration_id_tag =
52 
53  db::mutate<LinearSolver::Tags::IterationId,
54  orthogonalization_iteration_id_tag>(
55  make_not_null(&box),
57  iteration_id,
58  const gsl::not_null<
60  orthogonalization_iteration_id,
62  LinearSolver::Tags::IterationId>>& next_iteration_id) noexcept {
63  *iteration_id = next_iteration_id;
64  *orthogonalization_iteration_id = 0;
65  },
66  get<::Tags::Next<LinearSolver::Tags::IterationId>>(box));
67  return {std::move(box)};
68  }
69 };
70 
71 template <typename FieldsTag>
72 struct PerformStep {
73  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
74  typename ArrayIndex, typename ActionList,
75  typename ParallelComponent>
78  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
80  const ArrayIndex& array_index,
81  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
82  const ActionList /*meta*/,
83  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
84  const ParallelComponent* const /*meta*/) noexcept {
85  using fields_tag = FieldsTag;
86  using operand_tag =
88  using operator_tag =
90  using basis_history_tag =
92 
93  ASSERT(get<LinearSolver::Tags::IterationId>(box) !=
95  "Linear solve iteration ID is at initial state. Did you forget to "
96  "invoke 'PrepareStep'?");
97 
98  db::mutate<operand_tag>(
99  make_not_null(&box),
100  [](const gsl::not_null<db::item_type<operand_tag>*> operand,
101  const db::const_item_type<operator_tag>& operator_action) noexcept {
102  *operand = db::item_type<operand_tag>(operator_action);
103  },
104  get<operator_tag>(box));
105 
107  StoreOrthogonalization<FieldsTag, ParallelComponent>>(
108  Parallel::ReductionData<
110  get<basis_history_tag>(box)[0], get<operand_tag>(box))},
111  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
113  ResidualMonitor<Metavariables, FieldsTag>>(cache));
114 
115  // Terminate algorithm for now. The `ResidualMonitor` will receive the
116  // reduction that is performed above and then broadcast to the following
117  // action, which is responsible for restarting the algorithm.
118  return {std::move(box), true};
119  }
120 };
121 
122 template <typename FieldsTag>
123 struct OrthogonalizeOperand {
124  private:
125  using fields_tag = FieldsTag;
126  using operand_tag =
128  using orthogonalization_iteration_id_tag =
131  using basis_history_tag = LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>;
132 
133  public:
134  template <typename ParallelComponent, typename DbTagsList,
135  typename Metavariables, typename ArrayIndex,
136  typename DataBox = db::DataBox<DbTagsList>,
138  db::tag_is_retrievable_v<operand_tag, DataBox> and
140  orthogonalization_iteration_id_tag, DataBox> and
141  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
142  db::tag_is_retrievable_v<LinearSolver::Tags::IterationId,
143  DataBox>> = nullptr>
144  static void apply(db::DataBox<DbTagsList>& box,
146  const ArrayIndex& array_index,
147  const double orthogonalization) noexcept {
148  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
149  make_not_null(&box),
150  [orthogonalization](
152  const gsl::not_null<
154  orthogonalization_iteration_id,
156  basis_history) noexcept {
157  *operand -= orthogonalization *
158  gsl::at(basis_history, *orthogonalization_iteration_id);
159  (*orthogonalization_iteration_id)++;
160  },
161  get<basis_history_tag>(box));
162 
163  const auto& next_orthogonalization_iteration_id =
164  get<orthogonalization_iteration_id_tag>(box);
165  const auto& iteration_id = get<LinearSolver::Tags::IterationId>(box);
166 
167  if (next_orthogonalization_iteration_id <= iteration_id) {
169  StoreOrthogonalization<FieldsTag, ParallelComponent>>(
170  Parallel::ReductionData<
172  inner_product(gsl::at(get<basis_history_tag>(box),
173  next_orthogonalization_iteration_id),
174  get<operand_tag>(box))},
175  Parallel::get_parallel_component<ParallelComponent>(
176  cache)[array_index],
178  ResidualMonitor<Metavariables, FieldsTag>>(cache));
179  } else {
181  StoreFinalOrthogonalization<FieldsTag, ParallelComponent>>(
182  Parallel::ReductionData<
184  inner_product(get<operand_tag>(box), get<operand_tag>(box))},
185  Parallel::get_parallel_component<ParallelComponent>(
186  cache)[array_index],
188  ResidualMonitor<Metavariables, FieldsTag>>(cache));
189  }
190  }
191 };
192 
193 template <typename FieldsTag>
194 struct NormalizeOperandAndUpdateField {
195  private:
196  using fields_tag = FieldsTag;
197  using initial_fields_tag =
199  using operand_tag =
201  using basis_history_tag = LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>;
202 
203  public:
204  template <typename ParallelComponent, typename DbTagsList,
205  typename Metavariables, typename ArrayIndex,
206  typename DataBox = db::DataBox<DbTagsList>,
208  db::tag_is_retrievable_v<initial_fields_tag, DataBox> and
209  db::tag_is_retrievable_v<operand_tag, DataBox> and
210  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
212  DataBox>> = nullptr>
213  static void apply(db::DataBox<DbTagsList>& box,
215  const ArrayIndex& array_index, const double normalization,
216  const DenseVector<double>& minres,
218  has_converged) noexcept {
219  db::mutate<operand_tag, basis_history_tag, fields_tag,
220  LinearSolver::Tags::HasConverged>(
221  make_not_null(&box),
222  [
223  normalization, &minres, &has_converged
228  local_has_converged,
230  initial_field) noexcept {
231  *operand /= normalization;
232  basis_history->push_back(*operand);
233  *field = initial_field;
234  for (size_t i = 0; i < minres.size(); i++) {
235  *field += minres[i] * gsl::at(*basis_history, i);
236  }
237  *local_has_converged = has_converged;
238  },
239  get<initial_fields_tag>(box));
240 
241  // Proceed with algorithm
242  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
243  .perform_algorithm(true);
244  }
245 };
246 
247 } // namespace gmres_detail
248 } // 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.
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:1020
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.
A set of vectors that form a basis of the -th Krylov subspace .
Definition: Tags.hpp:194
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
The prefix for tags related to an orthogonalization procedurce.
Definition: Tags.hpp:158
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:223
Prefix indicating the value a quantity will take on the next iteration of the algorithm.
Definition: Prefixes.hpp:152
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:246
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
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