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 "IO/Logging/Tags.hpp"
16 #include "IO/Logging/Verbosity.hpp"
17 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
18 #include "NumericalAlgorithms/Convergence/Tags.hpp"
20 #include "Parallel/AlgorithmMetafunctions.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/GetOutput.hpp"
33 #include "Utilities/Gsl.hpp"
35 #include "Utilities/TMPL.hpp"
36 #include "Utilities/TaggedTuple.hpp"
37 
38 /// \cond
39 namespace NonlinearSolver::newton_raphson::detail {
40 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
41 struct ResidualMonitor;
42 template <typename FieldsTag, typename OptionsGroup, typename Label>
43 struct PrepareStep;
44 template <typename FieldsTag, typename OptionsGroup, typename Label>
45 struct CompleteStep;
46 } // namespace NonlinearSolver::newton_raphson::detail
47 /// \endcond
48 
49 namespace NonlinearSolver::newton_raphson::detail {
50 
51 using ResidualReductionData = Parallel::ReductionData<
52  // Iteration ID
54  // Globalization iteration ID
56  // Residual magnitude square
58  // Step length
60 
61 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
62 struct InitializeElement {
63  private:
64  using fields_tag = FieldsTag;
65  using source_tag = SourceTag;
66  using nonlinear_operator_applied_to_fields_tag =
68  using correction_tag =
70  using globalization_fields_tag =
72 
73  public:
74  using simple_tags =
75  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
77  nonlinear_operator_applied_to_fields_tag, correction_tag,
81  globalization_fields_tag>;
82  using compute_tags = tmpl::list<
84 
85  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
86  typename ArrayIndex, typename ActionList,
87  typename ParallelComponent>
88  static auto apply(db::DataBox<DbTagsList>& box,
89  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
90  const Parallel::GlobalCache<Metavariables>& /*cache*/,
91  const ArrayIndex& /*array_index*/,
92  const ActionList /*meta*/,
93  const ParallelComponent* const /*meta*/) noexcept {
95  tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
102  return std::make_tuple(std::move(box));
103  }
104 };
105 
106 // Reset to the initial state of the algorithm. To determine the initial
107 // residual magnitude and to check if the algorithm has already converged, we
108 // perform a global reduction to the `ResidualMonitor`.
109 template <typename FieldsTag, typename OptionsGroup, typename Label>
110 struct PrepareSolve {
111  private:
112  using fields_tag = FieldsTag;
113  using nonlinear_residual_tag =
115 
116  public:
117  using const_global_cache_tags =
118  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
119 
120  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
121  typename ArrayIndex, typename ActionList,
122  typename ParallelComponent>
123  static std::tuple<db::DataBox<DbTagsList>&&> apply(
124  db::DataBox<DbTagsList>& box,
125  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
127  const ArrayIndex& array_index, const ActionList /*meta*/,
128  const ParallelComponent* const /*meta*/) noexcept {
130  ::Verbosity::Debug)) {
131  Parallel::printf("%s %s: Prepare solve\n", get_output(array_index),
132  Options::name<OptionsGroup>());
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 DbTagsList, typename... InboxTags, typename Metavariables,
166  typename ArrayIndex, typename ActionList,
167  typename ParallelComponent>
169  size_t>
170  apply(db::DataBox<DbTagsList>& box,
172  const Parallel::GlobalCache<Metavariables>& /*cache*/,
173  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
174  const ParallelComponent* const /*meta*/) noexcept {
175  auto& inbox = get<Tags::GlobalizationResult<OptionsGroup>>(inboxes);
177  box)) == inbox.end()) {
178  return {std::move(box), Parallel::AlgorithmExecution::Retry,
180  }
181 
182  // Retrieve reduction data from inbox
183  auto globalization_result = std::move(
184  inbox
186  .mapped());
187  ASSERT(
188  std::holds_alternative<Convergence::HasConverged>(globalization_result),
189  "No globalization should occur for the initial residual. This is a "
190  "bug, so please file an issue.");
191  auto& has_converged = get<Convergence::HasConverged>(globalization_result);
192 
193  db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
194  make_not_null(&box),
195  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
196  local_has_converged) noexcept {
197  *local_has_converged = std::move(has_converged);
198  });
199 
200  // Skip steps entirely if the solve has already converged
201  constexpr size_t complete_step_index =
202  tmpl::index_of<ActionList,
203  CompleteStep<FieldsTag, OptionsGroup, Label>>::value;
204  constexpr size_t this_action_index =
205  tmpl::index_of<ActionList, ReceiveInitialHasConverged>::value;
206  return {std::move(box), Parallel::AlgorithmExecution::Continue,
207  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
208  ? complete_step_index
209  : (this_action_index + 1)};
210  }
211 };
212 
213 // Prepare the next Newton-Raphson step. In particular, we prepare the DataBox
214 // for the linear solver which will run after this action to solve the
215 // linearized problem for the `correction_tag`. The source for the linear
216 // solver is the residual, which is in the DataBox as a compute tag so needs no
217 // preparation. We only need to set the initial guess.
218 //
219 // We also prepare the line-search globalization here. Since we don't know if
220 // a step will be sufficient before taking it, we have to store the original
221 // field values.
222 //
223 // The algorithm jumps back to this action from `CompleteStep` to continue
224 // iterating the nonlinear solve.
225 template <typename FieldsTag, typename OptionsGroup, typename Label>
226 struct PrepareStep {
227  private:
228  using fields_tag = FieldsTag;
229  using correction_tag =
231  using linear_operator_applied_to_correction_tag =
233  using globalization_fields_tag =
235 
236  public:
237  using const_global_cache_tags =
238  tmpl::list<NonlinearSolver::Tags::DampingFactor<OptionsGroup>,
240 
241  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
242  typename ArrayIndex, typename ActionList,
243  typename ParallelComponent>
244  static std::tuple<db::DataBox<DbTagsList>&&> apply(
245  db::DataBox<DbTagsList>& box,
246  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
247  const Parallel::GlobalCache<Metavariables>& /*cache*/,
248  const ArrayIndex& array_index, const ActionList /*meta*/,
249  const ParallelComponent* const /*meta*/) noexcept {
251  ::Verbosity::Debug)) {
253  "%s %s(%zu): Prepare step\n", get_output(array_index),
254  Options::name<OptionsGroup>(),
256  }
257 
258  db::mutate<Convergence::Tags::IterationId<OptionsGroup>, correction_tag,
259  linear_operator_applied_to_correction_tag,
263  globalization_fields_tag>(
264  make_not_null(&box),
265  [](const gsl::not_null<size_t*> iteration_id, const auto correction,
266  const auto linear_operator_applied_to_correction,
267  const gsl::not_null<size_t*> globalization_iteration_id,
268  const gsl::not_null<double*> step_length,
269  const auto globalization_fields, const auto& fields,
270  const double damping_factor) noexcept {
271  ++(*iteration_id);
272  // Begin the linear solve with a zero initial guess
273  *correction =
274  make_with_value<typename correction_tag::type>(fields, 0.);
275  // Since the initial guess is zero, we don't need to apply the linear
276  // operator to it but can just set it to zero as well. Linear things
277  // are nice :)
278  *linear_operator_applied_to_correction = make_with_value<
279  typename linear_operator_applied_to_correction_tag::type>(fields,
280  0.);
281  // Prepare line search globalization
282  *globalization_iteration_id = 0;
283  *step_length = damping_factor;
284  *globalization_fields = fields;
285  },
286  db::get<fields_tag>(box),
287  db::get<NonlinearSolver::Tags::DampingFactor<OptionsGroup>>(box));
288  return {std::move(box)};
289  }
290 };
291 
292 // Between `PrepareStep` and this action the linear solver has run, so the
293 // `correction_tag` is now filled with the solution to the linearized problem.
294 // We just take a step in the direction of the correction.
295 //
296 // The `Globalize` action will jump back here in case the step turned out to not
297 // decrease the residual sufficiently. It will have adjusted the step length, so
298 // we can just try again with a step of that length.
299 template <typename FieldsTag, typename OptionsGroup, typename Label>
300 struct PerformStep {
301  private:
302  using fields_tag = FieldsTag;
303  using correction_tag =
305  using globalization_fields_tag =
307 
308  public:
309  using const_global_cache_tags =
310  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
311 
312  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
313  typename ArrayIndex, typename ActionList,
314  typename ParallelComponent>
315  static std::tuple<db::DataBox<DbTagsList>&&> apply(
316  db::DataBox<DbTagsList>& box,
317  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
318  const Parallel::GlobalCache<Metavariables>& /*cache*/,
319  const ArrayIndex& array_index, const ActionList /*meta*/,
320  const ParallelComponent* const /*meta*/) noexcept {
322  ::Verbosity::Debug)) {
324  "%s %s(%zu): Perform step with length: %g\n", get_output(array_index),
325  Options::name<OptionsGroup>(),
328  }
329 
330  // Apply the correction that the linear solve has determined to attempt
331  // improving the nonlinear solution
332  db::mutate<fields_tag>(
333  make_not_null(&box),
334  [](const auto fields, const auto& correction, const double step_length,
335  const auto& globalization_fields) {
336  *fields = globalization_fields + step_length * correction;
337  },
338  db::get<correction_tag>(box),
339  db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box),
340  db::get<globalization_fields_tag>(box));
341 
342  return {std::move(box)};
343  }
344 };
345 
346 // Between `PerformStep` and this action the nonlinear operator has been applied
347 // to the updated fields. The residual is up-to-date because it is a compute
348 // tag, so at this point we need to check if the step has sufficiently reduced
349 // the residual. We perform a global reduction to the `ResidualMonitor` for this
350 // purpose.
351 template <typename FieldsTag, typename OptionsGroup, typename Label>
352 struct ContributeToResidualMagnitudeReduction {
353  private:
354  using fields_tag = FieldsTag;
355  using nonlinear_residual_tag =
357 
358  public:
359  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
360  typename ArrayIndex, typename ActionList,
361  typename ParallelComponent>
362  static std::tuple<db::DataBox<DbTagsList>&&> apply(
363  db::DataBox<DbTagsList>& box,
364  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
366  const ArrayIndex& array_index, const ActionList /*meta*/,
367  const ParallelComponent* const /*meta*/) noexcept {
368  const auto& residual = db::get<nonlinear_residual_tag>(box);
369  const double local_residual_magnitude_square =
370  LinearSolver::inner_product(residual, residual);
371  ResidualReductionData reduction_data{
372  db::get<Convergence::Tags::IterationId<OptionsGroup>>(box),
375  local_residual_magnitude_square,
376  db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box)};
378  CheckResidualMagnitude<FieldsTag, OptionsGroup, ParallelComponent>>(
379  std::move(reduction_data),
380  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
382  ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
383  return {std::move(box)};
384  }
385 };
386 
387 // Wait for the `ResidualMonitor` to broadcast whether or not it has determined
388 // the step has sufficiently decreased the residual. If it has, we just proceed
389 // to `CompleteStep`. If it hasn't, the `ResidualMonitor` has also sent the new
390 // step length along, so we mutate it in the DataBox and jump back to
391 // `PerformStep` to try again with the updated step length.
392 template <typename FieldsTag, typename OptionsGroup, typename Label>
393 struct Globalize {
394  using const_global_cache_tags =
395  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
396  using inbox_tags = tmpl::list<Tags::GlobalizationResult<OptionsGroup>>;
397 
398  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
399  typename ArrayIndex, typename ActionList,
400  typename ParallelComponent>
402  size_t>
403  apply(db::DataBox<DbTagsList>& box,
405  const Parallel::GlobalCache<Metavariables>& /*cache*/,
406  const ArrayIndex& array_index, const ActionList /*meta*/,
407  const ParallelComponent* const /*meta*/) noexcept {
408  auto& inbox = get<Tags::GlobalizationResult<OptionsGroup>>(inboxes);
410  box)) == inbox.end()) {
411  return {std::move(box), Parallel::AlgorithmExecution::Retry,
413  }
414 
415  // Retrieve reduction data from inbox
416  auto globalization_result = std::move(
417  inbox
419  .mapped());
420 
421  if (std::holds_alternative<double>(globalization_result)) {
423  ::Verbosity::Debug)) {
425  "%s %s(%zu): Globalize(%zu)\n", get_output(array_index),
426  Options::name<OptionsGroup>(),
430  }
431 
432  // Update the step length
433  db::mutate<NonlinearSolver::Tags::StepLength<OptionsGroup>,
436  make_not_null(&box),
437  [&globalization_result](const gsl::not_null<double*> step_length,
439  globalization_iteration_id) noexcept {
440  *step_length = get<double>(globalization_result);
441  ++(*globalization_iteration_id);
442  });
443  // Continue globalization by taking a step with the updated step length
444  // and checking the residual again
445  constexpr size_t perform_step_index =
446  tmpl::index_of<ActionList,
447  PerformStep<FieldsTag, OptionsGroup, Label>>::value;
448  return {std::move(box), Parallel::AlgorithmExecution::Continue,
449  perform_step_index};
450  }
451 
452  // At this point globalization is complete, so we proceed with the algorithm
453  auto& has_converged = get<Convergence::HasConverged>(globalization_result);
454 
455  db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
456  make_not_null(&box),
457  [&has_converged](const gsl::not_null<Convergence::HasConverged*>
458  local_has_converged) noexcept {
459  *local_has_converged = std::move(has_converged);
460  });
461 
462  constexpr size_t this_action_index =
463  tmpl::index_of<ActionList, Globalize>::value;
464  return {std::move(box), Parallel::AlgorithmExecution::Continue,
465  this_action_index + 1};
466  }
467 };
468 
469 // Jump back to `PrepareStep` to continue iterating if the algorithm has not yet
470 // converged, or complete the solve and proceed with the action list if it has
471 // converged. This is a separate action because the user has the opportunity to
472 // insert actions before the step completes, for example to do observations.
473 template <typename FieldsTag, typename OptionsGroup, typename Label>
474 struct CompleteStep {
475  using const_global_cache_tags =
476  tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
477 
478  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
479  typename ArrayIndex, typename ActionList,
480  typename ParallelComponent>
481  static std::tuple<db::DataBox<DbTagsList>&&, bool, size_t> apply(
482  db::DataBox<DbTagsList>& box,
483  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
484  const Parallel::GlobalCache<Metavariables>& /*cache*/,
485  const ArrayIndex& array_index, const ActionList /*meta*/,
486  const ParallelComponent* const /*meta*/) noexcept {
488  ::Verbosity::Debug)) {
490  "%s %s(%zu): Complete step\n", get_output(array_index),
491  Options::name<OptionsGroup>(),
493  }
494 
495  // Repeat steps until the solve has converged
496  constexpr size_t prepare_step_index =
497  tmpl::index_of<ActionList,
498  PrepareStep<FieldsTag, OptionsGroup, Label>>::value;
499  constexpr size_t this_action_index =
500  tmpl::index_of<ActionList, CompleteStep>::value;
501  return {std::move(box), false,
502  get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
503  ? (this_action_index + 1)
504  : prepare_step_index};
505  }
506 };
507 
508 } // namespace NonlinearSolver::newton_raphson::detail
utility
Parallel::AlgorithmExecution::Retry
@ Retry
Temporarily stop executing iterable actions, but try the same action again after receiving data from ...
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
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:65
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:535
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:791
NonlinearSolver::Tags::StepLength
The length of nonlinear solver steps.
Definition: Tags.hpp:182
cmath
Tags.hpp
get_output
std::string get_output(const T &t) noexcept
Get the streamed output of t as a std::string
Definition: GetOutput.hpp:14
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:39
NonlinearSolver::Tags::ResidualCompute
Compute the residual from the SourceTag and the db::add_tag_prefix<NonlinearSolver::Tags::OperatorA...
Definition: Tags.hpp:154
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
Parallel::AlgorithmExecution::Continue
@ Continue
Leave the algorithm termination flag in its current state.
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)
Parallel::AlgorithmExecution
AlgorithmExecution
The possible options for altering the current execution of the algorithm, used in the return type of ...
Definition: AlgorithmMetafunctions.hpp:31
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
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
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)
Parallel::contribute_to_reduction
void contribute_to_reduction(ReductionData< Ts... > reduction_data, const SenderProxy &sender_component, const TargetProxy &target_component, [[maybe_unused]] const gsl::not_null< SectionType * > section=&no_section()) noexcept
Perform a reduction from the sender_component (typically your own parallel component) to the target_c...
Definition: Reduction.hpp:272
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: ReadSpecPiecewisePolynomial.hpp:13
variant