ResidualMonitorActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <cstddef>
8 #include <variant>
9 
12 #include "Informer/Tags.hpp"
13 #include "Informer/Verbosity.hpp"
14 #include "Parallel/GlobalCache.hpp"
15 #include "Parallel/Invoke.hpp"
16 #include "Parallel/Printf.hpp"
17 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/LineSearch.hpp"
18 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Tags/InboxTags.hpp"
19 #include "ParallelAlgorithms/NonlinearSolver/Observe.hpp"
21 #include "Utilities/EqualWithinRoundoff.hpp"
22 #include "Utilities/Functional.hpp"
23 
24 /// \cond
25 namespace tuples {
26 template <typename...>
27 class TaggedTuple;
28 } // namespace tuples
29 /// \endcond
30 
31 namespace NonlinearSolver::newton_raphson::detail {
32 
33 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
34 struct CheckResidualMagnitude {
35  using fields_tag = FieldsTag;
36  using residual_tag =
38  using residual_magnitude_square_tag =
40  using initial_residual_magnitude_tag =
42  using prev_residual_magnitude_square_tag =
44 
45  template <typename ParallelComponent, typename DataBox,
46  typename Metavariables, typename ArrayIndex, typename... Args>
47  static void apply(DataBox& box, Parallel::GlobalCache<Metavariables>& cache,
48  const ArrayIndex& /*array_index*/,
49  Args&&... args) noexcept {
50  if constexpr (db::tag_is_retrievable_v<residual_magnitude_square_tag,
51  DataBox>) {
52  apply_impl<ParallelComponent>(box, cache, std::forward<Args>(args)...);
53  } else {
54  ERROR(
55  "The residual monitor is not yet initialized. This is a bug, so "
56  "please file an issue.");
57  }
58  }
59 
60  private:
61  template <typename ParallelComponent, typename DbTagsList,
62  typename Metavariables>
63  static void apply_impl(db::DataBox<DbTagsList>& box,
65  const size_t iteration_id,
66  const size_t globalization_iteration_id,
67  const double next_residual_magnitude_square,
68  const double step_length) noexcept {
69  const double residual_magnitude = sqrt(next_residual_magnitude_square);
70 
71  NonlinearSolver::observe_detail::contribute_to_reduction_observer<
72  OptionsGroup, ParallelComponent>(
73  iteration_id, globalization_iteration_id, residual_magnitude,
74  step_length, cache);
75 
76  if (UNLIKELY(iteration_id == 0)) {
77  db::mutate<initial_residual_magnitude_tag>(
78  make_not_null(&box),
79  [residual_magnitude](const gsl::not_null<double*>
80  initial_residual_magnitude) noexcept {
81  *initial_residual_magnitude = sqrt(residual_magnitude);
82  });
83  } else {
84  // Make sure we are converging. Far away from the solution the correction
85  // determined by the linear solve might be bad, so we employ a
86  // globalization strategy to guide the solver towards the solution when
87  // the residual doesn't decrease sufficiently. See the `NewtonRaphson`
88  // class documentation for details.
89  const double sufficient_decrease =
90  get<NonlinearSolver::Tags::SufficientDecrease<OptionsGroup>>(box);
91  const double residual_magnitude_square =
92  get<residual_magnitude_square_tag>(box);
93  const double initial_residual_magnitude =
94  get<initial_residual_magnitude_tag>(box);
95  const double abs_tolerance =
96  get<Convergence::Tags::Criteria<OptionsGroup>>(box).absolute_residual;
97  const double rel_tolerance =
98  get<Convergence::Tags::Criteria<OptionsGroup>>(box).relative_residual;
99  // This is the directional derivative of the residual magnitude square
100  // f(x) = |r(x)|^2 in the descent direction
101  const double residual_magnitude_square_slope =
102  -2. * residual_magnitude_square;
103  // Check the sufficient decrease condition. Also make sure the residual
104  // didn't hit the tolerance.
105  if (residual_magnitude > abs_tolerance and
106  residual_magnitude / initial_residual_magnitude > rel_tolerance and
107  next_residual_magnitude_square >
108  residual_magnitude_square + sufficient_decrease * step_length *
109  residual_magnitude_square_slope) {
110  // The residual didn't sufficiently decrease. Perform a globalization
111  // step.
112  if (globalization_iteration_id <
114  box)) {
115  const double next_step_length = std::clamp(
117  globalization_iteration_id, step_length,
119  residual_magnitude_square, residual_magnitude_square_slope,
120  next_residual_magnitude_square,
121  get<prev_residual_magnitude_square_tag>(box)),
122  0.1 * step_length, 0.5 * step_length);
123  db::mutate<NonlinearSolver::Tags::StepLength<OptionsGroup>,
124  prev_residual_magnitude_square_tag>(
125  make_not_null(&box),
126  [step_length, next_residual_magnitude_square](
127  const gsl::not_null<double*> prev_step_length,
129  prev_residual_magnitude_square) noexcept {
130  *prev_step_length = step_length;
131  *prev_residual_magnitude_square =
132  next_residual_magnitude_square;
133  });
134  // Do some logging
136  ::Verbosity::Verbose)) {
138  "%s(%zu): Step with length %g didn't sufficiently decrease the "
139  "residual (possible overshoot). Residual: %e. Next step "
140  "length: %g.\n",
141  Options::name<OptionsGroup>(), iteration_id, step_length,
142  residual_magnitude, next_step_length);
144  ::Verbosity::Debug)) {
145  Parallel::printf("Residual magnitude slope: %e\n",
146  residual_magnitude_square_slope);
147  }
148  }
149  // Broadcast back to the elements signaling that they should perform a
150  // globalization step, then return early.
151  Parallel::receive_data<Tags::GlobalizationResult<OptionsGroup>>(
152  Parallel::get_parallel_component<BroadcastTarget>(cache),
153  iteration_id,
156  return;
158  ::Verbosity::Quiet)) {
160  "%s(%zu): WARNING: Failed to sufficiently decrease the residual "
161  "in %zu globalization steps. This is usually indicative of an "
162  "ill-posed problem, for example when the linearization of the "
163  "nonlinear operator is not computed correctly.",
164  Options::name<OptionsGroup>(), iteration_id,
165  globalization_iteration_id);
166  } // min_step_length
167  } // sufficient decrease condition
168  } // initial iteration
169 
170  db::mutate<residual_magnitude_square_tag>(
171  make_not_null(&box), [next_residual_magnitude_square](
173  local_residual_magnitude_square) noexcept {
174  *local_residual_magnitude_square = next_residual_magnitude_square;
175  });
176 
177  // At this point, the iteration is complete. We proceed with logging and
178  // checking convergence before broadcasting back to the elements.
179 
180  // Determine whether the nonlinear solver has converged
181  Convergence::HasConverged has_converged{
182  get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
183  residual_magnitude, get<initial_residual_magnitude_tag>(box)};
184 
185  // Do some logging
187  ::Verbosity::Quiet)) {
188  if (UNLIKELY(iteration_id == 0)) {
189  Parallel::printf("%s initialized with residual: %e\n",
190  Options::name<OptionsGroup>(), residual_magnitude);
191  } else {
193  "%s(%zu) iteration complete (%zu globalization steps, step length "
194  "%g). Remaining residual: %e\n",
195  Options::name<OptionsGroup>(), iteration_id,
196  globalization_iteration_id, step_length, residual_magnitude);
197  }
198  }
199  if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
200  box) >= ::Verbosity::Quiet)) {
201  if (UNLIKELY(iteration_id == 0)) {
202  Parallel::printf("%s has converged without any iterations: %s\n",
203  Options::name<OptionsGroup>(), has_converged);
204  } else {
205  Parallel::printf("%s has converged in %zu iterations: %s\n",
206  Options::name<OptionsGroup>(), iteration_id,
207  has_converged);
208  }
209  }
210 
211  Parallel::receive_data<Tags::GlobalizationResult<OptionsGroup>>(
212  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
214  // NOLINTNEXTLINE(performance-move-const-arg)
215  std::move(has_converged)));
216  }
217 };
218 
219 } // namespace NonlinearSolver::newton_raphson::detail
std::apply
T apply(T... args)
LinearSolver::Tags::MagnitudeSquare
The magnitude square w.r.t. the LinearSolver::inner_product
Definition: Tags.hpp:100
NonlinearSolver::newton_raphson::next_step_length
double next_step_length(size_t globalization_iteration_id, double step_length, double prev_step_length, double residual, double residual_slope, double next_residual, double prev_residual) noexcept
Find the next step length for the line-search globalization.
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:660
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
sqrt
auto sqrt(const TensorExpression< T, typename T::type, tmpl::list<>, tmpl::list<>, tmpl::list<>> &t)
Returns the tensor expression representing the square root of a tensor expression that evaluates to a...
Definition: SquareRoot.hpp:80
Tags.hpp
Tags::Initial
Prefix indicating the initial value of a quantity.
Definition: Prefixes.hpp:85
GlobalCache.hpp
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
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:68
NonlinearSolver::Tags::StepLength
The length of nonlinear solver steps.
Definition: Tags.hpp:182
algorithm
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
Printf.hpp
DataBox.hpp
cstddef
logging::Tags::Verbosity
Tag for putting Verbosity in a DataBox.
Definition: Tags.hpp:33
NonlinearSolver::Tags::MaxGlobalizationSteps
The maximum number of allowed globalization steps.
Definition: Tags.hpp:227
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:380
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
Prefixes.hpp
NonlinearSolver::Tags::Globalization
Prefix indicating the Tag is related to the globalization procedure.
Definition: Tags.hpp:240
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
variant