ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <tuple>
7 
11 #include "NumericalAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp"
16 #include "Parallel/Info.hpp"
17 #include "Parallel/Invoke.hpp"
18 #include "Parallel/Reduction.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>
31 struct ResidualMonitor;
32 } // namespace gmres_detail
33 } // namespace LinearSolver
34 /// \endcond
35 
36 namespace LinearSolver {
37 namespace gmres_detail {
38 
39 struct NormalizeInitialOperand {
40  template <typename... DbTags, typename... InboxTags, typename Metavariables,
41  typename ArrayIndex, typename ActionList,
42  typename ParallelComponent,
43  Requires<sizeof...(DbTags) != 0> = nullptr>
44  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
47  const ArrayIndex& /*array_index*/,
48  const ActionList /*meta*/,
49  const ParallelComponent* const /*meta*/,
50  const double residual_magnitude,
52  has_converged) noexcept {
53  using fields_tag = typename Metavariables::system::fields_tag;
54  using operand_tag =
56  using basis_history_tag =
58 
59  db::mutate<operand_tag, basis_history_tag,
61  make_not_null(&box),
62  [
63  residual_magnitude, &has_converged
67  local_has_converged) noexcept {
68  *operand /= residual_magnitude;
69  basis_history->push_back(*operand);
70  *local_has_converged = has_converged;
71  });
72  }
73 };
74 
75 struct PerformStep {
76  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
77  typename ArrayIndex, typename ActionList,
78  typename ParallelComponent>
79  static auto apply(db::DataBox<DbTagsList>& box,
82  const ArrayIndex& array_index, const ActionList /*meta*/,
83  const ParallelComponent* const /*meta*/) noexcept {
84  using fields_tag = typename Metavariables::system::fields_tag;
85  using operand_tag =
87  using operator_tag =
89  using basis_history_tag =
91 
92  db::mutate<operand_tag>(
93  make_not_null(&box),
95  const db::item_type<operator_tag>& operator_action) noexcept {
96  *operand = db::item_type<operand_tag>(operator_action);
97  },
98  get<operator_tag>(box));
99 
101  StoreOrthogonalization<ParallelComponent>>(
102  Parallel::ReductionData<
104  get<basis_history_tag>(box)[0], get<operand_tag>(box))},
105  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
106  Parallel::get_parallel_component<ResidualMonitor<Metavariables>>(
107  cache));
108 
109  return std::tuple<db::DataBox<DbTagsList>&&, bool>(std::move(box), true);
110  }
111 };
112 
113 struct OrthogonalizeOperand {
114  template <typename... DbTags, typename... InboxTags, typename Metavariables,
115  typename ArrayIndex, typename ActionList,
116  typename ParallelComponent,
117  Requires<sizeof...(DbTags) != 0> = nullptr>
118  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
121  const ArrayIndex& array_index, const ActionList /*meta*/,
122  const ParallelComponent* const /*meta*/,
123  const double orthogonalization) noexcept {
124  using fields_tag = typename Metavariables::system::fields_tag;
125  using operand_tag =
127  using orthogonalization_iteration_id_tag =
130  using basis_history_tag =
132 
133  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
134  make_not_null(&box),
135  [orthogonalization](
137  const gsl::not_null<IterationId*> orthogonalization_iteration_id,
138  const db::item_type<basis_history_tag>& basis_history) noexcept {
139  *operand -= orthogonalization *
140  gsl::at(basis_history,
141  orthogonalization_iteration_id->step_number);
142  orthogonalization_iteration_id->step_number++;
143  },
144  get<basis_history_tag>(box));
145 
146  const auto next_orthogonalization_iteration_id =
147  get<orthogonalization_iteration_id_tag>(box).step_number;
148  const auto iteration_id =
149  get<LinearSolver::Tags::IterationId>(box).step_number;
150 
151  if (next_orthogonalization_iteration_id <= iteration_id) {
153  StoreOrthogonalization<ParallelComponent>>(
154  Parallel::ReductionData<
156  inner_product(gsl::at(get<basis_history_tag>(box),
157  next_orthogonalization_iteration_id),
158  get<operand_tag>(box))},
159  Parallel::get_parallel_component<ParallelComponent>(
160  cache)[array_index],
161  Parallel::get_parallel_component<ResidualMonitor<Metavariables>>(
162  cache));
163  } else {
165  StoreFinalOrthogonalization<ParallelComponent>>(
166  Parallel::ReductionData<
168  inner_product(get<operand_tag>(box), get<operand_tag>(box))},
169  Parallel::get_parallel_component<ParallelComponent>(
170  cache)[array_index],
171  Parallel::get_parallel_component<ResidualMonitor<Metavariables>>(
172  cache));
173  }
174  }
175 };
176 
177 struct NormalizeOperandAndUpdateField {
178  template <typename... DbTags, typename... InboxTags, typename Metavariables,
179  typename ArrayIndex, typename ActionList,
180  typename ParallelComponent,
181  Requires<sizeof...(DbTags) != 0> = nullptr>
182  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
185  const ArrayIndex& array_index, const ActionList /*meta*/,
186  const ParallelComponent* const /*meta*/,
187  const double normalization,
188  const DenseVector<double>& minres,
190  has_converged) noexcept {
191  using fields_tag = typename Metavariables::system::fields_tag;
192  using initial_fields_tag =
194  using operand_tag =
196  using orthogonalization_iteration_id_tag =
199  using basis_history_tag =
201 
202  db::mutate<LinearSolver::Tags::IterationId,
204  orthogonalization_iteration_id_tag, operand_tag,
205  basis_history_tag, fields_tag, LinearSolver::Tags::HasConverged>(
206  make_not_null(&box),
207  [
208  normalization, &minres, &has_converged
209  ](const gsl::not_null<IterationId*> iteration_id,
210  const gsl::not_null<IterationId*> next_iteration_id,
211  const gsl::not_null<IterationId*> orthogonalization_iteration_id,
216  local_has_converged,
217  const db::item_type<initial_fields_tag>& initial_field) noexcept {
218  iteration_id->step_number++;
219  next_iteration_id->step_number = iteration_id->step_number + 1;
220  orthogonalization_iteration_id->step_number = 0;
221  *operand /= normalization;
222  basis_history->push_back(*operand);
223  *field = initial_field;
224  for (size_t i = 0; i < minres.size(); i++) {
225  *field += minres[i] * gsl::at(*basis_history, i);
226  }
227  *local_has_converged = has_converged;
228  },
229  get<initial_fields_tag>(box));
230 
231  // Proceed with algorithm
232  // We use `ckLocal()` here since this is essentially retrieving "self",
233  // which is guaranteed to be on the local processor. This ensures the calls
234  // are evaluated in order.
235  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
236  .ckLocal()
237  ->set_terminate(false);
238  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
239  .perform_algorithm();
240  }
241 };
242 
243 } // namespace gmres_detail
244 } // 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
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
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
A set of vectors that form a basis of the -th Krylov subspace .
Definition: Tags.hpp:179
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
The prefix for tags related to an orthogonalization procedurce.
Definition: Tags.hpp:143
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: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.
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