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 <string>
10 : #include <tuple>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 : #include "Domain/Structure/Element.hpp"
17 : #include "IO/Logging/Tags.hpp"
18 : #include "IO/Logging/Verbosity.hpp"
19 : #include "IO/Observer/Actions/RegisterWithObservers.hpp"
20 : #include "IO/Observer/GetSectionObservationKey.hpp"
21 : #include "IO/Observer/Helpers.hpp"
22 : #include "IO/Observer/ObservationId.hpp"
23 : #include "IO/Observer/ObserverComponent.hpp"
24 : #include "IO/Observer/Protocols/ReductionDataFormatter.hpp"
25 : #include "IO/Observer/ReductionActions.hpp"
26 : #include "IO/Observer/Tags.hpp"
27 : #include "IO/Observer/TypeOfObservation.hpp"
28 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
29 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
30 : #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
31 : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
32 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
33 : #include "Parallel/AlgorithmExecution.hpp"
34 : #include "Parallel/ArrayComponentId.hpp"
35 : #include "Parallel/GlobalCache.hpp"
36 : #include "Parallel/Invoke.hpp"
37 : #include "Parallel/Local.hpp"
38 : #include "Parallel/Printf/Printf.hpp"
39 : #include "Parallel/Reduction.hpp"
40 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
41 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
42 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
43 : #include "Utilities/Functional.hpp"
44 : #include "Utilities/GetOutput.hpp"
45 : #include "Utilities/Gsl.hpp"
46 : #include "Utilities/PrettyType.hpp"
47 : #include "Utilities/ProtocolHelpers.hpp"
48 : #include "Utilities/Requires.hpp"
49 : #include "Utilities/TMPL.hpp"
50 :
51 : /// \cond
52 : namespace tuples {
53 : template <typename...>
54 : class TaggedTuple;
55 : } // namespace tuples
56 : namespace LinearSolver::async_solvers {
57 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
58 : typename Label, typename ArraySectionIdTag,
59 : bool ObserveInitialResidual>
60 : struct CompleteStep;
61 : } // namespace LinearSolver::async_solvers
62 : /// \endcond
63 :
64 : /// Functionality shared between parallel linear solvers that have no global
65 : /// synchronization points
66 1 : namespace LinearSolver::async_solvers {
67 :
68 0 : using reduction_data = Parallel::ReductionData<
69 : // Iteration
70 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
71 : // Residual
72 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Sqrt<>>>;
73 :
74 : template <typename OptionsGroup>
75 0 : struct ResidualReductionFormatter
76 : : tt::ConformsTo<observers::protocols::ReductionDataFormatter> {
77 0 : using reduction_data = async_solvers::reduction_data;
78 0 : ResidualReductionFormatter() = default;
79 0 : ResidualReductionFormatter(std::string local_section_observation_key)
80 : : section_observation_key(std::move(local_section_observation_key)) {}
81 0 : std::string operator()(const size_t iteration_id,
82 : const double residual) const {
83 : if (iteration_id == 0) {
84 : return pretty_type::name<OptionsGroup>() + section_observation_key +
85 : " initialized with residual: " + get_output(residual) + "\n";
86 : } else {
87 : return pretty_type::name<OptionsGroup>() + section_observation_key + "(" +
88 : get_output(iteration_id) +
89 : ") iteration complete. Remaining residual: " +
90 : get_output(residual) + "\n";
91 : }
92 : }
93 : // NOLINTNEXTLINE(google-runtime-references)
94 0 : void pup(PUP::er& p) { p | section_observation_key; }
95 0 : std::string section_observation_key{};
96 : };
97 :
98 : template <typename OptionsGroup, typename ParallelComponent,
99 : typename Metavariables, typename ArrayIndex>
100 0 : void contribute_to_residual_observation(
101 : const size_t iteration_id, const double residual_magnitude_square,
102 : Parallel::GlobalCache<Metavariables>& cache, const ArrayIndex& array_index,
103 : const std::string& section_observation_key) {
104 : auto& local_observer = *Parallel::local_branch(
105 : Parallel::get_parallel_component<observers::Observer<Metavariables>>(
106 : cache));
107 : auto formatter =
108 : UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
109 : ::Verbosity::Quiet)
110 : ? std::make_optional(ResidualReductionFormatter<OptionsGroup>{
111 : section_observation_key})
112 : : std::nullopt;
113 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
114 : local_observer,
115 : observers::ObservationId(
116 : iteration_id,
117 : pretty_type::get_name<OptionsGroup>() + section_observation_key),
118 : Parallel::make_array_component_id<ParallelComponent>(array_index),
119 : std::string{"/" + pretty_type::name<OptionsGroup>() +
120 : section_observation_key + "Residuals"},
121 : std::vector<std::string>{"Iteration", "Residual"},
122 : reduction_data{iteration_id, residual_magnitude_square},
123 : std::move(formatter));
124 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
125 : ::Verbosity::Debug)) {
126 : if (iteration_id == 0) {
127 : Parallel::printf(
128 : "%s %s initialized with local residual: %e\n",
129 : get_output(array_index),
130 : pretty_type::name<OptionsGroup>() + section_observation_key,
131 : sqrt(residual_magnitude_square));
132 : } else {
133 : Parallel::printf(
134 : "%s %s(%zu) iteration complete. Remaining local residual: %e\n",
135 : get_output(array_index),
136 : pretty_type::name<OptionsGroup>() + section_observation_key,
137 : iteration_id, sqrt(residual_magnitude_square));
138 : }
139 : }
140 : }
141 :
142 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
143 0 : struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
144 : private:
145 0 : using fields_tag = FieldsTag;
146 0 : using operator_applied_to_fields_tag =
147 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, fields_tag>;
148 0 : using source_tag = SourceTag;
149 0 : using residual_tag =
150 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
151 0 : using residual_magnitude_square_tag =
152 : LinearSolver::Tags::MagnitudeSquare<residual_tag>;
153 :
154 : public: // Iterable action
155 0 : using const_global_cache_tags =
156 : tmpl::list<Convergence::Tags::Iterations<OptionsGroup>>;
157 :
158 0 : using simple_tags = tmpl::list<Convergence::Tags::IterationId<OptionsGroup>,
159 : Convergence::Tags::HasConverged<OptionsGroup>,
160 : operator_applied_to_fields_tag>;
161 0 : using compute_tags =
162 : tmpl::list<LinearSolver::Tags::ResidualCompute<fields_tag, source_tag>>;
163 :
164 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
165 : typename ArrayIndex, typename ActionList,
166 : typename ParallelComponent>
167 0 : static Parallel::iterable_action_return_t apply(
168 : db::DataBox<DbTagsList>& box,
169 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
170 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
171 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
172 : const ParallelComponent* const /*meta*/) {
173 : // The `PrepareSolve` action populates these tags with initial
174 : // values, except for `operator_applied_to_fields_tag` which is
175 : // expected to be updated in every iteration of the algorithm
176 : Initialization::mutate_assign<
177 : tmpl::list<Convergence::Tags::IterationId<OptionsGroup>>>(
178 : make_not_null(&box), std::numeric_limits<size_t>::max());
179 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
180 : }
181 :
182 : public: // amr::protocols::Projector
183 0 : using argument_tags = tmpl::list<>;
184 0 : using return_tags = simple_tags;
185 :
186 : // p-refinement
187 : template <size_t Dim>
188 0 : static void apply(
189 : const gsl::not_null<size_t*> /*iteration_id*/,
190 : const gsl::not_null<Convergence::HasConverged*> /*has_converged*/,
191 : const gsl::not_null<typename operator_applied_to_fields_tag::
192 : type*> /*operator_applied_to_fields*/,
193 : const std::pair<Mesh<Dim>, Element<Dim>>& /*old_mesh_and_element*/) {}
194 :
195 : // h-refinement
196 : template <typename... ParentTags>
197 0 : static void apply(
198 : const gsl::not_null<size_t*> iteration_id,
199 : const gsl::not_null<Convergence::HasConverged*> has_converged,
200 : const gsl::not_null<typename operator_applied_to_fields_tag::
201 : type*> /*operator_applied_to_fields*/,
202 : const tuples::TaggedTuple<ParentTags...>& parent_items) {
203 : *iteration_id =
204 : get<Convergence::Tags::IterationId<OptionsGroup>>(parent_items);
205 : *has_converged =
206 : get<Convergence::Tags::HasConverged<OptionsGroup>>(parent_items);
207 : }
208 :
209 : // h-coarsening
210 : template <size_t Dim, typename... ChildTags>
211 0 : static void apply(
212 : const gsl::not_null<size_t*> iteration_id,
213 : const gsl::not_null<Convergence::HasConverged*> has_converged,
214 : const gsl::not_null<typename operator_applied_to_fields_tag::
215 : type*> /*operator_applied_to_fields*/,
216 : const std::unordered_map<
217 : ElementId<Dim>, tuples::TaggedTuple<ChildTags...>>& children_items) {
218 : *iteration_id = get<Convergence::Tags::IterationId<OptionsGroup>>(
219 : children_items.begin()->second);
220 : *has_converged = get<Convergence::Tags::HasConverged<OptionsGroup>>(
221 : children_items.begin()->second);
222 : }
223 : };
224 :
225 : template <typename OptionsGroup, typename ArraySectionIdTag = void>
226 0 : struct RegisterObservers {
227 : template <typename ParallelComponent, typename DbTagsList,
228 : typename ArrayIndex>
229 : static std::pair<observers::TypeOfObservation, observers::ObservationKey>
230 0 : register_info(const db::DataBox<DbTagsList>& box,
231 : const ArrayIndex& /*array_index*/) {
232 : // Get the observation key, or "Unused" if the element does not belong
233 : // to a section with this tag. In the latter case, no observations will
234 : // ever be contributed.
235 : const std::optional<std::string> section_observation_key =
236 : observers::get_section_observation_key<ArraySectionIdTag>(box);
237 : ASSERT(section_observation_key != "Unused",
238 : "The identifier 'Unused' is reserved to indicate that no "
239 : "observations with this key will be contributed. Use a different "
240 : "key, or change the identifier 'Unused' to something else.");
241 : return {
242 : observers::TypeOfObservation::Reduction,
243 : observers::ObservationKey{pretty_type::get_name<OptionsGroup>() +
244 : section_observation_key.value_or("Unused")}};
245 : }
246 : };
247 :
248 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
249 : typename ArraySectionIdTag = void>
250 0 : using RegisterElement = observers::Actions::RegisterWithObservers<
251 : RegisterObservers<OptionsGroup, ArraySectionIdTag>>;
252 :
253 : /*!
254 : * \brief Prepare the asynchronous linear solver for a solve
255 : *
256 : * This action resets the asynchronous linear solver to its initial state, and
257 : * optionally observes the initial residual. If the initial residual should be
258 : * observed, both the `SourceTag` as well as the
259 : * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, FieldsTag>`
260 : * must be up-to-date at the time this action is invoked.
261 : *
262 : * This action also provides an anchor point in the action list for looping the
263 : * linear solver, in the sense that the algorithm jumps back to the action
264 : * immediately following this one when a step is complete and the solver hasn't
265 : * yet converged (see `LinearSolver::async_solvers::CompleteStep`).
266 : *
267 : * \tparam FieldsTag The data `x` in the linear equation `Ax = b` to be solved.
268 : * Should hold the initial guess `x_0` at this point in the algorithm.
269 : * \tparam OptionsGroup An options group identifying the linear solver
270 : * \tparam SourceTag The data `b` in `Ax = b`
271 : * \tparam Label An optional compile-time label for the solver to distinguish
272 : * different solves with the same solver in the action list
273 : * \tparam ArraySectionIdTag Observe the residual norm separately for each
274 : * array section identified by this tag (see `Parallel::Section`). Set to `void`
275 : * to observe the residual norm over all elements of the array (default). The
276 : * section only affects observations of residuals and has no effect on the
277 : * solver algorithm.
278 : * \tparam ObserveInitialResidual Whether or not to observe the initial residual
279 : * `b - A x_0` (default: `true`). Disable when `b` or `A x_0` are not yet
280 : * available at the preparation stage.
281 : */
282 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
283 : typename Label = OptionsGroup, typename ArraySectionIdTag = void,
284 : bool ObserveInitialResidual = true>
285 1 : struct PrepareSolve {
286 : private:
287 0 : using fields_tag = FieldsTag;
288 0 : using residual_tag =
289 : db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>;
290 :
291 : public:
292 0 : using const_global_cache_tags =
293 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
294 :
295 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
296 : typename ArrayIndex, typename ActionList,
297 : typename ParallelComponent>
298 0 : static Parallel::iterable_action_return_t apply(
299 : db::DataBox<DbTagsList>& box,
300 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
301 : Parallel::GlobalCache<Metavariables>& cache,
302 : const ArrayIndex& array_index, const ActionList /*meta*/,
303 : const ParallelComponent* const /*meta*/) {
304 : constexpr size_t iteration_id = 0;
305 :
306 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
307 : ::Verbosity::Debug)) {
308 : Parallel::printf("%s %s: Prepare solve\n", get_output(array_index),
309 : pretty_type::name<OptionsGroup>());
310 : }
311 :
312 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
313 : Convergence::Tags::HasConverged<OptionsGroup>>(
314 : [](const gsl::not_null<size_t*> local_iteration_id,
315 : const gsl::not_null<Convergence::HasConverged*> has_converged,
316 : const size_t num_iterations) {
317 : *local_iteration_id = iteration_id;
318 : *has_converged =
319 : Convergence::HasConverged{num_iterations, iteration_id};
320 : },
321 : make_not_null(&box),
322 : get<Convergence::Tags::Iterations<OptionsGroup>>(box));
323 :
324 : if constexpr (ObserveInitialResidual) {
325 : // Observe the initial residual even if no steps are going to be performed
326 : const std::optional<std::string> section_observation_key =
327 : observers::get_section_observation_key<ArraySectionIdTag>(box);
328 : if (section_observation_key.has_value()) {
329 : const auto& residual = get<residual_tag>(box);
330 : const double residual_magnitude_square = magnitude_square(residual);
331 : contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
332 : iteration_id, residual_magnitude_square, cache, array_index,
333 : *section_observation_key);
334 : }
335 : }
336 :
337 : // Skip steps entirely if the solve has already converged
338 : constexpr size_t step_end_index = tmpl::index_of<
339 : ActionList,
340 : CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
341 : ArraySectionIdTag, ObserveInitialResidual>>::value;
342 : constexpr size_t this_action_index =
343 : tmpl::index_of<ActionList, PrepareSolve>::value;
344 : return {Parallel::AlgorithmExecution::Continue,
345 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
346 : ? (step_end_index + 1)
347 : : (this_action_index + 1)};
348 : }
349 : };
350 :
351 : /*!
352 : * \brief Complete a step of the asynchronous linear solver
353 : *
354 : * This action prepares the next step of the asynchronous linear solver, and
355 : * observes the residual. To observe the correct residual, make sure the
356 : * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, FieldsTag>` is
357 : * up-to-date at the time this action is invoked.
358 : *
359 : * This action checks if the algorithm has converged, i.e. it has completed the
360 : * requested number of steps. If it hasn't, the algorithm jumps back to the
361 : * action immediately following the `LinearSolver::async_solvers::PrepareSolve`
362 : * to perform another iteration. Make sure both actions use the same template
363 : * parameters.
364 : *
365 : * \tparam FieldsTag The data `x` in the linear equation `Ax = b` to be solved.
366 : * \tparam OptionsGroup An options group identifying the linear solver
367 : * \tparam SourceTag The data `b` in `Ax = b`
368 : * \tparam Label An optional compile-time label for the solver to distinguish
369 : * different solves with the same solver in the action list
370 : * \tparam ArraySectionIdTag Observe the residual norm separately for each
371 : * array section identified by this tag (see `Parallel::Section`). Set to `void`
372 : * to observe the residual norm over all elements of the array (default). The
373 : * section only affects observations of residuals and has no effect on the
374 : * solver algorithm.
375 : * \tparam ObserveInitialResidual Whether or not to observe the _initial_
376 : * residual `b - A x_0`. This parameter should match the one passed to
377 : * `PrepareSolve`.
378 : */
379 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
380 : typename Label = OptionsGroup, typename ArraySectionIdTag = void,
381 : bool ObserveInitialResidual = true>
382 1 : struct CompleteStep {
383 : private:
384 0 : using fields_tag = FieldsTag;
385 0 : using residual_tag =
386 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
387 :
388 : public:
389 0 : using const_global_cache_tags =
390 : tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
391 :
392 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
393 : typename ArrayIndex, typename ActionList,
394 : typename ParallelComponent>
395 0 : static Parallel::iterable_action_return_t apply(
396 : db::DataBox<DbTagsList>& box,
397 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
398 : Parallel::GlobalCache<Metavariables>& cache,
399 : const ArrayIndex& array_index, const ActionList /*meta*/,
400 : const ParallelComponent* const /*meta*/) {
401 : // Prepare for next iteration
402 : db::mutate<Convergence::Tags::IterationId<OptionsGroup>,
403 : Convergence::Tags::HasConverged<OptionsGroup>>(
404 : [](const gsl::not_null<size_t*> iteration_id,
405 : const gsl::not_null<Convergence::HasConverged*> has_converged,
406 : const size_t num_iterations) {
407 : ++(*iteration_id);
408 : *has_converged =
409 : Convergence::HasConverged{num_iterations, *iteration_id};
410 : },
411 : make_not_null(&box),
412 : get<Convergence::Tags::Iterations<OptionsGroup>>(box));
413 :
414 : // Observe element-local residual magnitude
415 : const std::optional<std::string> section_observation_key =
416 : observers::get_section_observation_key<ArraySectionIdTag>(box);
417 : if (section_observation_key.has_value()) {
418 : const size_t completed_iterations =
419 : get<Convergence::Tags::IterationId<OptionsGroup>>(box);
420 : const auto& residual = get<residual_tag>(box);
421 : const double residual_magnitude_square = magnitude_square(residual);
422 : contribute_to_residual_observation<OptionsGroup, ParallelComponent>(
423 : completed_iterations, residual_magnitude_square, cache, array_index,
424 : *section_observation_key);
425 : }
426 :
427 : // Repeat steps until the solve has converged
428 : constexpr size_t step_begin_index =
429 : tmpl::index_of<
430 : ActionList,
431 : PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label,
432 : ArraySectionIdTag, ObserveInitialResidual>>::value +
433 : 1;
434 : constexpr size_t this_action_index =
435 : tmpl::index_of<ActionList, CompleteStep>::value;
436 : return {Parallel::AlgorithmExecution::Continue,
437 : get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
438 : ? (this_action_index + 1)
439 : : step_begin_index};
440 : }
441 : };
442 :
443 : } // namespace LinearSolver::async_solvers
|
