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"
22 #include "Utilities/Requires.hpp"
23 
24 /// \cond
25 namespace tuples {
26 template <typename...>
27 class TaggedTuple;
28 } // namespace tuples
29 namespace LinearSolver::gmres::detail {
30 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
31 struct ResidualMonitor;
32 } // namespace LinearSolver::gmres::detail
33 /// \endcond
34 
35 namespace LinearSolver::gmres::detail {
36 
37 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned>
38 struct PrepareSolve {
39  private:
40  using fields_tag = FieldsTag;
41  using initial_fields_tag =
44  using operator_applied_to_fields_tag =
46  using operand_tag =
48  using basis_history_tag =
50 
51  public:
52  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
53  typename ArrayIndex, typename ActionList,
54  typename ParallelComponent>
57  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
59  const ArrayIndex& array_index, const ActionList /*meta*/,
60  const ParallelComponent* const /*meta*/) noexcept {
61  db::mutate<LinearSolver::Tags::IterationId<OptionsGroup>, operand_tag,
62  initial_fields_tag, basis_history_tag>(
63  make_not_null(&box),
64  [](const gsl::not_null<size_t*> iteration_id, const auto operand,
65  const auto initial_fields, const auto basis_history,
66  const auto& source, const auto& operator_applied_to_fields,
67  const auto& fields) noexcept {
68  *iteration_id = 0;
69  *operand = source - operator_applied_to_fields;
70  *initial_fields = fields;
71  *basis_history = db::item_type<basis_history_tag>{};
72  },
73  get<source_tag>(box), get<operator_applied_to_fields_tag>(box),
74  get<fields_tag>(box));
75 
76  Parallel::contribute_to_reduction<InitializeResidualMagnitude<
77  FieldsTag, OptionsGroup, Preconditioned, ParallelComponent>>(
78  Parallel::ReductionData<
80  inner_product(get<operand_tag>(box), get<operand_tag>(box))},
81  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
83  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
84 
85  if constexpr (Preconditioned) {
86  using preconditioned_operand_tag =
88  using preconditioned_basis_history_tag =
90 
91  db::mutate<preconditioned_basis_history_tag>(
92  make_not_null(&box),
93  [](const auto preconditioned_basis_history) noexcept {
94  *preconditioned_basis_history =
96  });
97  }
98 
99  return {
100  std::move(box),
101  // Terminate algorithm for now. The `ResidualMonitor` will receive the
102  // reduction that is performed above and then broadcast to the following
103  // action, which is responsible for restarting the algorithm.
104  true};
105  }
106 };
107 
108 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned>
109 struct NormalizeInitialOperand {
110  private:
111  using fields_tag = FieldsTag;
112  using operand_tag =
114  using basis_history_tag =
116 
117  public:
118  template <
119  typename ParallelComponent, typename DbTagsList, typename Metavariables,
120  typename ArrayIndex, typename DataBox = db::DataBox<DbTagsList>,
122  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
123  db::tag_is_retrievable_v<
125  nullptr>
126  static void apply(db::DataBox<DbTagsList>& box,
128  const ArrayIndex& array_index,
129  const double residual_magnitude,
130  const Convergence::HasConverged& has_converged) noexcept {
131  db::mutate<operand_tag, basis_history_tag,
133  make_not_null(&box), [residual_magnitude, &has_converged](
134  const auto operand, const auto basis_history,
136  local_has_converged) noexcept {
137  *operand /= residual_magnitude;
138  basis_history->push_back(*operand);
139  *local_has_converged = has_converged;
140  });
141 
142  // Proceed with algorithm
143  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
144  .perform_algorithm(true);
145  }
146 };
147 
148 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned>
149 struct PrepareStep {
150  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
151  typename ArrayIndex, typename ActionList,
152  typename ParallelComponent>
155  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
157  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
158  const ParallelComponent* const /*meta*/) noexcept {
159  if constexpr (Preconditioned) {
160  using fields_tag = FieldsTag;
161  using operand_tag =
163  using preconditioned_operand_tag =
165 
166  db::mutate<preconditioned_operand_tag>(
167  make_not_null(&box),
168  [](const auto preconditioned_operand, const auto& operand) noexcept {
169  // Start the preconditioner at zero because we have no reason to
170  // expect the remaining residual to have a particular form.
171  // Another possibility would be to start the preconditioner with an
172  // initial guess equal to its source, so not running the
173  // preconditioner at all means it is the identity, but that approach
174  // appears to yield worse results.
175  *preconditioned_operand =
176  make_with_value<db::item_type<preconditioned_operand_tag>>(
177  operand, 0.);
178  },
179  get<operand_tag>(box));
180  }
181  return {std::move(box)};
182  }
183 };
184 
185 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned>
186 struct PerformStep {
187  private:
188  using fields_tag = FieldsTag;
189  using operand_tag =
191  using preconditioned_operand_tag =
193 
194  public:
195  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
196  typename ArrayIndex, typename ActionList,
197  typename ParallelComponent>
200  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
202  const ArrayIndex& array_index,
203  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
204  const ActionList /*meta*/,
205  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
206  const ParallelComponent* const /*meta*/) noexcept {
207  using operator_tag = db::add_tag_prefix<
209  std::conditional_t<Preconditioned, preconditioned_operand_tag,
210  operand_tag>>;
211  using orthogonalization_iteration_id_tag =
214  using basis_history_tag =
216 
217  if constexpr (Preconditioned) {
218  using preconditioned_basis_history_tag =
220 
221  db::mutate<preconditioned_basis_history_tag>(
222  make_not_null(&box),
223  [](const auto preconditioned_basis_history,
224  const auto& preconditioned_operand) noexcept {
225  preconditioned_basis_history->push_back(preconditioned_operand);
226  },
227  get<preconditioned_operand_tag>(box));
228  }
229 
230  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
231  make_not_null(&box),
232  [](const auto operand,
233  const gsl::not_null<size_t*> orthogonalization_iteration_id,
234  const auto& operator_action) noexcept {
235  *operand = db::item_type<operand_tag>(operator_action);
236  *orthogonalization_iteration_id = 0;
237  },
238  get<operator_tag>(box));
239 
240  Parallel::contribute_to_reduction<StoreOrthogonalization<
241  FieldsTag, OptionsGroup, Preconditioned, ParallelComponent>>(
242  Parallel::ReductionData<
244  get<basis_history_tag>(box)[0], get<operand_tag>(box))},
245  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
247  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
248 
249  // Terminate algorithm for now. The `ResidualMonitor` will receive the
250  // reduction that is performed above and then broadcast to the following
251  // action, which is responsible for restarting the algorithm.
252  return {std::move(box), true};
253  }
254 };
255 
256 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned>
257 struct OrthogonalizeOperand {
258  private:
259  using fields_tag = FieldsTag;
260  using operand_tag =
262  using orthogonalization_iteration_id_tag =
265  using basis_history_tag =
267 
268  public:
269  template <
270  typename ParallelComponent, typename DbTagsList, typename Metavariables,
271  typename ArrayIndex, typename DataBox = db::DataBox<DbTagsList>,
273  db::tag_is_retrievable_v<operand_tag, DataBox> and
274  db::tag_is_retrievable_v<orthogonalization_iteration_id_tag,
275  DataBox> and
276  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
277  db::tag_is_retrievable_v<
279  nullptr>
280  static void apply(db::DataBox<DbTagsList>& box,
282  const ArrayIndex& array_index,
283  const double orthogonalization) noexcept {
284  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
285  make_not_null(&box),
286  [orthogonalization](
287  const auto operand,
288  const gsl::not_null<size_t*> orthogonalization_iteration_id,
289  const auto& basis_history) noexcept {
290  *operand -= orthogonalization *
291  gsl::at(basis_history, *orthogonalization_iteration_id);
292  ++(*orthogonalization_iteration_id);
293  },
294  get<basis_history_tag>(box));
295 
296  const auto& next_orthogonalization_iteration_id =
297  get<orthogonalization_iteration_id_tag>(box);
298  const auto& iteration_id =
299  get<LinearSolver::Tags::IterationId<OptionsGroup>>(box);
300 
301  if (next_orthogonalization_iteration_id <= iteration_id) {
302  Parallel::contribute_to_reduction<StoreOrthogonalization<
303  FieldsTag, OptionsGroup, Preconditioned, ParallelComponent>>(
304  Parallel::ReductionData<
306  inner_product(gsl::at(get<basis_history_tag>(box),
307  next_orthogonalization_iteration_id),
308  get<operand_tag>(box))},
309  Parallel::get_parallel_component<ParallelComponent>(
310  cache)[array_index],
312  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
313  } else {
314  Parallel::contribute_to_reduction<StoreFinalOrthogonalization<
315  FieldsTag, OptionsGroup, Preconditioned, ParallelComponent>>(
316  Parallel::ReductionData<
318  inner_product(get<operand_tag>(box), get<operand_tag>(box))},
319  Parallel::get_parallel_component<ParallelComponent>(
320  cache)[array_index],
322  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
323  }
324  }
325 };
326 
327 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned>
328 struct NormalizeOperandAndUpdateField {
329  private:
330  using fields_tag = FieldsTag;
331  using initial_fields_tag =
333  using operand_tag =
335  using preconditioned_operand_tag =
337  using basis_history_tag =
339  using preconditioned_basis_history_tag =
341  Preconditioned, preconditioned_operand_tag, operand_tag>>;
342 
343  public:
344  template <
345  typename ParallelComponent, typename DbTagsList, typename Metavariables,
346  typename ArrayIndex, typename DataBox = db::DataBox<DbTagsList>,
348  db::tag_is_retrievable_v<initial_fields_tag, DataBox> and
349  db::tag_is_retrievable_v<operand_tag, DataBox> and
350  db::tag_is_retrievable_v<basis_history_tag, DataBox> and
351  db::tag_is_retrievable_v<
353  db::tag_is_retrievable_v<
355  nullptr>
356  static void apply(db::DataBox<DbTagsList>& box,
358  const ArrayIndex& array_index, const double normalization,
359  const DenseVector<double>& minres,
360  const Convergence::HasConverged& has_converged) noexcept {
361  db::mutate<operand_tag, basis_history_tag, fields_tag,
364  make_not_null(&box),
365  [normalization, &minres, &has_converged](
366  const auto operand, const auto basis_history, const auto field,
367  const gsl::not_null<size_t*> iteration_id,
368  const gsl::not_null<Convergence::HasConverged*> local_has_converged,
369  const auto& initial_field,
370  const auto& preconditioned_basis_history) noexcept {
371  // Avoid an FPE if the new operand norm is exactly zero. In that case
372  // the problem is solved and the algorithm will terminate (see
373  // Proposition 9.3 in \cite Saad2003). Since there will be no next
374  // iteration we don't need to normalize the operand.
375  if (LIKELY(normalization > 0.)) {
376  *operand /= normalization;
377  }
378  basis_history->push_back(*operand);
379  *field = initial_field;
380  for (size_t i = 0; i < minres.size(); i++) {
381  *field += minres[i] * gsl::at(preconditioned_basis_history, i);
382  }
383  ++(*iteration_id);
384  *local_has_converged = has_converged;
385  },
386  get<initial_fields_tag>(box),
387  get<preconditioned_basis_history_tag>(box));
388 
389  // Proceed with algorithm
390  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
391  .perform_algorithm(true);
392  }
393 };
394 
395 } // namespace LinearSolver::gmres::detail
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
gsl::at
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:125
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64
InnerProduct.hpp
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:70
LinearSolver::Tags::KrylovSubspaceBasis
A set of vectors that form a basis of the -th Krylov subspace .
Definition: Tags.hpp:220
DenseVector.hpp
LinearSolver::Tags::IterationId
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:84
tuple
DenseVector< double >
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:59
Info.hpp
db::add_tag_prefix
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Definition: PrefixHelpers.hpp:145
Tags.hpp
db::mutate
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
funcl::Sqrt
Functional for computing sqrt on an object.
Definition: Functional.hpp:273
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:246
db::apply
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:1443
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
MakeWithValue.hpp
db::item_type
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
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
LinearSolver::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the linear solver has converged,...
Definition: Tags.hpp:242
Gsl.hpp
Parallel::get_parallel_component
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
LinearSolver::Tags::Orthogonalization
The prefix for tags related to an orthogonalization procedurce.
Definition: Tags.hpp:185
Requires.hpp
LIKELY
#define LIKELY(x)
Definition: Gsl.hpp:67
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
Requires
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
std::conditional_t
Prefixes.hpp
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
ConstGlobalCache.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183