Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 : #include <optional>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
16 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
18 : #include "Parallel/AlgorithmExecution.hpp"
19 : #include "Parallel/GetSection.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "Parallel/Invoke.hpp"
22 : #include "Parallel/Printf/Printf.hpp"
23 : #include "Parallel/Reduction.hpp"
24 : #include "Parallel/Tags/Section.hpp"
25 : #include "ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp"
26 : #include "ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp"
27 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
28 : #include "Utilities/ErrorHandling/Assert.hpp"
29 : #include "Utilities/Functional.hpp"
30 : #include "Utilities/GetOutput.hpp"
31 : #include "Utilities/Gsl.hpp"
32 : #include "Utilities/MakeWithValue.hpp"
33 : #include "Utilities/PrettyType.hpp"
34 : #include "Utilities/Requires.hpp"
35 : #include "Utilities/TMPL.hpp"
36 :
37 : /// \cond
38 : namespace tuples {
39 : template <typename...>
40 : class TaggedTuple;
41 : } // namespace tuples
42 : namespace LinearSolver::gmres::detail {
43 : template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
44 : struct ResidualMonitor;
45 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
46 : typename Label, typename ArraySectionIdTag>
47 : struct PrepareStep;
48 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
49 : typename Label, typename ArraySectionIdTag>
50 : struct NormalizeOperandAndUpdateField;
51 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
52 : typename Label, typename ArraySectionIdTag>
53 : struct CompleteStep;
54 : } // namespace LinearSolver::gmres::detail
55 : /// \endcond
56 :
57 : namespace LinearSolver::gmres::detail {
58 :
59 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
60 : typename Label, typename SourceTag, typename ArraySectionIdTag>
61 : struct PrepareSolve {
62 : private:
63 : using fields_tag = FieldsTag;
64 : using initial_fields_tag = db::add_tag_prefix<::Tags::Initial, fields_tag>;
65 : using source_tag = SourceTag;
66 : using operator_applied_to_fields_tag =
67 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, fields_tag>;
68 : using operand_tag =
69 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
70 : using basis_history_tag =
71 : LinearSolver::Tags::KrylovSubspaceBasis<operand_tag>;
72 :
73 : public:
74 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
75 : typename ArrayIndex, typename ActionList,
76 : typename ParallelComponent>
77 : static Parallel::iterable_action_return_t apply(
78 : db::DataBox<DbTagsList>& box,
79 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
80 : Parallel::GlobalCache<Metavariables>& cache,
81 : const ArrayIndex& array_index, const ActionList /*meta*/,
82 : const ParallelComponent* const /*meta*/) {
83 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>>(
84 : [](const gsl::not_null<size_t*> iteration_id) { *iteration_id = 0; },
85 : make_not_null(&box));
86 :
87 : // Skip the initial reduction on elements that are not part of the section
88 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
89 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
90 : ArraySectionIdTag>>(box)
91 : .has_value()) {
92 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
93 : }
94 : }
95 :
96 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
97 : ::Verbosity::Debug)) {
98 : Parallel::printf("%s %s: Prepare solve\n", get_output(array_index),
99 : pretty_type::name<OptionsGroup>());
100 : }
101 :
102 : db::mutate<operand_tag, initial_fields_tag, basis_history_tag>(
103 : [](const auto operand, const auto initial_fields,
104 : const auto basis_history, const auto& source,
105 : const auto& operator_applied_to_fields, const auto& fields) {
106 : *operand = source - operator_applied_to_fields;
107 : *initial_fields = fields;
108 : *basis_history = typename basis_history_tag::type{};
109 : },
110 : make_not_null(&box), get<source_tag>(box),
111 : get<operator_applied_to_fields_tag>(box), get<fields_tag>(box));
112 :
113 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
114 : make_not_null(&box));
115 : Parallel::contribute_to_reduction<InitializeResidualMagnitude<
116 : FieldsTag, OptionsGroup, ParallelComponent>>(
117 : Parallel::ReductionData<
118 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Sqrt<>>>{
119 : inner_product(get<operand_tag>(box), get<operand_tag>(box))},
120 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
121 : Parallel::get_parallel_component<
122 : ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache),
123 : make_not_null(§ion));
124 :
125 : if constexpr (Preconditioned) {
126 : using preconditioned_operand_tag =
127 : db::add_tag_prefix<LinearSolver::Tags::Preconditioned, operand_tag>;
128 : using preconditioned_basis_history_tag =
129 : LinearSolver::Tags::KrylovSubspaceBasis<preconditioned_operand_tag>;
130 :
131 : db::mutate<preconditioned_basis_history_tag>(
132 : [](const auto preconditioned_basis_history) {
133 : *preconditioned_basis_history =
134 : typename preconditioned_basis_history_tag::type{};
135 : },
136 : make_not_null(&box));
137 : }
138 :
139 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
140 : }
141 : };
142 :
143 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
144 : typename Label, typename ArraySectionIdTag>
145 : struct NormalizeInitialOperand {
146 : private:
147 : using fields_tag = FieldsTag;
148 : using operand_tag =
149 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
150 : using basis_history_tag =
151 : LinearSolver::Tags::KrylovSubspaceBasis<operand_tag>;
152 :
153 : public:
154 : using const_global_cache_tags =
155 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
156 : using inbox_tags = tmpl::list<Tags::InitialOrthogonalization<OptionsGroup>>;
157 :
158 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
159 : typename ArrayIndex, typename ActionList,
160 : typename ParallelComponent>
161 : static Parallel::iterable_action_return_t apply(
162 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
163 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
164 : const ArrayIndex& array_index, const ActionList /*meta*/,
165 : const ParallelComponent* const /*meta*/) {
166 : const size_t iteration_id =
167 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
168 : auto& inbox = get<Tags::InitialOrthogonalization<OptionsGroup>>(inboxes);
169 : if (inbox.find(iteration_id) == inbox.end()) {
170 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
171 : }
172 :
173 : auto received_data = std::move(inbox.extract(iteration_id).mapped());
174 : const double residual_magnitude = get<0>(received_data);
175 : auto& has_converged = get<1>(received_data);
176 : db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
177 : [&has_converged](const gsl::not_null<Convergence::HasConverged*>
178 : local_has_converged) {
179 : *local_has_converged = std::move(has_converged);
180 : },
181 : make_not_null(&box));
182 :
183 : // Skip steps entirely if the solve has already converged
184 : constexpr size_t step_end_index =
185 : tmpl::index_of<ActionList,
186 : CompleteStep<FieldsTag, OptionsGroup, Preconditioned,
187 : Label, ArraySectionIdTag>>::value;
188 : if (get<Convergence::Tags::HasConverged<OptionsGroup>>(box)) {
189 : return {Parallel::AlgorithmExecution::Continue, step_end_index + 1};
190 : }
191 :
192 : // Skip the solve entirely on elements that are not part of the section. To
193 : // do so, we skip ahead to the `ApplyOperatorActions` action list between
194 : // `PrepareStep` and `PerformStep`. Those actions need a chance to run even
195 : // on elements that are not part of the section, because they may take part
196 : // in preconditioning (see Multigrid preconditioner).
197 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
198 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
199 : ArraySectionIdTag>>(box)
200 : .has_value()) {
201 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
202 : }
203 : }
204 :
205 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
206 : ::Verbosity::Debug)) {
207 : Parallel::printf("%s %s(%zu): Normalize initial operand\n",
208 : get_output(array_index),
209 : pretty_type::name<OptionsGroup>(), iteration_id);
210 : }
211 :
212 : db::mutate<operand_tag, basis_history_tag>(
213 : [residual_magnitude](const auto operand, const auto basis_history) {
214 : *operand /= residual_magnitude;
215 : basis_history->push_back(*operand);
216 : },
217 : make_not_null(&box));
218 :
219 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
220 : }
221 : };
222 :
223 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
224 : typename Label, typename ArraySectionIdTag>
225 : struct PrepareStep {
226 : using const_global_cache_tags =
227 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
228 :
229 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
230 : typename ArrayIndex, typename ActionList,
231 : typename ParallelComponent>
232 : static Parallel::iterable_action_return_t apply(
233 : db::DataBox<DbTagsList>& box,
234 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
235 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
236 : const ArrayIndex& array_index, const ActionList /*meta*/,
237 : const ParallelComponent* const /*meta*/) {
238 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>>(
239 : [](const gsl::not_null<size_t*> iteration_id) { ++(*iteration_id); },
240 : make_not_null(&box));
241 :
242 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
243 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
244 : ArraySectionIdTag>>(box)
245 : .has_value()) {
246 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
247 : }
248 : }
249 :
250 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
251 : ::Verbosity::Debug)) {
252 : Parallel::printf(
253 : "%s %s(%zu): Prepare step\n", get_output(array_index),
254 : pretty_type::name<OptionsGroup>(),
255 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box));
256 : }
257 :
258 : if constexpr (Preconditioned) {
259 : using fields_tag = FieldsTag;
260 : using operand_tag =
261 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
262 : using preconditioned_operand_tag =
263 : db::add_tag_prefix<LinearSolver::Tags::Preconditioned, operand_tag>;
264 : using operator_tag = db::add_tag_prefix<
265 : LinearSolver::Tags::OperatorAppliedTo,
266 : std::conditional_t<Preconditioned, preconditioned_operand_tag,
267 : operand_tag>>;
268 :
269 : db::mutate<preconditioned_operand_tag, operator_tag>(
270 : [](const auto preconditioned_operand,
271 : const auto operator_applied_to_operand, const auto& operand) {
272 : // Start the preconditioner at zero because we have no reason to
273 : // expect the remaining residual to have a particular form.
274 : // Another possibility would be to start the preconditioner with an
275 : // initial guess equal to its source, so not running the
276 : // preconditioner at all means it is the identity, but that approach
277 : // appears to yield worse results.
278 : *preconditioned_operand =
279 : make_with_value<typename preconditioned_operand_tag::type>(
280 : operand, 0.);
281 : // Also set the operator applied to the initial preconditioned
282 : // operand to zero because it's linear. This may save the
283 : // preconditioner an operator application if it's optimized for
284 : // this.
285 : *operator_applied_to_operand =
286 : make_with_value<typename operator_tag::type>(operand, 0.);
287 : },
288 : make_not_null(&box), get<operand_tag>(box));
289 : }
290 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
291 : }
292 : };
293 :
294 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
295 : typename Label, typename ArraySectionIdTag>
296 : struct PerformStep {
297 : private:
298 : using fields_tag = FieldsTag;
299 : using operand_tag =
300 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
301 : using preconditioned_operand_tag =
302 : db::add_tag_prefix<LinearSolver::Tags::Preconditioned, operand_tag>;
303 :
304 : public:
305 : using const_global_cache_tags =
306 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
307 :
308 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
309 : typename ArrayIndex, typename ActionList,
310 : typename ParallelComponent>
311 : static Parallel::iterable_action_return_t apply(
312 : db::DataBox<DbTagsList>& box,
313 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
314 : Parallel::GlobalCache<Metavariables>& cache,
315 : const ArrayIndex& array_index, const ActionList /*meta*/,
316 : const ParallelComponent* const /*meta*/) {
317 : // Skip to the end of the step on elements that are not part of the section
318 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
319 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
320 : ArraySectionIdTag>>(box)
321 : .has_value()) {
322 : constexpr size_t step_end_index =
323 : tmpl::index_of<ActionList,
324 : NormalizeOperandAndUpdateField<
325 : FieldsTag, OptionsGroup, Preconditioned, Label,
326 : ArraySectionIdTag>>::value;
327 : return {Parallel::AlgorithmExecution::Continue, step_end_index};
328 : }
329 : }
330 :
331 : const size_t iteration_id =
332 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
333 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
334 : ::Verbosity::Debug)) {
335 : Parallel::printf("%s %s(%zu): Perform step\n", get_output(array_index),
336 : pretty_type::name<OptionsGroup>(), iteration_id);
337 : }
338 :
339 : using operator_tag = db::add_tag_prefix<
340 : LinearSolver::Tags::OperatorAppliedTo,
341 : std::conditional_t<Preconditioned, preconditioned_operand_tag,
342 : operand_tag>>;
343 : using orthogonalization_iteration_id_tag =
344 : LinearSolver::Tags::Orthogonalization<
345 : Convergence::Tags::IterationId<OptionsGroup>>;
346 : using basis_history_tag =
347 : LinearSolver::Tags::KrylovSubspaceBasis<operand_tag>;
348 :
349 : if constexpr (Preconditioned) {
350 : using preconditioned_basis_history_tag =
351 : LinearSolver::Tags::KrylovSubspaceBasis<preconditioned_operand_tag>;
352 :
353 : db::mutate<preconditioned_basis_history_tag>(
354 : [](const auto preconditioned_basis_history,
355 : const auto& preconditioned_operand) {
356 : preconditioned_basis_history->push_back(preconditioned_operand);
357 : },
358 : make_not_null(&box), get<preconditioned_operand_tag>(box));
359 : }
360 :
361 : db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
362 : [](const auto operand,
363 : const gsl::not_null<size_t*> orthogonalization_iteration_id,
364 : const auto& operator_action) {
365 : *operand = typename operand_tag::type(operator_action);
366 : *orthogonalization_iteration_id = 0;
367 : },
368 : make_not_null(&box), get<operator_tag>(box));
369 :
370 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
371 : make_not_null(&box));
372 : Parallel::contribute_to_reduction<
373 : StoreOrthogonalization<FieldsTag, OptionsGroup, ParallelComponent>>(
374 : Parallel::ReductionData<
375 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
376 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
377 : Parallel::ReductionDatum<double, funcl::Plus<>>>{
378 : get<Convergence::Tags::IterationId<OptionsGroup>>(box),
379 : get<orthogonalization_iteration_id_tag>(box),
380 : inner_product(get<basis_history_tag>(box)[0],
381 : get<operand_tag>(box))},
382 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
383 : Parallel::get_parallel_component<
384 : ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache),
385 : make_not_null(§ion));
386 :
387 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
388 : }
389 : };
390 :
391 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
392 : typename Label, typename ArraySectionIdTag>
393 : struct OrthogonalizeOperand {
394 : private:
395 : using fields_tag = FieldsTag;
396 : using operand_tag =
397 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
398 : using orthogonalization_iteration_id_tag =
399 : LinearSolver::Tags::Orthogonalization<
400 : Convergence::Tags::IterationId<OptionsGroup>>;
401 : using basis_history_tag =
402 : LinearSolver::Tags::KrylovSubspaceBasis<operand_tag>;
403 :
404 : public:
405 : using inbox_tags = tmpl::list<Tags::Orthogonalization<OptionsGroup>>;
406 :
407 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
408 : typename ArrayIndex, typename ActionList,
409 : typename ParallelComponent>
410 : static Parallel::iterable_action_return_t apply(
411 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
412 : Parallel::GlobalCache<Metavariables>& cache,
413 : const ArrayIndex& array_index, const ActionList /*meta*/,
414 : const ParallelComponent* const /*meta*/) {
415 : const size_t iteration_id =
416 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
417 : auto& inbox = get<Tags::Orthogonalization<OptionsGroup>>(inboxes);
418 : if (inbox.find(iteration_id) == inbox.end()) {
419 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
420 : }
421 :
422 : const double orthogonalization =
423 : std::move(inbox.extract(iteration_id).mapped());
424 :
425 : db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
426 : [orthogonalization](
427 : const auto operand,
428 : const gsl::not_null<size_t*> orthogonalization_iteration_id,
429 : const auto& basis_history) {
430 : *operand -= orthogonalization *
431 : gsl::at(basis_history, *orthogonalization_iteration_id);
432 : ++(*orthogonalization_iteration_id);
433 : },
434 : make_not_null(&box), get<basis_history_tag>(box));
435 :
436 : const auto& next_orthogonalization_iteration_id =
437 : get<orthogonalization_iteration_id_tag>(box);
438 : const bool orthogonalization_complete =
439 : next_orthogonalization_iteration_id == iteration_id;
440 : const double local_orthogonalization =
441 : inner_product(orthogonalization_complete
442 : ? get<operand_tag>(box)
443 : : gsl::at(get<basis_history_tag>(box),
444 : next_orthogonalization_iteration_id),
445 : get<operand_tag>(box));
446 :
447 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
448 : make_not_null(&box));
449 : Parallel::contribute_to_reduction<
450 : StoreOrthogonalization<FieldsTag, OptionsGroup, ParallelComponent>>(
451 : Parallel::ReductionData<
452 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
453 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
454 : Parallel::ReductionDatum<double, funcl::Plus<>>>{
455 : iteration_id, next_orthogonalization_iteration_id,
456 : local_orthogonalization},
457 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
458 : Parallel::get_parallel_component<
459 : ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache),
460 : make_not_null(§ion));
461 :
462 : // Repeat this action until orthogonalization is complete
463 : constexpr size_t this_action_index =
464 : tmpl::index_of<ActionList, OrthogonalizeOperand>::value;
465 : return {Parallel::AlgorithmExecution::Continue,
466 : orthogonalization_complete ? (this_action_index + 1)
467 : : this_action_index};
468 : }
469 : };
470 :
471 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
472 : typename Label, typename ArraySectionIdTag>
473 : struct NormalizeOperandAndUpdateField {
474 : private:
475 : using fields_tag = FieldsTag;
476 : using initial_fields_tag = db::add_tag_prefix<::Tags::Initial, fields_tag>;
477 : using operand_tag =
478 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
479 : using preconditioned_operand_tag =
480 : db::add_tag_prefix<LinearSolver::Tags::Preconditioned, operand_tag>;
481 : using basis_history_tag =
482 : LinearSolver::Tags::KrylovSubspaceBasis<operand_tag>;
483 : using preconditioned_basis_history_tag =
484 : LinearSolver::Tags::KrylovSubspaceBasis<std::conditional_t<
485 : Preconditioned, preconditioned_operand_tag, operand_tag>>;
486 :
487 : public:
488 : using const_global_cache_tags =
489 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
490 : using inbox_tags = tmpl::list<Tags::FinalOrthogonalization<OptionsGroup>>;
491 :
492 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
493 : typename ArrayIndex, typename ActionList,
494 : typename ParallelComponent>
495 : static Parallel::iterable_action_return_t apply(
496 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
497 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
498 : const ArrayIndex& array_index, const ActionList /*meta*/,
499 : const ParallelComponent* const /*meta*/) {
500 : const size_t iteration_id =
501 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
502 : auto& inbox = get<Tags::FinalOrthogonalization<OptionsGroup>>(inboxes);
503 : if (inbox.find(iteration_id) == inbox.end()) {
504 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
505 : }
506 :
507 : // Retrieve reduction data from inbox
508 : auto received_data = std::move(inbox.extract(iteration_id).mapped());
509 : const double normalization = get<0>(received_data);
510 : const auto& minres = get<1>(received_data);
511 : db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
512 : [&received_data](
513 : const gsl::not_null<Convergence::HasConverged*> has_converged) {
514 : *has_converged = std::move(get<2>(received_data));
515 : },
516 : make_not_null(&box));
517 :
518 : // Elements that are not part of the section jump directly to the
519 : // `ApplyOperationActions` for the next step.
520 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
521 : constexpr size_t complete_step_index =
522 : tmpl::index_of<ActionList,
523 : CompleteStep<FieldsTag, OptionsGroup, Preconditioned,
524 : Label, ArraySectionIdTag>>::value;
525 : constexpr size_t prepare_step_index =
526 : tmpl::index_of<ActionList,
527 : PrepareStep<FieldsTag, OptionsGroup, Preconditioned,
528 : Label, ArraySectionIdTag>>::value;
529 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
530 : ArraySectionIdTag>>(box)
531 : .has_value()) {
532 : return {Parallel::AlgorithmExecution::Continue,
533 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
534 : ? (complete_step_index + 1)
535 : : prepare_step_index};
536 : }
537 : }
538 :
539 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
540 : ::Verbosity::Debug)) {
541 : Parallel::printf("%s %s(%zu): Update field\n", get_output(array_index),
542 : pretty_type::name<OptionsGroup>(), iteration_id);
543 : }
544 :
545 : db::mutate<operand_tag, basis_history_tag, fields_tag>(
546 : [normalization, &minres](const auto operand, const auto basis_history,
547 : const auto field, const auto& initial_field,
548 : const auto& preconditioned_basis_history,
549 : const auto& has_converged) {
550 : // Avoid an FPE if the new operand norm is exactly zero. In that case
551 : // the problem is solved and the algorithm will terminate (see
552 : // Proposition 9.3 in \cite Saad2003). Since there will be no next
553 : // iteration we don't need to normalize the operand.
554 : if (LIKELY(normalization > 0.)) {
555 : *operand /= normalization;
556 : }
557 : basis_history->push_back(*operand);
558 : // Don't update the solution if an error occurred
559 : if (not(has_converged and
560 : has_converged.reason() == Convergence::Reason::Error)) {
561 : *field = initial_field;
562 : for (size_t i = 0; i < minres.size(); i++) {
563 : *field += minres[i] * gsl::at(preconditioned_basis_history, i);
564 : }
565 : }
566 : },
567 : make_not_null(&box), get<initial_fields_tag>(box),
568 : get<preconditioned_basis_history_tag>(box),
569 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box));
570 :
571 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
572 : }
573 : };
574 :
575 : // Jump back to `PrepareStep` to continue iterating if the algorithm has not yet
576 : // converged, or complete the solve and proceed with the action list if it has
577 : // converged. This is a separate action because the user has the opportunity to
578 : // insert actions before the step completes, for example to do observations.
579 : template <typename FieldsTag, typename OptionsGroup, bool Preconditioned,
580 : typename Label, typename ArraySectionIdTag>
581 : struct CompleteStep {
582 : using const_global_cache_tags =
583 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
584 :
585 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
586 : typename ArrayIndex, typename ActionList,
587 : typename ParallelComponent>
588 : static Parallel::iterable_action_return_t apply(
589 : db::DataBox<DbTagsList>& box,
590 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
591 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
592 : const ArrayIndex& array_index, const ActionList /*meta*/,
593 : const ParallelComponent* const /*meta*/) {
594 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
595 : ::Verbosity::Debug)) {
596 : Parallel::printf(
597 : "%s %s(%zu): Complete step\n", get_output(array_index),
598 : pretty_type::name<OptionsGroup>(),
599 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box));
600 : }
601 :
602 : // Repeat steps until the solve has converged
603 : constexpr size_t prepare_step_index =
604 : tmpl::index_of<ActionList,
605 : PrepareStep<FieldsTag, OptionsGroup, Preconditioned,
606 : Label, ArraySectionIdTag>>::value;
607 : constexpr size_t this_action_index =
608 : tmpl::index_of<ActionList, CompleteStep>::value;
609 : return {Parallel::AlgorithmExecution::Continue,
610 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
611 : ? (this_action_index + 1)
612 : : prepare_step_index};
613 : }
614 : };
615 :
616 : } // namespace LinearSolver::gmres::detail
|
