Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cmath>
8 : #include <cstddef>
9 : #include <string>
10 : #include <utility>
11 : #include <variant>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 : #include "DataStructures/DataBox/Prefixes.hpp"
16 : #include "IO/Logging/Verbosity.hpp"
17 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
18 : #include "NumericalAlgorithms/Convergence/Reason.hpp"
19 : #include "Parallel/GlobalCache.hpp"
20 : #include "Parallel/Invoke.hpp"
21 : #include "Parallel/Printf/Printf.hpp"
22 : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/LineSearch.hpp"
23 : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Tags/InboxTags.hpp"
24 : #include "ParallelAlgorithms/NonlinearSolver/Observe.hpp"
25 : #include "Utilities/Gsl.hpp"
26 : #include "Utilities/PrettyType.hpp"
27 :
28 : /// \cond
29 : namespace Convergence::Tags {
30 : template <typename OptionsGroup>
31 : struct Criteria;
32 : } // namespace Convergence::Tags
33 : namespace LinearSolver::Tags {
34 : template <typename Tag>
35 : struct Magnitude;
36 : template <typename Tag>
37 : struct MagnitudeSquare;
38 : } // namespace LinearSolver::Tags
39 : namespace NonlinearSolver::Tags {
40 : template <typename Tag>
41 : struct Globalization;
42 : template <typename OptionsGroup>
43 : struct MaxGlobalizationSteps;
44 : template <typename Tag>
45 : struct Residual;
46 : template <typename OptionsGroup>
47 : struct StepLength;
48 : template <typename OptionsGroup>
49 : struct SufficientDecrease;
50 : } // namespace NonlinearSolver::Tags
51 : namespace logging::Tags {
52 : template <typename OptionsGroup>
53 : struct Verbosity;
54 : } // namespace logging::Tags
55 : namespace tuples {
56 : template <typename...>
57 : class TaggedTuple;
58 : } // namespace tuples
59 : /// \endcond
60 :
61 : namespace NonlinearSolver::newton_raphson::detail {
62 :
63 : template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
64 : struct CheckResidualMagnitude {
65 : using fields_tag = FieldsTag;
66 : using residual_tag =
67 : db::add_tag_prefix<NonlinearSolver::Tags::Residual, fields_tag>;
68 : using residual_magnitude_square_tag =
69 : LinearSolver::Tags::MagnitudeSquare<residual_tag>;
70 : using initial_residual_magnitude_tag =
71 : ::Tags::Initial<LinearSolver::Tags::Magnitude<residual_tag>>;
72 : using prev_residual_magnitude_square_tag =
73 : NonlinearSolver::Tags::Globalization<residual_magnitude_square_tag>;
74 :
75 : template <typename ParallelComponent, typename DataBox,
76 : typename Metavariables, typename ArrayIndex, typename... Args>
77 : static void apply(DataBox& box, Parallel::GlobalCache<Metavariables>& cache,
78 : const ArrayIndex& /*array_index*/,
79 : const size_t iteration_id,
80 : const size_t globalization_iteration_id,
81 : const double next_residual_magnitude_square,
82 : const double step_length) {
83 : const double residual_magnitude = sqrt(next_residual_magnitude_square);
84 :
85 : NonlinearSolver::observe_detail::contribute_to_reduction_observer<
86 : OptionsGroup, ParallelComponent>(
87 : iteration_id, globalization_iteration_id, residual_magnitude,
88 : step_length, cache);
89 :
90 : if (UNLIKELY(iteration_id == 0)) {
91 : db::mutate<initial_residual_magnitude_tag>(
92 : [residual_magnitude](
93 : const gsl::not_null<double*> initial_residual_magnitude) {
94 : *initial_residual_magnitude = residual_magnitude;
95 : },
96 : make_not_null(&box));
97 : } else {
98 : // Make sure we are converging. Far away from the solution the correction
99 : // determined by the linear solve might be bad, so we employ a
100 : // globalization strategy to guide the solver towards the solution when
101 : // the residual doesn't decrease sufficiently. See the `NewtonRaphson`
102 : // class documentation for details.
103 : const double sufficient_decrease =
104 : get<NonlinearSolver::Tags::SufficientDecrease<OptionsGroup>>(box);
105 : const double residual_magnitude_square =
106 : get<residual_magnitude_square_tag>(box);
107 : const double initial_residual_magnitude =
108 : get<initial_residual_magnitude_tag>(box);
109 : const double abs_tolerance =
110 : get<Convergence::Tags::Criteria<OptionsGroup>>(box).absolute_residual;
111 : const double rel_tolerance =
112 : get<Convergence::Tags::Criteria<OptionsGroup>>(box).relative_residual;
113 : // This is the directional derivative of the residual magnitude square
114 : // f(x) = |r(x)|^2 in the descent direction
115 : const double residual_magnitude_square_slope =
116 : -2. * residual_magnitude_square;
117 : // Check the sufficient decrease condition. Also make sure the residual
118 : // didn't hit the tolerance.
119 : if (residual_magnitude > abs_tolerance and
120 : residual_magnitude / initial_residual_magnitude > rel_tolerance and
121 : next_residual_magnitude_square >
122 : residual_magnitude_square + sufficient_decrease * step_length *
123 : residual_magnitude_square_slope) {
124 : // The residual didn't sufficiently decrease. Perform a globalization
125 : // step.
126 : if (globalization_iteration_id <
127 : get<NonlinearSolver::Tags::MaxGlobalizationSteps<OptionsGroup>>(
128 : box)) {
129 : const double next_step_length = std::clamp(
130 : NonlinearSolver::newton_raphson::next_step_length(
131 : globalization_iteration_id, step_length,
132 : get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box),
133 : residual_magnitude_square, residual_magnitude_square_slope,
134 : next_residual_magnitude_square,
135 : get<prev_residual_magnitude_square_tag>(box)),
136 : 0.1 * step_length, 0.5 * step_length);
137 : db::mutate<NonlinearSolver::Tags::StepLength<OptionsGroup>,
138 : prev_residual_magnitude_square_tag>(
139 : [step_length, next_residual_magnitude_square](
140 : const gsl::not_null<double*> prev_step_length,
141 : const gsl::not_null<double*> prev_residual_magnitude_square) {
142 : *prev_step_length = step_length;
143 : *prev_residual_magnitude_square =
144 : next_residual_magnitude_square;
145 : },
146 : make_not_null(&box));
147 : // Do some logging
148 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
149 : ::Verbosity::Verbose)) {
150 : Parallel::printf(
151 : "%s(%zu): Step with length %g didn't sufficiently decrease the "
152 : "residual (possible overshoot). Residual: %e. Next step "
153 : "length: %g.\n",
154 : pretty_type::name<OptionsGroup>(), iteration_id, step_length,
155 : residual_magnitude, next_step_length);
156 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
157 : ::Verbosity::Debug)) {
158 : Parallel::printf("Residual magnitude slope: %e\n",
159 : residual_magnitude_square_slope);
160 : }
161 : }
162 : // Broadcast back to the elements signaling that they should perform a
163 : // globalization step, then return early.
164 : Parallel::receive_data<Tags::GlobalizationResult<OptionsGroup>>(
165 : Parallel::get_parallel_component<BroadcastTarget>(cache),
166 : iteration_id,
167 : std::variant<double, Convergence::HasConverged>{
168 : next_step_length});
169 : return;
170 : } else {
171 : // We have performed the maximum number of globalization steps without
172 : // sufficiently decreasing the residual. The step length is so small
173 : // now that the algorithm won't be converging anymore. Treat this as
174 : // an error and stop the Newton-Raphson algorithm.
175 : Convergence::HasConverged convergence_error{
176 : Convergence::Reason::Error,
177 : "Failed to sufficiently decrease the residual in " +
178 : std::to_string(globalization_iteration_id) +
179 : " globalization steps. This is usually indicative of an "
180 : "ill-posed problem, for example when the linearization of "
181 : "the nonlinear operator is not computed correctly.",
182 : iteration_id};
183 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
184 : ::Verbosity::Quiet)) {
185 : Parallel::printf("%s(%zu): WARNING: %s\n",
186 : pretty_type::name<OptionsGroup>(), iteration_id,
187 : convergence_error.error_message());
188 : }
189 : Parallel::receive_data<Tags::GlobalizationResult<OptionsGroup>>(
190 : Parallel::get_parallel_component<BroadcastTarget>(cache),
191 : iteration_id,
192 : std::variant<double, Convergence::HasConverged>{
193 : std::move(convergence_error)});
194 : return;
195 : } // min_step_length
196 : } // sufficient decrease condition
197 : } // initial iteration
198 :
199 : db::mutate<residual_magnitude_square_tag>(
200 : [next_residual_magnitude_square](
201 : const gsl::not_null<double*> local_residual_magnitude_square) {
202 : *local_residual_magnitude_square = next_residual_magnitude_square;
203 : },
204 : make_not_null(&box));
205 :
206 : // At this point, the iteration is complete. We proceed with logging and
207 : // checking convergence before broadcasting back to the elements.
208 :
209 : // Determine whether the nonlinear solver has converged
210 : Convergence::HasConverged has_converged{
211 : get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
212 : residual_magnitude, get<initial_residual_magnitude_tag>(box)};
213 :
214 : // Do some logging
215 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
216 : ::Verbosity::Quiet)) {
217 : if (UNLIKELY(iteration_id == 0)) {
218 : Parallel::printf("%s initialized with residual: %e\n",
219 : pretty_type::name<OptionsGroup>(), residual_magnitude);
220 : } else {
221 : Parallel::printf(
222 : "%s(%zu) iteration complete (%zu globalization steps, step length "
223 : "%g). Remaining residual: %e\n",
224 : pretty_type::name<OptionsGroup>(), iteration_id,
225 : globalization_iteration_id, step_length, residual_magnitude);
226 : }
227 : }
228 : if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
229 : box) >= ::Verbosity::Quiet)) {
230 : if (UNLIKELY(iteration_id == 0)) {
231 : Parallel::printf("%s has converged without any iterations: %s\n",
232 : pretty_type::name<OptionsGroup>(), has_converged);
233 : } else {
234 : Parallel::printf("%s has converged in %zu iterations: %s\n",
235 : pretty_type::name<OptionsGroup>(), iteration_id,
236 : has_converged);
237 : }
238 : }
239 :
240 : Parallel::receive_data<Tags::GlobalizationResult<OptionsGroup>>(
241 : Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
242 : std::variant<double, Convergence::HasConverged>(
243 : // NOLINTNEXTLINE(performance-move-const-arg)
244 : std::move(has_converged)));
245 : }
246 : };
247 :
248 : } // namespace NonlinearSolver::newton_raphson::detail
|
