Line data Source code
1 0 : // 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 <optional>
10 : #include <tuple>
11 : #include <utility>
12 : #include <variant>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 : #include "IO/Logging/Tags.hpp"
17 : #include "IO/Logging/Verbosity.hpp"
18 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
19 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
20 : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
21 : #include "Parallel/AlgorithmExecution.hpp"
22 : #include "Parallel/GetSection.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "Parallel/Invoke.hpp"
25 : #include "Parallel/Printf/Printf.hpp"
26 : #include "Parallel/Reduction.hpp"
27 : #include "Parallel/Tags/Section.hpp"
28 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
29 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
30 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
31 : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/ResidualMonitorActions.hpp"
32 : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Tags/InboxTags.hpp"
33 : #include "ParallelAlgorithms/NonlinearSolver/Tags.hpp"
34 : #include "Utilities/ErrorHandling/Assert.hpp"
35 : #include "Utilities/Functional.hpp"
36 : #include "Utilities/GetOutput.hpp"
37 : #include "Utilities/Gsl.hpp"
38 : #include "Utilities/MakeWithValue.hpp"
39 : #include "Utilities/PrettyType.hpp"
40 : #include "Utilities/ProtocolHelpers.hpp"
41 : #include "Utilities/TMPL.hpp"
42 : #include "Utilities/TaggedTuple.hpp"
43 :
44 : /// \cond
45 : namespace NonlinearSolver::newton_raphson::detail {
46 : template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
47 : struct ResidualMonitor;
48 : template <typename FieldsTag, typename OptionsGroup, typename Label,
49 : typename ArraySectionIdTag>
50 : struct ReceiveInitialHasConverged;
51 : template <typename FieldsTag, typename OptionsGroup, typename Label,
52 : typename ArraySectionIdTag>
53 : struct PrepareStep;
54 : template <typename FieldsTag, typename OptionsGroup, typename Label,
55 : typename ArraySectionIdTag>
56 : struct Globalize;
57 : template <typename FieldsTag, typename OptionsGroup, typename Label,
58 : typename ArraySectionIdTag>
59 : struct CompleteStep;
60 : } // namespace NonlinearSolver::newton_raphson::detail
61 : /// \endcond
62 :
63 1 : namespace NonlinearSolver::newton_raphson::detail {
64 :
65 : using ResidualReductionData = Parallel::ReductionData<
66 : // Iteration ID
67 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
68 : // Globalization iteration ID
69 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
70 : // Residual magnitude square
71 : Parallel::ReductionDatum<double, funcl::Plus<>>,
72 : // Step length
73 : Parallel::ReductionDatum<double, funcl::AssertEqual<>>>;
74 :
75 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
76 : struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
77 : private:
78 : using fields_tag = FieldsTag;
79 : using source_tag = SourceTag;
80 : using nonlinear_operator_applied_to_fields_tag =
81 : db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, fields_tag>;
82 : using correction_tag =
83 : db::add_tag_prefix<NonlinearSolver::Tags::Correction, fields_tag>;
84 : using globalization_fields_tag =
85 : db::add_tag_prefix<NonlinearSolver::Tags::Globalization, fields_tag>;
86 :
87 : public: // Iterable action
88 : using simple_tags =
89 : tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
90 : Convergence::Tags::HasConverged<OptionsGroup>,
91 : nonlinear_operator_applied_to_fields_tag, correction_tag,
92 : NonlinearSolver::Tags::Globalization<
93 : Convergence::Tags::IterationId<OptionsGroup>>,
94 : NonlinearSolver::Tags::StepLength<OptionsGroup>,
95 : globalization_fields_tag>;
96 : using compute_tags = tmpl::list<
97 : NonlinearSolver::Tags::ResidualCompute<fields_tag, source_tag>>;
98 :
99 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
100 : typename ArrayIndex, typename ActionList,
101 : typename ParallelComponent>
102 : static Parallel::iterable_action_return_t apply(
103 : db::DataBox<DbTagsList>& box,
104 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
105 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
106 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
107 : const ParallelComponent* const /*meta*/) {
108 : ::Initialization::mutate_assign<
109 : tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
110 : NonlinearSolver::Tags::Globalization<
111 : Convergence::Tags::IterationId<OptionsGroup>>,
112 : NonlinearSolver::Tags::StepLength<OptionsGroup>>>(
113 : make_not_null(&box), std::numeric_limits<size_t>::max(),
114 : std::numeric_limits<size_t>::max(),
115 : std::numeric_limits<double>::signaling_NaN());
116 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
117 : }
118 :
119 : public: // amr::protocols::Projector
120 : using argument_tags = tmpl::list<>;
121 : using return_tags = simple_tags;
122 :
123 : template <typename... AmrData>
124 : static void apply(const gsl::not_null<size_t*> /*unused*/,
125 : const AmrData&... /*all_items*/) {
126 : // No need to reset or initialize any of the items during AMR because they
127 : // will be set in `PrepareSolve` below. AMR can't happen _during_ a solve.
128 : }
129 : };
130 :
131 : // Reset to the initial state of the algorithm.
132 : template <typename FieldsTag, typename OptionsGroup, typename Label,
133 : typename ArraySectionIdTag>
134 : struct PrepareSolve {
135 : private:
136 : using fields_tag = FieldsTag;
137 : using nonlinear_residual_tag =
138 : db::add_tag_prefix<NonlinearSolver::Tags::Residual, fields_tag>;
139 :
140 : public:
141 : using const_global_cache_tags =
142 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
143 :
144 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
145 : typename ArrayIndex, typename ActionList,
146 : typename ParallelComponent>
147 : static Parallel::iterable_action_return_t apply(
148 : db::DataBox<DbTagsList>& box,
149 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
150 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
151 : const ArrayIndex& array_index, const ActionList /*meta*/,
152 : const ParallelComponent* const /*meta*/) {
153 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>>(
154 : [](const gsl::not_null<size_t*> iteration_id) { *iteration_id = 0; },
155 : make_not_null(&box));
156 :
157 : // Skip the initial reduction on elements that are not part of the section
158 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
159 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
160 : ArraySectionIdTag>>(box)
161 : .has_value()) {
162 : constexpr size_t receive_initial_index = tmpl::index_of<
163 : ActionList,
164 : ReceiveInitialHasConverged<FieldsTag, OptionsGroup, Label,
165 : ArraySectionIdTag>>::value;
166 : return {Parallel::AlgorithmExecution::Continue, receive_initial_index};
167 : }
168 : }
169 :
170 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
171 : ::Verbosity::Debug)) {
172 : Parallel::printf("%s %s: Prepare solve\n", get_output(array_index),
173 : pretty_type::name<OptionsGroup>());
174 : }
175 :
176 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
177 : }
178 : };
179 :
180 : // To determine the initial residual magnitude and to check if the algorithm has
181 : // already converged, we perform a global reduction to the `ResidualMonitor`.
182 : template <typename FieldsTag, typename OptionsGroup, typename Label,
183 : typename ArraySectionIdTag>
184 : struct SendInitialResidualMagnitude {
185 : private:
186 : using fields_tag = FieldsTag;
187 : using nonlinear_residual_tag =
188 : db::add_tag_prefix<NonlinearSolver::Tags::Residual, fields_tag>;
189 :
190 : public:
191 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
192 : typename ArrayIndex, typename ActionList,
193 : typename ParallelComponent>
194 : static Parallel::iterable_action_return_t apply(
195 : db::DataBox<DbTagsList>& box,
196 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
197 : Parallel::GlobalCache<Metavariables>& cache,
198 : const ArrayIndex& array_index, const ActionList /*meta*/,
199 : const ParallelComponent* const /*meta*/) {
200 : // Perform a global reduction to compute the initial residual magnitude
201 : const auto& residual = db::get<nonlinear_residual_tag>(box);
202 : const double local_residual_magnitude_square =
203 : LinearSolver::inner_product(residual, residual);
204 : ResidualReductionData reduction_data{0, 0, local_residual_magnitude_square,
205 : 1.};
206 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
207 : make_not_null(&box));
208 : Parallel::contribute_to_reduction<
209 : CheckResidualMagnitude<FieldsTag, OptionsGroup, ParallelComponent>>(
210 : std::move(reduction_data),
211 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
212 : Parallel::get_parallel_component<
213 : ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache),
214 : make_not_null(§ion));
215 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
216 : }
217 : };
218 :
219 : // Wait for the broadcast from the `ResidualMonitor` to complete the preparation
220 : // for the solve. We skip the solve altogether if the algorithm has already
221 : // converged.
222 : template <typename FieldsTag, typename OptionsGroup, typename Label,
223 : typename ArraySectionIdTag>
224 : struct ReceiveInitialHasConverged {
225 : using inbox_tags = tmpl::list<Tags::GlobalizationResult<OptionsGroup>>;
226 :
227 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
228 : typename ArrayIndex, typename ActionList,
229 : typename ParallelComponent>
230 : static Parallel::iterable_action_return_t apply(
231 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
232 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
233 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
234 : const ParallelComponent* const /*meta*/) {
235 : const size_t iteration_id =
236 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
237 : auto& inbox = get<Tags::GlobalizationResult<OptionsGroup>>(inboxes);
238 : if (inbox.find(db::get<Convergence::Tags::IterationId<OptionsGroup>>(
239 : box)) == inbox.end()) {
240 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
241 : }
242 :
243 : // Retrieve reduction data from inbox
244 : auto globalization_result = std::move(inbox.extract(iteration_id).mapped());
245 : ASSERT(
246 : std::holds_alternative<Convergence::HasConverged>(globalization_result),
247 : "No globalization should occur for the initial residual. This is a "
248 : "bug, so please file an issue.");
249 : auto& has_converged = get<Convergence::HasConverged>(globalization_result);
250 : db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
251 : [&has_converged](const gsl::not_null<Convergence::HasConverged*>
252 : local_has_converged) {
253 : *local_has_converged = std::move(has_converged);
254 : },
255 : make_not_null(&box));
256 :
257 : // Skip steps entirely if the solve has already converged
258 : constexpr size_t complete_step_index =
259 : tmpl::index_of<ActionList, CompleteStep<FieldsTag, OptionsGroup, Label,
260 : ArraySectionIdTag>>::value;
261 : if (get<Convergence::Tags::HasConverged<OptionsGroup>>(box)) {
262 : return {Parallel::AlgorithmExecution::Continue, complete_step_index + 1};
263 : }
264 :
265 : // Skip the solve entirely on elements that are not part of the section. To
266 : // do so, we skip ahead to the `SolveLinearization` action list between
267 : // `PrepareStep` and `PerformStep`. Those actions need a chance to run even
268 : // on elements that are not part of the section, because they may take part
269 : // in preconditioning (see Multigrid preconditioner).
270 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
271 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
272 : ArraySectionIdTag>>(box)
273 : .has_value()) {
274 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>>(
275 : [](const gsl::not_null<size_t*> local_iteration_id) {
276 : *local_iteration_id = 1;
277 : },
278 : make_not_null(&box));
279 : constexpr size_t prepare_step_index =
280 : tmpl::index_of<ActionList,
281 : PrepareStep<FieldsTag, OptionsGroup, Label,
282 : ArraySectionIdTag>>::value;
283 : return {Parallel::AlgorithmExecution::Continue, prepare_step_index + 1};
284 : }
285 : }
286 :
287 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
288 : }
289 : };
290 :
291 : // Prepare the next Newton-Raphson step. In particular, we prepare the DataBox
292 : // for the linear solver which will run after this action to solve the
293 : // linearized problem for the `correction_tag`. The source for the linear
294 : // solver is the residual, which is in the DataBox as a compute tag so needs no
295 : // preparation. We only need to set the initial guess.
296 : //
297 : // We also prepare the line-search globalization here. Since we don't know if
298 : // a step will be sufficient before taking it, we have to store the original
299 : // field values.
300 : //
301 : // The algorithm jumps back to this action from `CompleteStep` to continue
302 : // iterating the nonlinear solve.
303 : template <typename FieldsTag, typename OptionsGroup, typename Label,
304 : typename ArraySectionIdTag>
305 : struct PrepareStep {
306 : private:
307 : using fields_tag = FieldsTag;
308 : using correction_tag =
309 : db::add_tag_prefix<NonlinearSolver::Tags::Correction, fields_tag>;
310 : using linear_operator_applied_to_correction_tag =
311 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, correction_tag>;
312 : using globalization_fields_tag =
313 : db::add_tag_prefix<NonlinearSolver::Tags::Globalization, fields_tag>;
314 :
315 : public:
316 : using const_global_cache_tags =
317 : tmpl::list<NonlinearSolver::Tags::DampingFactor<OptionsGroup>,
318 : logging::Tags::Verbosity<OptionsGroup>>;
319 :
320 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
321 : typename ArrayIndex, typename ActionList,
322 : typename ParallelComponent>
323 : static Parallel::iterable_action_return_t apply(
324 : db::DataBox<DbTagsList>& box,
325 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
326 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
327 : const ArrayIndex& array_index, const ActionList /*meta*/,
328 : const ParallelComponent* const /*meta*/) {
329 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
330 : ::Verbosity::Debug)) {
331 : Parallel::printf(
332 : "%s %s(%zu): Prepare step\n", get_output(array_index),
333 : pretty_type::name<OptionsGroup>(),
334 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box) + 1);
335 : }
336 :
337 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>, correction_tag,
338 : linear_operator_applied_to_correction_tag,
339 : NonlinearSolver::Tags::Globalization<
340 : Convergence::Tags::IterationId<OptionsGroup>>,
341 : NonlinearSolver::Tags::StepLength<OptionsGroup>,
342 : globalization_fields_tag>(
343 : [](const gsl::not_null<size_t*> iteration_id, const auto correction,
344 : const auto linear_operator_applied_to_correction,
345 : const gsl::not_null<size_t*> globalization_iteration_id,
346 : const gsl::not_null<double*> step_length,
347 : const auto globalization_fields, const auto& fields,
348 : const double damping_factor) {
349 : ++(*iteration_id);
350 : // Begin the linear solve with a zero initial guess
351 : *correction =
352 : make_with_value<typename correction_tag::type>(fields, 0.);
353 : // Since the initial guess is zero, we don't need to apply the linear
354 : // operator to it but can just set it to zero as well. Linear things
355 : // are nice :)
356 : *linear_operator_applied_to_correction = make_with_value<
357 : typename linear_operator_applied_to_correction_tag::type>(fields,
358 : 0.);
359 : // Prepare line search globalization
360 : *globalization_iteration_id = 0;
361 : *step_length = damping_factor;
362 : *globalization_fields = fields;
363 : },
364 : make_not_null(&box), db::get<fields_tag>(box),
365 : db::get<NonlinearSolver::Tags::DampingFactor<OptionsGroup>>(box));
366 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
367 : }
368 : };
369 :
370 : // Between `PrepareStep` and this action the linear solver has run, so the
371 : // `correction_tag` is now filled with the solution to the linearized problem.
372 : // We just take a step in the direction of the correction.
373 : //
374 : // The `Globalize` action will jump back here in case the step turned out to not
375 : // decrease the residual sufficiently. It will have adjusted the step length, so
376 : // we can just try again with a step of that length.
377 : template <typename FieldsTag, typename OptionsGroup, typename Label,
378 : typename ArraySectionIdTag>
379 : struct PerformStep {
380 : private:
381 : using fields_tag = FieldsTag;
382 : using correction_tag =
383 : db::add_tag_prefix<NonlinearSolver::Tags::Correction, fields_tag>;
384 : using globalization_fields_tag =
385 : db::add_tag_prefix<NonlinearSolver::Tags::Globalization, fields_tag>;
386 :
387 : public:
388 : using const_global_cache_tags =
389 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
390 :
391 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
392 : typename ArrayIndex, typename ActionList,
393 : typename ParallelComponent>
394 : static Parallel::iterable_action_return_t apply(
395 : db::DataBox<DbTagsList>& box,
396 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
397 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
398 : const ArrayIndex& array_index, const ActionList /*meta*/,
399 : const ParallelComponent* const /*meta*/) {
400 : // Skip to the end of the step on elements that are not part of the section
401 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
402 : if (not db::get<
403 : Parallel::Tags::Section<ParallelComponent, ArraySectionIdTag>>(
404 : box)) {
405 : constexpr size_t globalize_index =
406 : tmpl::index_of<ActionList, Globalize<FieldsTag, OptionsGroup, Label,
407 : ArraySectionIdTag>>::value;
408 : return {Parallel::AlgorithmExecution::Continue, globalize_index};
409 : }
410 : }
411 :
412 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
413 : ::Verbosity::Debug)) {
414 : Parallel::printf(
415 : "%s %s(%zu): Perform step with length: %g\n", get_output(array_index),
416 : pretty_type::name<OptionsGroup>(),
417 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box),
418 : db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box));
419 : }
420 :
421 : // Apply the correction that the linear solve has determined to attempt
422 : // improving the nonlinear solution
423 : db::mutate<fields_tag>(
424 : [](const auto fields, const auto& correction, const double step_length,
425 : const auto& globalization_fields) {
426 : *fields = globalization_fields + step_length * correction;
427 : },
428 : make_not_null(&box), db::get<correction_tag>(box),
429 : db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box),
430 : db::get<globalization_fields_tag>(box));
431 :
432 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
433 : }
434 : };
435 :
436 : // Between `PerformStep` and this action the nonlinear operator has been applied
437 : // to the updated fields. The residual is up-to-date because it is a compute
438 : // tag, so at this point we need to check if the step has sufficiently reduced
439 : // the residual. We perform a global reduction to the `ResidualMonitor` for this
440 : // purpose.
441 : template <typename FieldsTag, typename OptionsGroup, typename Label,
442 : typename ArraySectionIdTag>
443 : struct ContributeToResidualMagnitudeReduction {
444 : private:
445 : using fields_tag = FieldsTag;
446 : using nonlinear_residual_tag =
447 : db::add_tag_prefix<NonlinearSolver::Tags::Residual, fields_tag>;
448 :
449 : public:
450 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
451 : typename ArrayIndex, typename ActionList,
452 : typename ParallelComponent>
453 : static Parallel::iterable_action_return_t apply(
454 : db::DataBox<DbTagsList>& box,
455 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
456 : Parallel::GlobalCache<Metavariables>& cache,
457 : const ArrayIndex& array_index, const ActionList /*meta*/,
458 : const ParallelComponent* const /*meta*/) {
459 : const auto& residual = db::get<nonlinear_residual_tag>(box);
460 : const double local_residual_magnitude_square =
461 : LinearSolver::inner_product(residual, residual);
462 : ResidualReductionData reduction_data{
463 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box),
464 : db::get<NonlinearSolver::Tags::Globalization<
465 : Convergence::Tags::IterationId<OptionsGroup>>>(box),
466 : local_residual_magnitude_square,
467 : db::get<NonlinearSolver::Tags::StepLength<OptionsGroup>>(box)};
468 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
469 : make_not_null(&box));
470 : Parallel::contribute_to_reduction<
471 : CheckResidualMagnitude<FieldsTag, OptionsGroup, ParallelComponent>>(
472 : std::move(reduction_data),
473 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
474 : Parallel::get_parallel_component<
475 : ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache),
476 : make_not_null(§ion));
477 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
478 : }
479 : };
480 :
481 : // Wait for the `ResidualMonitor` to broadcast whether or not it has determined
482 : // the step has sufficiently decreased the residual. If it has, we just proceed
483 : // to `CompleteStep`. If it hasn't, the `ResidualMonitor` has also sent the new
484 : // step length along, so we mutate it in the DataBox and jump back to
485 : // `PerformStep` to try again with the updated step length.
486 : template <typename FieldsTag, typename OptionsGroup, typename Label,
487 : typename ArraySectionIdTag>
488 : struct Globalize {
489 : using const_global_cache_tags =
490 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
491 : using inbox_tags = tmpl::list<Tags::GlobalizationResult<OptionsGroup>>;
492 :
493 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
494 : typename ArrayIndex, typename ActionList,
495 : typename ParallelComponent>
496 : static Parallel::iterable_action_return_t apply(
497 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
498 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
499 : const ArrayIndex& array_index, const ActionList /*meta*/,
500 : const ParallelComponent* const /*meta*/) {
501 : auto& inbox = get<Tags::GlobalizationResult<OptionsGroup>>(inboxes);
502 : const size_t iteration_id =
503 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
504 : if (inbox.find(iteration_id) == inbox.end()) {
505 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
506 : }
507 :
508 : // Retrieve reduction data from inbox
509 : auto globalization_result = std::move(inbox.extract(iteration_id).mapped());
510 :
511 : // Elements that are not part of the section jump directly to the
512 : // `SolveLinearization` for the next step.
513 : constexpr size_t complete_step_index =
514 : tmpl::index_of<ActionList, CompleteStep<FieldsTag, OptionsGroup, Label,
515 : ArraySectionIdTag>>::value;
516 : constexpr size_t prepare_step_index =
517 : tmpl::index_of<ActionList, PrepareStep<FieldsTag, OptionsGroup, Label,
518 : ArraySectionIdTag>>::value;
519 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
520 : if (not db::get<
521 : Parallel::Tags::Section<ParallelComponent, ArraySectionIdTag>>(
522 : box)) {
523 : const bool globalization_is_complete =
524 : std::holds_alternative<Convergence::HasConverged>(
525 : globalization_result);
526 : // Wait until globalization is complete
527 : if (not globalization_is_complete) {
528 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
529 : }
530 : auto& has_converged =
531 : get<Convergence::HasConverged>(globalization_result);
532 :
533 : db::mutate<Convergence::Tags::HasConverged<OptionsGroup>,
534 : Convergence::Tags::IterationId<OptionsGroup>>(
535 : [&has_converged](const gsl::not_null<Convergence::HasConverged*>
536 : local_has_converged,
537 : const gsl::not_null<size_t*> local_iteration_id) {
538 : *local_has_converged = std::move(has_converged);
539 : ++(*local_iteration_id);
540 : },
541 : make_not_null(&box));
542 :
543 : return {Parallel::AlgorithmExecution::Continue,
544 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
545 : ? (complete_step_index + 1)
546 : : (prepare_step_index + 1)};
547 : }
548 : }
549 :
550 : if (std::holds_alternative<double>(globalization_result)) {
551 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
552 : ::Verbosity::Debug)) {
553 : Parallel::printf(
554 : "%s %s(%zu): Globalize(%zu)\n", get_output(array_index),
555 : pretty_type::name<OptionsGroup>(),
556 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box),
557 : db::get<NonlinearSolver::Tags::Globalization<
558 : Convergence::Tags::IterationId<OptionsGroup>>>(box));
559 : }
560 :
561 : // Update the step length
562 : db::mutate<NonlinearSolver::Tags::StepLength<OptionsGroup>,
563 : NonlinearSolver::Tags::Globalization<
564 : Convergence::Tags::IterationId<OptionsGroup>>>(
565 : [&globalization_result](
566 : const gsl::not_null<double*> step_length,
567 : const gsl::not_null<size_t*> globalization_iteration_id) {
568 : *step_length = get<double>(globalization_result);
569 : ++(*globalization_iteration_id);
570 : },
571 : make_not_null(&box));
572 : // Continue globalization by taking a step with the updated step length
573 : // and checking the residual again
574 : constexpr size_t perform_step_index =
575 : tmpl::index_of<ActionList, PerformStep<FieldsTag, OptionsGroup, Label,
576 : ArraySectionIdTag>>::value;
577 : return {Parallel::AlgorithmExecution::Continue, perform_step_index};
578 : }
579 :
580 : // At this point globalization is complete, so we proceed with the algorithm
581 : auto& has_converged = get<Convergence::HasConverged>(globalization_result);
582 :
583 : db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
584 : [&has_converged](const gsl::not_null<Convergence::HasConverged*>
585 : local_has_converged) {
586 : *local_has_converged = std::move(has_converged);
587 : },
588 : make_not_null(&box));
589 :
590 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
591 : }
592 : };
593 :
594 : // Jump back to `PrepareStep` to continue iterating if the algorithm has not yet
595 : // converged, or complete the solve and proceed with the action list if it has
596 : // converged. This is a separate action because the user has the opportunity to
597 : // insert actions before the step completes, for example to do observations.
598 : template <typename FieldsTag, typename OptionsGroup, typename Label,
599 : typename ArraySectionIdTag>
600 : struct CompleteStep {
601 : using const_global_cache_tags =
602 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
603 :
604 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
605 : typename ArrayIndex, typename ActionList,
606 : typename ParallelComponent>
607 : static Parallel::iterable_action_return_t apply(
608 : db::DataBox<DbTagsList>& box,
609 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
610 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
611 : const ArrayIndex& array_index, const ActionList /*meta*/,
612 : const ParallelComponent* const /*meta*/) {
613 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
614 : ::Verbosity::Debug)) {
615 : Parallel::printf(
616 : "%s %s(%zu): Complete step\n", get_output(array_index),
617 : pretty_type::name<OptionsGroup>(),
618 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box));
619 : }
620 :
621 : // Repeat steps until the solve has converged
622 : constexpr size_t prepare_step_index =
623 : tmpl::index_of<ActionList, PrepareStep<FieldsTag, OptionsGroup, Label,
624 : ArraySectionIdTag>>::value;
625 : constexpr size_t this_action_index =
626 : tmpl::index_of<ActionList, CompleteStep>::value;
627 : return {Parallel::AlgorithmExecution::Continue,
628 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
629 : ? (this_action_index + 1)
630 : : prepare_step_index};
631 : }
632 : };
633 :
634 : } // namespace NonlinearSolver::newton_raphson::detail
|
