ElementActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cmath>
7 #include <cstddef>
8 #include <limits>
9 #include <tuple>
10 #include <utility>
11 #include <variant>
12 
14 #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 #include "Informer/Tags.hpp"
16 #include "Informer/Verbosity.hpp"
17 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
18 #include "NumericalAlgorithms/Convergence/Tags.hpp"
20 #include "Parallel/Actions/SetupDataBox.hpp"
21 #include "Parallel/GlobalCache.hpp"
22 #include "Parallel/Invoke.hpp"
23 #include "Parallel/Printf.hpp"
24 #include "Parallel/Reduction.hpp"
25 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
27 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/ResidualMonitorActions.hpp"
28 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Tags/InboxTags.hpp"
31 #include "Utilities/Functional.hpp"
32 #include "Utilities/Gsl.hpp"
34 #include "Utilities/TMPL.hpp"
35 #include "Utilities/TaggedTuple.hpp"
36 
37 /// \cond
38 namespace NonlinearSolver::newton_raphson::detail {
39 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
40 struct ResidualMonitor;
41 template <typename FieldsTag, typename OptionsGroup, typename Label>
42 struct PrepareStep;
43 template <typename FieldsTag, typename OptionsGroup, typename Label>
44 struct Globalize;
45 } // namespace NonlinearSolver::newton_raphson::detail
46 /// \endcond
47 
48 namespace NonlinearSolver::newton_raphson::detail {
49 
50 using ResidualReductionData = Parallel::ReductionData<
51  // Iteration ID
53  // Globalization iteration ID
55  // Residual magnitude square
57  // Step length
59 
60 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
61 struct InitializeElement {
62  private:
63  using fields_tag = FieldsTag;
64  using source_tag = SourceTag;
65  using nonlinear_operator_applied_to_fields_tag =
67  using correction_tag =
69  using globalization_fields_tag =
71 
72  public:
73  using simple_tags =
74  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
76  nonlinear_operator_applied_to_fields_tag, correction_tag,
80  globalization_fields_tag>;
81  using compute_tags = tmpl::list<
83 
84  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
85  typename ArrayIndex, typename ActionList,
86  typename ParallelComponent>
87  static auto apply(db::DataBox<DbTagsList>& box,
88  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
89  const Parallel::GlobalCache<Metavariables>& /*cache*/,
90  const ArrayIndex& /*array_index*/,
91  const ActionList /*meta*/,
92  const ParallelComponent* const /*meta*/) noexcept {
94  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
101  return std::make_tuple(std::move(box));
102  }
103 };
104 
105 // Reset to the initial state of the algorithm. To determine the initial
106 // residual magnitude and to check if the algorithm has already converged, we
107 // perform a global reduction to the `ResidualMonitor`.
108 template <typename FieldsTag, typename OptionsGroup, typename Label>
109 struct PrepareSolve {
110  private:
111  using fields_tag = FieldsTag;
112  using nonlinear_residual_tag =
114 
115  public:
116  using const_global_cache_tags =
117  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
118 
119  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
120  typename ArrayIndex, typename ActionList,
121  typename ParallelComponent>
122  static std::tuple<db::DataBox<DbTagsList>&&> apply(
123  db::DataBox<DbTagsList>& box,
124  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
126  const ArrayIndex& array_index, const ActionList /*meta*/,
127  const ParallelComponent* const /*meta*/) noexcept {
129  ::Verbosity::Debug)) {
131  "%s " + Options::name<OptionsGroup>() + ": Prepare solve\n",
132  array_index);
133  }
134 
135  db::mutate<Convergence::Tags::IterationId<OptionsGroup>>(
136  make_not_null(&box),
137  [](const gsl::not_null<size_t*> iteration_id) noexcept {
138  *iteration_id = 0;
139  });
140 
141  // Perform a global reduction to compute the initial residual magnitude
142  const auto& residual = db::get<nonlinear_residual_tag>(box);
143  const double local_residual_magnitude_square =
144  LinearSolver::inner_product(residual, residual);
145  ResidualReductionData reduction_data{0, 0, local_residual_magnitude_square,
146  1.};
148  CheckResidualMagnitude<FieldsTag, OptionsGroup, ParallelComponent>>(
149  std::move(reduction_data),
150  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
152  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
153 
154  return {std::move(box)};
155  }
156 };
157 
158 // Wait for the broadcast from the `ResidualMonitor` to complete the preparation
159 // for the solve. We skip the solve altogether if the algorithm has already
160 // converged.
161 template <typename FieldsTag, typename OptionsGroup, typename Label>
162 struct ReceiveInitialHasConverged {
163  using inbox_tags = tmpl::list<Tags::GlobalizationResult<OptionsGroup>>;
164 
165  template <typename DbTags, typename... InboxTags, typename Metavariables,
166  typename ArrayIndex>
167  static bool is_ready(const db::DataBox<DbTags>& box,
168  const tuples::TaggedTuple<InboxTags...>& inboxes,
169  const Parallel::GlobalCache<Metavariables>& /*cache*/,
170  const ArrayIndex& /*array_index*/) noexcept {
171  const auto& inbox = get<Tags::GlobalizationResult<OptionsGroup>>(inboxes);
173  box)) != inbox.end();
174  }
175 
176  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
177  typename ArrayIndex, typename ActionList,
178  typename ParallelComponent>
179  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
180  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
181  const Parallel::GlobalCache<Metavariables>& /*cache*/,
182  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
183  const ParallelComponent* const /*meta*/) noexcept {
184  // Retrieve reduction data from inbox
185  auto globalization_result = std::move(
186  tuples::get<Tags::GlobalizationResult<OptionsGroup>>(inboxes)
188  .mapped());
189  ASSERT(
190  std::holds_alternative<Convergence::HasConverged>(globalization_result),
191  "No globalization should occur for the initial residual. This is a "
192  "bug, so please file an issue.");
193  auto& has_converged = get<Convergence::HasConverged>(globalization_result);
194 
195  db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
196  make_not_null(&box),
197  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
198  local_has_converged) noexcept {
199  *local_has_converged = std::move(has_converged);
200  });
201 
202  // Skip steps entirely if the solve has already converged
203  constexpr size_t complete_step_index =
204  tmpl::index_of<ActionList,
205  Globalize<FieldsTag, OptionsGroup, Label>>::value +
206  1;
207  constexpr size_t this_action_index =
208  tmpl::index_of<ActionList, ReceiveInitialHasConverged>::value;
209  return {std::move(box), false,
210  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
211  ? complete_step_index
212  : (this_action_index + 1)};
213  }
214 };
215 
216 // Prepare the next Newton-Raphson step. In particular, we prepare the DataBox
217 // for the linear solver which will run after this action to solve the
218 // linearized problem for the `correction_tag`. The source for the linear
219 // solver is the residual, which is in the DataBox as a compute tag so needs no
220 // preparation. We only need to set the initial guess.
221 //
222 // We also prepare the line-search globalization here. Since we don't know if
223 // a step will be sufficient before taking it, we have to store the original
224 // field values.
225 //
226 // The algorithm jumps back to this action from `CompleteStep` to continue
227 // iterating the nonlinear solve.
228 template <typename FieldsTag, typename OptionsGroup, typename Label>
229 struct PrepareStep {
230  private:
231  using fields_tag = FieldsTag;
232  using correction_tag =
234  using linear_operator_applied_to_correction_tag =
236  using globalization_fields_tag =
238 
239  public:
240  using const_global_cache_tags =
241  tmpl::list<NonlinearSolver::Tags::DampingFactor<OptionsGroup>,
243 
244  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
245  typename ArrayIndex, typename ActionList,
246  typename ParallelComponent>
247  static std::tuple<db::DataBox<DbTagsList>&&> apply(
248  db::DataBox<DbTagsList>& box,
249  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
250  const Parallel::GlobalCache<Metavariables>& /*cache*/,
251  const ArrayIndex& array_index, const ActionList /*meta*/,
252  const ParallelComponent* const /*meta*/) noexcept {
254  ::Verbosity::Debug)) {
256  "%s " + Options::name<OptionsGroup>() + "(%zu): Prepare step\n",
257  array_index,
259  }
260 
261  db::mutate<Convergence::Tags::IterationId<OptionsGroup>, correction_tag,
262  linear_operator_applied_to_correction_tag,
266  globalization_fields_tag>(
267  make_not_null(&box),
268  [](const gsl::not_null<size_t*> iteration_id, const auto correction,
269  const auto linear_operator_applied_to_correction,
270  const gsl::not_null<size_t*> globalization_iteration_id,
271  const gsl::not_null<double*> step_length,
272  const auto globalization_fields, const auto& fields,
273  const double damping_factor) noexcept {
274  ++(*iteration_id);
275  // Begin the linear solve with a zero initial guess
276  *correction =
277  make_with_value<typename correction_tag::type>(fields, 0.);
278  // Since the initial guess is zero, we don't need to apply the linear
279  // operator to it but can just set it to zero as well. Linear things
280  // are nice :)
281  *linear_operator_applied_to_correction = make_with_value<
282  typename linear_operator_applied_to_correction_tag::type>(fields,
283  0.);
284  // Prepare line search globalization
285  *globalization_iteration_id = 0;
286  *step_length = damping_factor;
287  *globalization_fields = fields;
288  },
289  db::get<fields_tag>(box),
290  db::get<NonlinearSolver::Tags::DampingFactor<OptionsGroup>>(box));
291  return {std::move(box)};
292  }
293 };
294 
295 // Between `PrepareStep` and this action the linear solver has run, so the
296 // `correction_tag` is now filled with the solution to the linearized problem.
297 // We just take a step in the direction of the correction.
298 //
299 // The `Globalize` action will jump back here in case the step turned out to not
300 // decrease the residual sufficiently. It will have adjusted the step length, so
301 // we can just try again with a step of that length.
302 template <typename FieldsTag, typename OptionsGroup, typename Label>
303 struct PerformStep {
304  private:
305  using fields_tag = FieldsTag;
306  using correction_tag =
308  using globalization_fields_tag =
310 
311  public:
312  using const_global_cache_tags =
313  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
314 
315  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
316  typename ArrayIndex, typename ActionList,
317  typename ParallelComponent>
318  static std::tuple<db::DataBox<DbTagsList>&&> apply(
319  db::DataBox<DbTagsList>& box,
320  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
321  const Parallel::GlobalCache<Metavariables>& /*cache*/,
322  const ArrayIndex& array_index, const ActionList /*meta*/,
323  const ParallelComponent* const /*meta*/) noexcept {
325  ::Verbosity::Debug)) {
327  "%s " + Options::name<OptionsGroup>() +
328  "(%zu): Perform step with length %f\n",
329  array_index,
332  }
333 
334  // Apply the correction that the linear solve has determined to attempt
335  // improving the nonlinear solution
336  db::mutate<fields_tag>(
337  make_not_null(&box),
338  [](const auto fields, const auto& correction, const double step_length,
339  const auto& globalization_fields) {
340  *fields = globalization_fields + step_length * correction;
341  },
342  db::get<correction_tag>(box),
343  db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box),
344  db::get<globalization_fields_tag>(box));
345 
346  return {std::move(box)};
347  }
348 };
349 
350 // Between `PerformStep` and this action the nonlinear operator has been applied
351 // to the updated fields. The residual is up-to-date because it is a compute
352 // tag, so at this point we need to check if the step has sufficiently reduced
353 // the residual. We perform a global reduction to the `ResidualMonitor` for this
354 // purpose.
355 template <typename FieldsTag, typename OptionsGroup, typename Label>
356 struct ContributeToResidualMagnitudeReduction {
357  private:
358  using fields_tag = FieldsTag;
359  using nonlinear_residual_tag =
361 
362  public:
363  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
364  typename ArrayIndex, typename ActionList,
365  typename ParallelComponent>
366  static std::tuple<db::DataBox<DbTagsList>&&> apply(
367  db::DataBox<DbTagsList>& box,
368  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
370  const ArrayIndex& array_index, const ActionList /*meta*/,
371  const ParallelComponent* const /*meta*/) noexcept {
372  const auto& residual = db::get<nonlinear_residual_tag>(box);
373  const double local_residual_magnitude_square =
374  LinearSolver::inner_product(residual, residual);
375  ResidualReductionData reduction_data{
376  db::get<Convergence::Tags::IterationId<OptionsGroup>>(box),
379  local_residual_magnitude_square,
380  db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box)};
382  CheckResidualMagnitude<FieldsTag, OptionsGroup, ParallelComponent>>(
383  std::move(reduction_data),
384  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
386  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
387  return {std::move(box)};
388  }
389 };
390 
391 // Wait for the `ResidualMonitor` to broadcast whether or not it has determined
392 // the step has sufficiently decreased the residual. If it has, we just proceed
393 // to `CompleteStep`. If it hasn't, the `ResidualMonitor` has also sent the new
394 // step length along, so we mutate it in the DataBox and jump back to
395 // `PerformStep` to try again with the updated step length.
396 template <typename FieldsTag, typename OptionsGroup, typename Label>
397 struct Globalize {
398  using const_global_cache_tags =
399  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
400  using inbox_tags = tmpl::list<Tags::GlobalizationResult<OptionsGroup>>;
401 
402  template <typename DbTags, typename... InboxTags, typename Metavariables,
403  typename ArrayIndex>
404  static bool is_ready(const db::DataBox<DbTags>& box,
405  const tuples::TaggedTuple<InboxTags...>& inboxes,
406  const Parallel::GlobalCache<Metavariables>& /*cache*/,
407  const ArrayIndex& /*array_index*/) noexcept {
408  const auto& inbox = get<Tags::GlobalizationResult<OptionsGroup>>(inboxes);
410  box)) != inbox.end();
411  }
412 
413  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
414  typename ArrayIndex, typename ActionList,
415  typename ParallelComponent>
416  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
417  db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
418  const Parallel::GlobalCache<Metavariables>& /*cache*/,
419  const ArrayIndex& array_index, const ActionList /*meta*/,
420  const ParallelComponent* const /*meta*/) noexcept {
421  // Retrieve reduction data from inbox
422  auto globalization_result = std::move(
423  tuples::get<Tags::GlobalizationResult<OptionsGroup>>(inboxes)
425  .mapped());
426 
427  if (std::holds_alternative<double>(globalization_result)) {
429  ::Verbosity::Debug)) {
431  "%s " + Options::name<OptionsGroup>() + "(%zu): Globalize(%zu)\n",
432  array_index,
436  }
437 
438  // Update the step length
439  db::mutate<NonlinearSolver::Tags::StepLength<OptionsGroup>,
442  make_not_null(&box),
443  [&globalization_result](const gsl::not_null<double*> step_length,
445  globalization_iteration_id) noexcept {
446  *step_length = get<double>(globalization_result);
447  ++(*globalization_iteration_id);
448  });
449  // Continue globalization by taking a step with the updated step length
450  // and checking the residual again
451  constexpr size_t perform_step_index =
452  tmpl::index_of<ActionList,
453  PerformStep<FieldsTag, OptionsGroup, Label>>::value;
454  return {std::move(box), false, perform_step_index};
455  }
456 
457  // At this point globalization is complete, so we proceed with the algorithm
458  auto& has_converged = get<Convergence::HasConverged>(globalization_result);
459 
460  db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
461  make_not_null(&box),
462  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
463  local_has_converged) noexcept {
464  *local_has_converged = std::move(has_converged);
465  });
466 
467  constexpr size_t this_action_index =
468  tmpl::index_of<ActionList, Globalize>::value;
469  return {std::move(box), false, this_action_index + 1};
470  }
471 };
472 
473 // Jump back to `PrepareStep` to continue iterating if the algorithm has not yet
474 // converged, or complete the solve and proceed with the action list if it has
475 // converged. This is a separate action because the user has the opportunity to
476 // insert actions before the step completes, for example to do observations.
477 template <typename FieldsTag, typename OptionsGroup, typename Label>
478 struct CompleteStep {
479  using const_global_cache_tags =
480  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
481 
482  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
483  typename ArrayIndex, typename ActionList,
484  typename ParallelComponent>
485  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
486  db::DataBox<DbTagsList>& box,
487  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
488  const Parallel::GlobalCache<Metavariables>& /*cache*/,
489  const ArrayIndex& array_index, const ActionList /*meta*/,
490  const ParallelComponent* const /*meta*/) noexcept {
492  ::Verbosity::Debug)) {
494  "%s " + Options::name<OptionsGroup>() + "(%zu): Complete step\n",
495  array_index,
497  }
498 
499  // Repeat steps until the solve has converged
500  constexpr size_t prepare_step_index =
501  tmpl::index_of<ActionList,
502  PrepareStep<FieldsTag, OptionsGroup, Label>>::value;
503  constexpr size_t this_action_index =
504  tmpl::index_of<ActionList, CompleteStep>::value;
505  return {std::move(box), false,
506  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
507  ? (this_action_index + 1)
508  : prepare_step_index};
509  }
510 };
511 
512 } // namespace NonlinearSolver::newton_raphson::detail
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:63
Tags.hpp
InnerProduct.hpp
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:521
Parallel::printf
void printf(const std::string &format, Args &&... args)
Print an atomic message to stdout with C printf usage.
Definition: Printf.hpp:103
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:969
NonlinearSolver::Tags::StepLength
The length of nonlinear solver steps.
Definition: Tags.hpp:182
cmath
Tags.hpp
Initialization::mutate_assign
constexpr void mutate_assign(const gsl::not_null< db::DataBox< BoxTags > * > box, Args &&... args) noexcept
Perform a mutation to the DataBox box, assigning the args to the tags in MutateTagList in order.
Definition: MutateAssign.hpp:40
NonlinearSolver::Tags::ResidualCompute
Compute the residual from the SourceTag and the db::add_tag_prefix<NonlinearSolver::Tags::OperatorA...
Definition: Tags.hpp:154
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
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: MockRuntimeSystemFreeFunctions.hpp:155
Printf.hpp
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
make_with_value
std::remove_const_t< R > make_with_value(const T &input, const ValueType &value) noexcept
Given an object of type T, create an object of type R whose elements are initialized to value.
Definition: MakeWithValue.hpp:88
cstddef
logging::Tags::Verbosity
Tag for putting Verbosity in a DataBox.
Definition: Tags.hpp:33
Assert.hpp
MakeWithValue.hpp
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Convergence::Tags::IterationId
Identifies a step in an iterative algorithm.
Definition: Tags.hpp:74
limits
Gsl.hpp
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::numeric_limits::max
T max(T... args)
Convergence::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the iterative algorithm has converged,...
Definition: Tags.hpp:86
NonlinearSolver::Tags::Globalization
Prefix indicating the Tag is related to the globalization procedure.
Definition: Tags.hpp:240
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
variant