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 #include "Utilities/GetOutput.hpp"
24 
25 /// \cond
26 namespace tuples {
27 template <typename...>
28 class TaggedTuple;
29 } // namespace tuples
30 /// \endcond
31 
32 namespace NonlinearSolver::newton_raphson::detail {
33 
34 template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
35 struct CheckResidualMagnitude {
36  using fields_tag = FieldsTag;
37  using residual_tag =
39  using residual_magnitude_square_tag =
41  using initial_residual_magnitude_tag =
43  using prev_residual_magnitude_square_tag =
45 
46  template <typename ParallelComponent, typename DataBox,
47  typename Metavariables, typename ArrayIndex, typename... Args>
48  static void apply(DataBox& box, Parallel::GlobalCache<Metavariables>& cache,
49  const ArrayIndex& /*array_index*/,
50  Args&&... args) noexcept {
51  if constexpr (db::tag_is_retrievable_v<residual_magnitude_square_tag,
52  DataBox>) {
53  apply_impl(box, cache, std::forward<Args>(args)...);
54  } else {
55  ERROR(
56  "The residual monitor is not yet initialized. This is a bug, so "
57  "please file an issue.");
58  }
59  }
60 
61  private:
62  template <typename DbTagsList, 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>(iteration_id, globalization_iteration_id,
73  residual_magnitude, step_length, cache);
74 
75  if (UNLIKELY(iteration_id == 0)) {
76  db::mutate<initial_residual_magnitude_tag>(
77  make_not_null(&box),
78  [residual_magnitude](const gsl::not_null<double*>
79  initial_residual_magnitude) noexcept {
80  *initial_residual_magnitude = sqrt(residual_magnitude);
81  });
82  } else {
83  // Make sure we are converging. Far away from the solution the correction
84  // determined by the linear solve might be bad, so we employ a
85  // globalization strategy to guide the solver towards the solution when
86  // the residual doesn't decrease sufficiently. See the `NewtonRaphson`
87  // class documentation for details.
88  const double sufficient_decrease =
89  get<NonlinearSolver::Tags::SufficientDecrease<OptionsGroup>>(box);
90  const double residual_magnitude_square =
91  get<residual_magnitude_square_tag>(box);
92  const double initial_residual_magnitude =
93  get<initial_residual_magnitude_tag>(box);
94  const double abs_tolerance =
95  get<Convergence::Tags::Criteria<OptionsGroup>>(box).absolute_residual;
96  const double rel_tolerance =
97  get<Convergence::Tags::Criteria<OptionsGroup>>(box).relative_residual;
98  // This is the directional derivative of the residual magnitude square
99  // f(x) = |r(x)|^2 in the descent direction
100  const double residual_magnitude_square_slope =
101  -2. * residual_magnitude_square;
102  // Check the sufficient decrease condition. Also make sure the residual
103  // didn't hit the tolerance.
104  if (residual_magnitude > abs_tolerance and
105  residual_magnitude / initial_residual_magnitude > rel_tolerance and
106  next_residual_magnitude_square >
107  residual_magnitude_square + sufficient_decrease * step_length *
108  residual_magnitude_square_slope) {
109  // The residual didn't sufficiently decrease. Perform a globalization
110  // step.
111  if (globalization_iteration_id <
113  box)) {
114  const double next_step_length = std::clamp(
116  globalization_iteration_id, step_length,
118  residual_magnitude_square, residual_magnitude_square_slope,
119  next_residual_magnitude_square,
120  get<prev_residual_magnitude_square_tag>(box)),
121  0.1 * step_length, 0.5 * step_length);
122  db::mutate<NonlinearSolver::Tags::StepLength<OptionsGroup>,
123  prev_residual_magnitude_square_tag>(
124  make_not_null(&box),
125  [step_length, next_residual_magnitude_square](
126  const gsl::not_null<double*> prev_step_length,
128  prev_residual_magnitude_square) noexcept {
129  *prev_step_length = step_length;
130  *prev_residual_magnitude_square =
131  next_residual_magnitude_square;
132  });
133  // Do some logging
135  ::Verbosity::Verbose)) {
137  "Step with length %s didn't sufficiently decrease the '" +
138  Options::name<OptionsGroup>() +
139  "' iteration %zu residual (possible overshoot). Residual: "
140  "%s. Next step length: %s.\n",
141  get_output(step_length), iteration_id,
142  get_output(residual_magnitude), get_output(next_step_length));
144  ::Verbosity::Debug)) {
145  Parallel::printf("Residual magnitude slope: %s\n",
146  get_output(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  "WARNING: Failed to sufficiently decrease the '" +
161  Options::name<OptionsGroup>() +
162  "' iteration %zu residual in %zu globalization steps. This "
163  "is usually indicative of an ill-posed problem, for "
164  "example when the linearization of the nonlinear operator "
165  "is not computed correctly.",
166  iteration_id, globalization_iteration_id);
167  } // min_step_length
168  } // sufficient decrease condition
169  } // initial iteration
170 
171  db::mutate<residual_magnitude_square_tag>(
172  make_not_null(&box), [next_residual_magnitude_square](
174  local_residual_magnitude_square) noexcept {
175  *local_residual_magnitude_square = next_residual_magnitude_square;
176  });
177 
178  // At this point, the iteration is complete. We proceed with logging and
179  // checking convergence before broadcasting back to the elements.
180 
181  // Determine whether the nonlinear solver has converged
182  Convergence::HasConverged has_converged{
183  get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
184  residual_magnitude, get<initial_residual_magnitude_tag>(box)};
185 
186  // Do some logging
188  ::Verbosity::Quiet)) {
189  if (UNLIKELY(iteration_id == 0)) {
190  Parallel::printf("Nonlinear solver '" + Options::name<OptionsGroup>() +
191  "' initialized with residual %e.\n",
192  residual_magnitude);
193  } else {
194  Parallel::printf("Nonlinear solver '" + Options::name<OptionsGroup>() +
195  "' iteration %zu done. Remaining residual: %e\n",
196  iteration_id, 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("The nonlinear solver '" +
203  Options::name<OptionsGroup>() +
204  "' has converged without any iterations: %s\n",
205  has_converged);
206  } else {
207  Parallel::printf("The nonlinear solver '" +
208  Options::name<OptionsGroup>() +
209  "' has converged in %zu iterations: %s\n",
210  iteration_id, has_converged);
211  }
212  }
213 
214  Parallel::receive_data<Tags::GlobalizationResult<OptionsGroup>>(
215  Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
217  // NOLINTNEXTLINE(performance-move-const-arg)
218  std::move(has_converged)));
219  }
220 };
221 
222 } // 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(const size_t globalization_iteration_id, const double step_length, const double prev_step_length, const double residual, const double residual_slope, const double next_residual, const double prev_residual) noexcept
Find the next step length for the line-search globalization.
Definition: LineSearch.cpp:12
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
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
std::sqrt
T sqrt(T... args)
get_output
std::string get_output(const T &t) noexcept
Get the streamed output of t as a std::string
Definition: GetOutput.hpp:14
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
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
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