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 <type_traits>
9 #include <utility>
10 
12 #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
14 #include "NumericalAlgorithms/Convergence/Tags.hpp"
16 #include "Parallel/GlobalCache.hpp"
17 #include "Parallel/Invoke.hpp"
18 #include "Parallel/Reduction.hpp"
19 #include "ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp"
20 #include "ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp"
22 #include "Utilities/Functional.hpp"
23 #include "Utilities/Gsl.hpp"
25 #include "Utilities/Requires.hpp"
26 #include "Utilities/TMPL.hpp"
27 
28 /// \cond
29 namespace tuples {
30 template <typename...>
31 class TaggedTuple;
32 } // namespace tuples
33 namespace LinearSolver::gmres::detail {
34 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
35 struct ResidualMonitor;
36 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
37  typename Label>
38 struct NormalizeOperandAndUpdateField;
39 } // namespace LinearSolver::gmres::detail
40 /// \endcond
41 
42 namespace LinearSolver::gmres::detail {
43 
44 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
45  typename Label, typename SourceTag>
46 struct PrepareSolve {
47  private:
48  using fields_tag = FieldsTag;
49  using initial_fields_tag = db::add_tag_prefix<::Tags::Initial, fields_tag>;
50  using source_tag = SourceTag;
51  using operator_applied_to_fields_tag =
53  using operand_tag =
55  using basis_history_tag =
57 
58  public:
59  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
60  typename ArrayIndex, typename ActionList,
61  typename ParallelComponent>
64  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
66  const ArrayIndex& array_index, const ActionList /*meta*/,
67  const ParallelComponent* const /*meta*/) noexcept {
68  db::mutate<Convergence::Tags::IterationId<OptionsGroup>, operand_tag,
69  initial_fields_tag, basis_history_tag>(
70  make_not_null(&box),
71  [](const gsl::not_null<size_t*> iteration_id, const auto operand,
72  const auto initial_fields, const auto basis_history,
73  const auto& source, const auto& operator_applied_to_fields,
74  const auto& fields) noexcept {
75  *iteration_id = 0;
76  *operand = source - operator_applied_to_fields;
77  *initial_fields = fields;
78  *basis_history = typename basis_history_tag::type{};
79  },
80  get<source_tag>(box), get<operator_applied_to_fields_tag>(box),
81  get<fields_tag>(box));
82 
83  Parallel::contribute_to_reduction<InitializeResidualMagnitude<
84  FieldsTag, OptionsGroup, ParallelComponent>>(
85  Parallel::ReductionData<
87  inner_product(get<operand_tag>(box), get<operand_tag>(box))},
88  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
90  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
91 
92  if constexpr (Preconditioned) {
93  using preconditioned_operand_tag =
95  using preconditioned_basis_history_tag =
97 
98  db::mutate<preconditioned_basis_history_tag>(
99  make_not_null(&box),
100  [](const auto preconditioned_basis_history) noexcept {
101  *preconditioned_basis_history =
102  typename preconditioned_basis_history_tag::type{};
103  });
104  }
105 
106  return {std::move(box)};
107  }
108 };
109 
110 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
111  typename Label>
112 struct NormalizeInitialOperand {
113  private:
114  using fields_tag = FieldsTag;
115  using operand_tag =
117  using basis_history_tag =
119 
120  public:
121  using inbox_tags = tmpl::list<Tags::InitialOrthogonalization<OptionsGroup>>;
122 
123  template <typename DbTags, typename... InboxTags, typename Metavariables,
124  typename ArrayIndex>
125  static bool is_ready(const db::DataBox<DbTags>& box,
126  const tuples::TaggedTuple<InboxTags...>& inboxes,
127  const Parallel::GlobalCache<Metavariables>& /*cache*/,
128  const ArrayIndex& /*array_index*/) noexcept {
129  const auto& inbox =
130  get<Tags::InitialOrthogonalization<OptionsGroup>>(inboxes);
132  box)) != inbox.end();
133  }
134 
135  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
136  typename ArrayIndex, typename ActionList,
137  typename ParallelComponent>
138  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
140  const Parallel::GlobalCache<Metavariables>& /*cache*/,
141  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
142  const ParallelComponent* const /*meta*/) noexcept {
143  auto received_data = std::move(
144  tuples::get<Tags::InitialOrthogonalization<OptionsGroup>>(inboxes)
146  .mapped());
147  const double residual_magnitude = get<0>(received_data);
148  auto& has_converged = get<1>(received_data);
149 
150  db::mutate<operand_tag, basis_history_tag,
152  make_not_null(&box), [residual_magnitude, &has_converged](
153  const auto operand, const auto basis_history,
155  local_has_converged) noexcept {
156  *operand /= residual_magnitude;
157  basis_history->push_back(*operand);
158  *local_has_converged = std::move(has_converged);
159  });
160 
161  // Skip steps entirely if the solve has already converged
162  constexpr size_t step_end_index = tmpl::index_of<
163  ActionList, NormalizeOperandAndUpdateField<
164  FieldsTag, OptionsGroup, Preconditioned, Label>>::value;
165  constexpr size_t this_action_index =
166  tmpl::index_of<ActionList, NormalizeInitialOperand>::value;
167  return {std::move(box), false,
168  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
169  ? (step_end_index + 1)
170  : (this_action_index + 1)};
171  }
172 };
173 
174 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
175  typename Label>
176 struct PrepareStep {
177  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
178  typename ArrayIndex, typename ActionList,
179  typename ParallelComponent>
182  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
183  const Parallel::GlobalCache<Metavariables>& /*cache*/,
184  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
185  const ParallelComponent* const /*meta*/) noexcept {
186  if constexpr (Preconditioned) {
187  using fields_tag = FieldsTag;
188  using operand_tag =
190  using preconditioned_operand_tag =
192 
193  db::mutate<preconditioned_operand_tag>(
194  make_not_null(&box),
195  [](const auto preconditioned_operand, const auto& operand) noexcept {
196  // Start the preconditioner at zero because we have no reason to
197  // expect the remaining residual to have a particular form.
198  // Another possibility would be to start the preconditioner with an
199  // initial guess equal to its source, so not running the
200  // preconditioner at all means it is the identity, but that approach
201  // appears to yield worse results.
202  *preconditioned_operand =
203  make_with_value<typename preconditioned_operand_tag::type>(
204  operand, 0.);
205  },
206  get<operand_tag>(box));
207  }
208  return {std::move(box)};
209  }
210 };
211 
212 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
213  typename Label>
214 struct PerformStep {
215  private:
216  using fields_tag = FieldsTag;
217  using operand_tag =
219  using preconditioned_operand_tag =
221 
222  public:
223  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
224  typename ArrayIndex, typename ActionList,
225  typename ParallelComponent>
228  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
230  const ArrayIndex& array_index, const ActionList /*meta*/,
231  const ParallelComponent* const /*meta*/) noexcept {
232  using operator_tag = db::add_tag_prefix<
234  std::conditional_t<Preconditioned, preconditioned_operand_tag,
235  operand_tag>>;
236  using orthogonalization_iteration_id_tag =
239  using basis_history_tag =
241 
242  if constexpr (Preconditioned) {
243  using preconditioned_basis_history_tag =
245 
246  db::mutate<preconditioned_basis_history_tag>(
247  make_not_null(&box),
248  [](const auto preconditioned_basis_history,
249  const auto& preconditioned_operand) noexcept {
250  preconditioned_basis_history->push_back(preconditioned_operand);
251  },
252  get<preconditioned_operand_tag>(box));
253  }
254 
255  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
256  make_not_null(&box),
257  [](const auto operand,
258  const gsl::not_null<size_t*> orthogonalization_iteration_id,
259  const auto& operator_action) noexcept {
260  *operand = typename operand_tag::type(operator_action);
261  *orthogonalization_iteration_id = 0;
262  },
263  get<operator_tag>(box));
264 
266  StoreOrthogonalization<FieldsTag, OptionsGroup, ParallelComponent>>(
267  Parallel::ReductionData<
271  get<Convergence::Tags::IterationId<OptionsGroup>>(box),
272  get<orthogonalization_iteration_id_tag>(box),
273  inner_product(get<basis_history_tag>(box)[0],
274  get<operand_tag>(box))},
275  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
277  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
278 
279  return {std::move(box)};
280  }
281 };
282 
283 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
284  typename Label>
285 struct OrthogonalizeOperand {
286  private:
287  using fields_tag = FieldsTag;
288  using operand_tag =
290  using orthogonalization_iteration_id_tag =
293  using basis_history_tag =
295 
296  public:
297  using inbox_tags = tmpl::list<Tags::Orthogonalization<OptionsGroup>>;
298 
299  template <typename DbTags, typename... InboxTags, typename Metavariables,
300  typename ArrayIndex>
301  static bool is_ready(const db::DataBox<DbTags>& box,
302  const tuples::TaggedTuple<InboxTags...>& inboxes,
303  const Parallel::GlobalCache<Metavariables>& /*cache*/,
304  const ArrayIndex& /*array_index*/) noexcept {
305  const auto& inbox = get<Tags::Orthogonalization<OptionsGroup>>(inboxes);
307  box)) != inbox.end();
308  }
309 
310  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
311  typename ArrayIndex, typename ActionList,
312  typename ParallelComponent>
313  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
316  const ArrayIndex& array_index, const ActionList /*meta*/,
317  const ParallelComponent* const /*meta*/) noexcept {
318  const double orthogonalization = std::move(
319  tuples::get<Tags::Orthogonalization<OptionsGroup>>(inboxes)
321  .mapped());
322 
323  db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
324  make_not_null(&box),
325  [orthogonalization](
326  const auto operand,
327  const gsl::not_null<size_t*> orthogonalization_iteration_id,
328  const auto& basis_history) noexcept {
329  *operand -= orthogonalization *
330  gsl::at(basis_history, *orthogonalization_iteration_id);
331  ++(*orthogonalization_iteration_id);
332  },
333  get<basis_history_tag>(box));
334 
335  const auto& next_orthogonalization_iteration_id =
336  get<orthogonalization_iteration_id_tag>(box);
337  const auto& iteration_id =
338  get<Convergence::Tags::IterationId<OptionsGroup>>(box);
339  const bool orthogonalization_complete =
340  next_orthogonalization_iteration_id == iteration_id + 1;
341  const double local_orthogonalization =
342  inner_product(orthogonalization_complete
343  ? get<operand_tag>(box)
344  : gsl::at(get<basis_history_tag>(box),
345  next_orthogonalization_iteration_id),
346  get<operand_tag>(box));
347 
349  StoreOrthogonalization<FieldsTag, OptionsGroup, ParallelComponent>>(
350  Parallel::ReductionData<
354  iteration_id, next_orthogonalization_iteration_id,
355  local_orthogonalization},
356  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
358  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
359 
360  // Repeat this action until orthogonalization is complete
361  constexpr size_t this_action_index =
362  tmpl::index_of<ActionList, OrthogonalizeOperand>::value;
363  return {std::move(box), false,
364  orthogonalization_complete ? (this_action_index + 1)
365  : this_action_index};
366  }
367 };
368 
369 template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
370  typename Label>
371 struct NormalizeOperandAndUpdateField {
372  private:
373  using fields_tag = FieldsTag;
374  using initial_fields_tag = db::add_tag_prefix<::Tags::Initial, fields_tag>;
375  using operand_tag =
377  using preconditioned_operand_tag =
379  using basis_history_tag =
381  using preconditioned_basis_history_tag =
383  Preconditioned, preconditioned_operand_tag, operand_tag>>;
384 
385  public:
386  using inbox_tags = tmpl::list<Tags::FinalOrthogonalization<OptionsGroup>>;
387 
388  template <typename DbTags, typename... InboxTags, typename Metavariables,
389  typename ArrayIndex>
390  static bool is_ready(const db::DataBox<DbTags>& box,
391  const tuples::TaggedTuple<InboxTags...>& inboxes,
392  const Parallel::GlobalCache<Metavariables>& /*cache*/,
393  const ArrayIndex& /*array_index*/) noexcept {
394  const auto& inbox =
395  get<Tags::FinalOrthogonalization<OptionsGroup>>(inboxes);
397  box)) != inbox.end();
398  }
399 
400  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
401  typename ArrayIndex, typename ActionList,
402  typename ParallelComponent>
403  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
405  const Parallel::GlobalCache<Metavariables>& /*cache*/,
406  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
407  const ParallelComponent* const /*meta*/) noexcept {
408  // Retrieve reduction data from inbox
409  auto received_data = std::move(
410  tuples::get<Tags::FinalOrthogonalization<OptionsGroup>>(inboxes)
412  .mapped());
413  const double normalization = get<0>(received_data);
414  const auto& minres = get<1>(received_data);
415  auto& has_converged = get<2>(received_data);
416 
417  db::mutate<operand_tag, basis_history_tag, fields_tag,
420  make_not_null(&box),
421  [normalization, &minres, &has_converged](
422  const auto operand, const auto basis_history, const auto field,
423  const gsl::not_null<size_t*> iteration_id,
424  const gsl::not_null<Convergence::HasConverged*> local_has_converged,
425  const auto& initial_field,
426  const auto& preconditioned_basis_history) noexcept {
427  // Avoid an FPE if the new operand norm is exactly zero. In that case
428  // the problem is solved and the algorithm will terminate (see
429  // Proposition 9.3 in \cite Saad2003). Since there will be no next
430  // iteration we don't need to normalize the operand.
431  if (LIKELY(normalization > 0.)) {
432  *operand /= normalization;
433  }
434  basis_history->push_back(*operand);
435  *field = initial_field;
436  for (size_t i = 0; i < minres.size(); i++) {
437  *field += minres[i] * gsl::at(preconditioned_basis_history, i);
438  }
439  ++(*iteration_id);
440  *local_has_converged = std::move(has_converged);
441  },
442  get<initial_fields_tag>(box),
443  get<preconditioned_basis_history_tag>(box));
444 
445  // Repeat steps until the solve has converged
446  constexpr size_t this_action_index =
447  tmpl::index_of<ActionList, NormalizeOperandAndUpdateField>::value;
448  constexpr size_t prepare_step_index =
449  tmpl::index_of<ActionList, PrepareStep<FieldsTag, OptionsGroup,
450  Preconditioned, Label>>::value;
451  return {std::move(box), false,
452  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
453  ? (this_action_index + 1)
454  : prepare_step_index};
455  }
456 };
457 
458 } // namespace LinearSolver::gmres::detail
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
utility
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
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
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:55
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:223
LinearSolver::Tags::KrylovSubspaceBasis
A set of vectors that form a basis of the -th Krylov subspace .
Definition: Tags.hpp:163
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:52
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:1116
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:1003
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:245
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:1424
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: ActionTesting.hpp:2053
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
MakeWithValue.hpp
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
Gsl.hpp
LinearSolver::Tags::Orthogonalization
The prefix for tags related to an orthogonalization procedure.
Definition: Tags.hpp:128
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
std::conditional_t
Convergence::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the iterative algorithm has converged,...
Definition: Tags.hpp:86
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
type_traits
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183