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 <map>
8 : #include <optional>
9 : #include <string>
10 : #include <tuple>
11 : #include <unordered_map>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 : #include "DataStructures/DataBox/Tag.hpp"
18 : #include "DataStructures/Index.hpp"
19 : #include "Domain/Structure/Direction.hpp"
20 : #include "Domain/Structure/DirectionMap.hpp"
21 : #include "Domain/Structure/Element.hpp"
22 : #include "Domain/Structure/ElementId.hpp"
23 : #include "Domain/Structure/OrientationMapHelpers.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "Domain/Tags/Faces.hpp"
26 : #include "IO/Logging/Tags.hpp"
27 : #include "IO/Logging/Verbosity.hpp"
28 : #include "IO/Observer/Actions/RegisterWithObservers.hpp"
29 : #include "IO/Observer/GetSectionObservationKey.hpp"
30 : #include "IO/Observer/ObservationId.hpp"
31 : #include "IO/Observer/Protocols/ReductionDataFormatter.hpp"
32 : #include "IO/Observer/ReductionActions.hpp"
33 : #include "IO/Observer/Tags.hpp"
34 : #include "IO/Observer/TypeOfObservation.hpp"
35 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
36 : #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
37 : #include "NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp"
38 : #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
39 : #include "NumericalAlgorithms/Spectral/CollocationPoints.hpp"
40 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
41 : #include "Parallel/AlgorithmExecution.hpp"
42 : #include "Parallel/ArrayComponentId.hpp"
43 : #include "Parallel/GlobalCache.hpp"
44 : #include "Parallel/InboxInserters.hpp"
45 : #include "Parallel/Invoke.hpp"
46 : #include "Parallel/Local.hpp"
47 : #include "Parallel/ParallelComponentHelpers.hpp"
48 : #include "Parallel/Printf/Printf.hpp"
49 : #include "Parallel/Reduction.hpp"
50 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
51 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
52 : #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
53 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp"
54 : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
55 : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
56 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
57 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Weighting.hpp"
58 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
59 : #include "Utilities/Gsl.hpp"
60 : #include "Utilities/PrettyType.hpp"
61 : #include "Utilities/ProtocolHelpers.hpp"
62 : #include "Utilities/SetNumberOfGridPoints.hpp"
63 : #include "Utilities/TMPL.hpp"
64 : #include "Utilities/TaggedTuple.hpp"
65 : #include "Utilities/TypeTraits/IsA.hpp"
66 :
67 : namespace LinearSolver::Schwarz::detail {
68 :
69 : using reduction_data = Parallel::ReductionData<
70 : // Iteration
71 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
72 : // Number of subdomains (= number of elements)
73 : Parallel::ReductionDatum<size_t, funcl::Plus<>>,
74 : // Average number of subdomain solver iterations
75 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
76 : std::index_sequence<1>>,
77 : // Minimum number of subdomain solver iterations
78 : Parallel::ReductionDatum<size_t, funcl::Min<>>,
79 : // Maximum number of subdomain solver iterations
80 : Parallel::ReductionDatum<size_t, funcl::Max<>>,
81 : // Total number of subdomain solver iterations
82 : Parallel::ReductionDatum<size_t, funcl::Plus<>>>;
83 :
84 : template <typename OptionsGroup>
85 : struct SubdomainStatsFormatter
86 : : tt::ConformsTo<observers::protocols::ReductionDataFormatter> {
87 : using reduction_data = Schwarz::detail::reduction_data;
88 : SubdomainStatsFormatter() = default;
89 : SubdomainStatsFormatter(std::string local_section_observation_key)
90 : : section_observation_key(std::move(local_section_observation_key)) {}
91 : std::string operator()(const size_t iteration_id, const size_t num_subdomains,
92 : const double avg_subdomain_its,
93 : const size_t min_subdomain_its,
94 : const size_t max_subdomain_its,
95 : const size_t total_subdomain_its) const {
96 : return pretty_type::name<OptionsGroup>() + section_observation_key + "(" +
97 : get_output(iteration_id) + ") completed all " +
98 : get_output(num_subdomains) +
99 : " subdomain solves. Average of number of iterations: " +
100 : get_output(avg_subdomain_its) + " (min " +
101 : get_output(min_subdomain_its) + ", max " +
102 : get_output(max_subdomain_its) + ", total " +
103 : get_output(total_subdomain_its) + ").\n";
104 : }
105 : // NOLINTNEXTLINE(google-runtime-references)
106 : void pup(PUP::er& p) { p | section_observation_key; }
107 : std::string section_observation_key{};
108 : };
109 :
110 : template <typename OptionsGroup, typename ArraySectionIdTag>
111 : struct RegisterObservers {
112 : template <typename ParallelComponent, typename DbTagsList,
113 : typename ArrayIndex>
114 : static std::pair<observers::TypeOfObservation, observers::ObservationKey>
115 : register_info(const db::DataBox<DbTagsList>& box,
116 : const ArrayIndex& /*array_index*/) {
117 : // Get the observation key, or "Unused" if the element does not belong
118 : // to a section with this tag. In the latter case, no observations will
119 : // ever be contributed.
120 : const std::optional<std::string> section_observation_key =
121 : observers::get_section_observation_key<ArraySectionIdTag>(box);
122 : ASSERT(section_observation_key != "Unused",
123 : "The identifier 'Unused' is reserved to indicate that no "
124 : "observations with this key will be contributed. Use a different "
125 : "key, or change the identifier 'Unused' to something else.");
126 : return {
127 : observers::TypeOfObservation::Reduction,
128 : observers::ObservationKey{pretty_type::get_name<OptionsGroup>() +
129 : section_observation_key.value_or("Unused") +
130 : "SubdomainSolves"}};
131 : }
132 : };
133 :
134 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
135 : typename ArraySectionIdTag>
136 : using RegisterElement = observers::Actions::RegisterWithObservers<
137 : RegisterObservers<OptionsGroup, ArraySectionIdTag>>;
138 :
139 : template <typename OptionsGroup, typename ParallelComponent,
140 : typename Metavariables, typename ArrayIndex>
141 : void contribute_to_subdomain_stats_observation(
142 : const size_t iteration_id, const size_t subdomain_solve_num_iterations,
143 : Parallel::GlobalCache<Metavariables>& cache, const ArrayIndex& array_index,
144 : const std::string& section_observation_key, const bool observe_per_core) {
145 : auto& local_observer = *Parallel::local_branch(
146 : Parallel::get_parallel_component<observers::Observer<Metavariables>>(
147 : cache));
148 :
149 : auto formatter =
150 : UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
151 : ::Verbosity::Verbose)
152 : ? std::make_optional(
153 : SubdomainStatsFormatter<OptionsGroup>{section_observation_key})
154 : : std::nullopt;
155 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
156 : local_observer,
157 : observers::ObservationId(iteration_id,
158 : pretty_type::get_name<OptionsGroup>() +
159 : section_observation_key + "SubdomainSolves"),
160 : Parallel::make_array_component_id<ParallelComponent>(array_index),
161 : std::string{"/" + pretty_type::name<OptionsGroup>() +
162 : section_observation_key + "SubdomainSolves"},
163 : std::vector<std::string>{"Iteration", "NumSubdomains", "AvgNumIterations",
164 : "MinNumIterations", "MaxNumIterations",
165 : "TotalNumIterations"},
166 : reduction_data{
167 : iteration_id, 1, static_cast<double>(subdomain_solve_num_iterations),
168 : subdomain_solve_num_iterations, subdomain_solve_num_iterations,
169 : subdomain_solve_num_iterations},
170 : std::move(formatter), observe_per_core);
171 : }
172 :
173 : template <typename SubdomainDataType, typename OptionsGroup>
174 : struct SubdomainDataBufferTag : db::SimpleTag {
175 : static std::string name() {
176 : return "SubdomainData(" + pretty_type::name<OptionsGroup>() + ")";
177 : }
178 : using type = SubdomainDataType;
179 : };
180 :
181 : // Allow factory-creating any of these serial linear solvers for use as
182 : // subdomain solver
183 : template <typename FieldsTag, typename SubdomainOperator,
184 : typename SubdomainPreconditioners,
185 : typename SubdomainData = ElementCenteredSubdomainData<
186 : SubdomainOperator::volume_dim,
187 : typename db::add_tag_prefix<LinearSolver::Tags::Residual,
188 : FieldsTag>::tags_list>>
189 : using subdomain_solver = LinearSolver::Serial::LinearSolver<tmpl::append<
190 : tmpl::list<::LinearSolver::Serial::Registrars::Gmres<SubdomainData>,
191 : ::LinearSolver::Serial::Registrars::ExplicitInverse<
192 : typename SubdomainData::value_type>>,
193 : SubdomainPreconditioners>>;
194 :
195 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
196 : typename SubdomainPreconditioners>
197 : struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
198 : private:
199 : using fields_tag = FieldsTag;
200 : using residual_tag =
201 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
202 : static constexpr size_t Dim = SubdomainOperator::volume_dim;
203 : using SubdomainData =
204 : ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
205 : using SubdomainSolver =
206 : subdomain_solver<FieldsTag, SubdomainOperator, SubdomainPreconditioners>;
207 : using subdomain_solver_tag =
208 : Tags::SubdomainSolver<SubdomainSolver, OptionsGroup>;
209 : using volume_data_tag =
210 : Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>;
211 :
212 : public: // Iterable action
213 : using simple_tags_from_options = tmpl::list<subdomain_solver_tag>;
214 :
215 : using const_global_cache_tags =
216 : tmpl::list<Tags::MaxOverlap<OptionsGroup>,
217 : LinearSolver::Tags::OutputVolumeData<OptionsGroup>>;
218 :
219 : using simple_tags = tmpl::list<
220 : Tags::IntrudingExtents<Dim, OptionsGroup>, Tags::Weight<OptionsGroup>,
221 : domain::Tags::Faces<Dim, Tags::Weight<OptionsGroup>>,
222 : SubdomainDataBufferTag<SubdomainData, OptionsGroup>,
223 : LinearSolver::Tags::ObservationId<OptionsGroup>, volume_data_tag>;
224 : using compute_tags = tmpl::list<>;
225 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
226 : typename ActionList, typename ParallelComponent>
227 : static Parallel::iterable_action_return_t apply(
228 : db::DataBox<DbTagsList>& box,
229 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
230 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
231 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
232 : const ParallelComponent* const /*meta*/) {
233 : db::mutate_apply<InitializeElement>(make_not_null(&box));
234 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
235 : }
236 :
237 : public: // amr::protocols::Projector
238 : using return_tags = tmpl::append<simple_tags, simple_tags_from_options>;
239 : using argument_tags =
240 : tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>,
241 : domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
242 : Tags::MaxOverlap<OptionsGroup>>;
243 :
244 : template <typename... AmrData>
245 : static void apply(
246 : const gsl::not_null<std::array<size_t, Dim>*> intruding_extents,
247 : const gsl::not_null<Scalar<DataVector>*> element_weight,
248 : const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
249 : intruding_overlap_weights,
250 : const gsl::not_null<SubdomainData*> subdomain_data,
251 : const gsl::not_null<size_t*> observation_id,
252 : const gsl::not_null<SubdomainData*> /*volume_data_for_output*/,
253 : const gsl::not_null<std::unique_ptr<SubdomainSolver>*> subdomain_solver,
254 : const Element<Dim>& element, const Mesh<Dim>& mesh,
255 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>& logical_coords,
256 : const std::optional<size_t> max_overlap, const AmrData&... amr_data) {
257 : const size_t num_points = mesh.number_of_grid_points();
258 :
259 : // Intruding overlaps
260 : std::array<double, Dim> intruding_overlap_widths{};
261 : for (size_t d = 0; d < Dim; ++d) {
262 : gsl::at(*intruding_extents, d) =
263 : LinearSolver::Schwarz::overlap_extent(mesh.extents(d), max_overlap);
264 : const auto& collocation_points =
265 : Spectral::collocation_points(mesh.slice_through(d));
266 : gsl::at(intruding_overlap_widths, d) =
267 : LinearSolver::Schwarz::overlap_width(gsl::at(*intruding_extents, d),
268 : collocation_points);
269 : }
270 :
271 : // Element weight
272 : // For max_overlap > 0 all overlaps will have non-zero extents on an LGL
273 : // mesh (because it has at least 2 points per dimension), so we don't need
274 : // to check their extents are non-zero individually
275 : if (LIKELY(max_overlap != 0)) {
276 : LinearSolver::Schwarz::element_weight(element_weight, logical_coords,
277 : intruding_overlap_widths,
278 : element.external_boundaries());
279 : } else {
280 : set_number_of_grid_points(element_weight, num_points);
281 : get(*element_weight) = 1.;
282 : }
283 :
284 : // Intruding overlap weights
285 : intruding_overlap_weights->clear();
286 : for (const auto& direction : element.internal_boundaries()) {
287 : const size_t dim = direction.dimension();
288 : if (gsl::at(*intruding_extents, dim) > 0) {
289 : const auto intruding_logical_coords =
290 : LinearSolver::Schwarz::data_on_overlap(
291 : logical_coords, mesh.extents(),
292 : gsl::at(*intruding_extents, dim), direction);
293 : (*intruding_overlap_weights)[direction] =
294 : LinearSolver::Schwarz::intruding_weight(
295 : intruding_logical_coords, direction, intruding_overlap_widths,
296 : element.neighbors().at(direction).size(),
297 : element.external_boundaries());
298 : }
299 : }
300 :
301 : // Subdomain data buffer
302 : *subdomain_data = SubdomainData{num_points};
303 :
304 : // Subdomain solver
305 : // The subdomain solver initially gets created from options on each element.
306 : // Then we have to copy it around during AMR.
307 : if constexpr (sizeof...(AmrData) == 0) {
308 : *observation_id = 0;
309 : (void)subdomain_solver;
310 : } else {
311 : if constexpr (tt::is_a_v<tuples::TaggedTuple, AmrData...>) {
312 : // h-refinement: copy from the parent
313 : *observation_id =
314 : get<LinearSolver::Tags::ObservationId<OptionsGroup>>(amr_data...);
315 : *subdomain_solver = get<subdomain_solver_tag>(amr_data...)->get_clone();
316 : } else if constexpr (tt::is_a_v<std::unordered_map, AmrData...>) {
317 : // h-coarsening, copy from one of the children (doesn't matter which)
318 : *observation_id = get<LinearSolver::Tags::ObservationId<OptionsGroup>>(
319 : amr_data.begin()->second...);
320 : *subdomain_solver =
321 : get<subdomain_solver_tag>(amr_data.begin()->second...)->get_clone();
322 : } else {
323 : (void)observation_id;
324 : }
325 : (*subdomain_solver)->reset();
326 : }
327 : }
328 : };
329 :
330 : // Restrict the residual to neighboring subdomains that overlap with this
331 : // element and send the data to those elements
332 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator>
333 : using SendOverlapData = LinearSolver::Schwarz::Actions::SendOverlapFields<
334 : tmpl::list<db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>>,
335 : OptionsGroup, true>;
336 :
337 : template <size_t Dim, typename OptionsGroup, typename OverlapSolution>
338 : struct OverlapSolutionInboxTag
339 : : public Parallel::InboxInserters::Map<
340 : OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapSolution>> {
341 : using temporal_id = size_t;
342 : using type = std::map<temporal_id, OverlapMap<Dim, OverlapSolution>>;
343 : };
344 :
345 : // Wait for the residual data on regions of this element's subdomain that
346 : // overlap with other elements. Once the residual data is available on all
347 : // overlaps, solve the restricted problem for this element-centered subdomain.
348 : // Apply the weighted solution on this element directly and send the solution on
349 : // overlap regions to the neighbors that they overlap with.
350 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
351 : typename SubdomainPreconditioners, typename ArraySectionIdTag>
352 : struct SolveSubdomain {
353 : private:
354 : using fields_tag = FieldsTag;
355 : using residual_tag =
356 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
357 : static constexpr size_t Dim = SubdomainOperator::volume_dim;
358 : using overlap_residuals_inbox_tag =
359 : Actions::detail::OverlapFieldsTag<Dim, tmpl::list<residual_tag>,
360 : OptionsGroup>;
361 : using SubdomainData =
362 : ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
363 : using SubdomainSolver =
364 : subdomain_solver<FieldsTag, SubdomainOperator, SubdomainPreconditioners>;
365 : using subdomain_solver_tag =
366 : Tags::SubdomainSolver<SubdomainSolver, OptionsGroup>;
367 : using OverlapData = typename SubdomainData::OverlapData;
368 : using overlap_solution_inbox_tag =
369 : OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapData>;
370 :
371 : public:
372 : using const_global_cache_tags =
373 : tmpl::list<Tags::MaxOverlap<OptionsGroup>,
374 : logging::Tags::Verbosity<OptionsGroup>,
375 : Tags::ObservePerCoreReductions<OptionsGroup>>;
376 : using inbox_tags = tmpl::list<overlap_residuals_inbox_tag>;
377 :
378 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
379 : typename ActionList, typename ParallelComponent>
380 : static Parallel::iterable_action_return_t apply(
381 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
382 : Parallel::GlobalCache<Metavariables>& cache,
383 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
384 : const ParallelComponent* const /*meta*/) {
385 : const size_t iteration_id =
386 : get<Convergence::Tags::IterationId<OptionsGroup>>(box);
387 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
388 : const auto& max_overlap = db::get<Tags::MaxOverlap<OptionsGroup>>(box);
389 :
390 : // Wait for communicated overlap data
391 : const bool has_overlap_data =
392 : max_overlap != 0 and element.number_of_neighbors() > 0;
393 : if (LIKELY(has_overlap_data) and
394 : not dg::has_received_from_all_mortars<overlap_residuals_inbox_tag>(
395 : iteration_id, element, inboxes)) {
396 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
397 : }
398 :
399 : // Do some logging
400 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
401 : ::Verbosity::Debug)) {
402 : Parallel::printf("%s %s(%zu): Solve subdomain\n", element_id,
403 : pretty_type::name<OptionsGroup>(), iteration_id);
404 : }
405 :
406 : // Assemble the subdomain data from the data on the element and the
407 : // communicated overlap data
408 : db::mutate<SubdomainDataBufferTag<SubdomainData, OptionsGroup>>(
409 : [&inboxes, &iteration_id, &has_overlap_data](
410 : const gsl::not_null<SubdomainData*> subdomain_data,
411 : const auto& residual) {
412 : subdomain_data->element_data = residual;
413 : // Nothing was communicated if the overlaps are empty
414 : if (LIKELY(has_overlap_data)) {
415 : subdomain_data->overlap_data =
416 : std::move(tuples::get<overlap_residuals_inbox_tag>(inboxes)
417 : .extract(iteration_id)
418 : .mapped());
419 : }
420 : },
421 : make_not_null(&box), db::get<residual_tag>(box));
422 : const auto& subdomain_residual =
423 : db::get<SubdomainDataBufferTag<SubdomainData, OptionsGroup>>(box);
424 :
425 : // Allocate workspace memory for repeatedly applying the subdomain operator
426 : const SubdomainOperator subdomain_operator{};
427 :
428 : // Solve the subdomain problem
429 : const auto& subdomain_solver = get<subdomain_solver_tag>(box);
430 : auto subdomain_solve_initial_guess_in_solution_out =
431 : make_with_value<SubdomainData>(subdomain_residual, 0.);
432 : const auto subdomain_solve_has_converged = subdomain_solver.solve(
433 : make_not_null(&subdomain_solve_initial_guess_in_solution_out),
434 : subdomain_operator, subdomain_residual, std::forward_as_tuple(box));
435 : // Re-naming the solution buffer for the code below
436 : auto& subdomain_solution = subdomain_solve_initial_guess_in_solution_out;
437 :
438 : // Do some logging and observing
439 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
440 : ::Verbosity::Quiet)) {
441 : if (not subdomain_solve_has_converged or
442 : subdomain_solve_has_converged.reason() ==
443 : Convergence::Reason::MaxIterations) {
444 : Parallel::printf(
445 : "%s %s(%zu): WARNING: Subdomain solver did not converge in %zu "
446 : "iterations: %e -> %e\n",
447 : element_id, pretty_type::name<OptionsGroup>(), iteration_id,
448 : subdomain_solve_has_converged.num_iterations(),
449 : subdomain_solve_has_converged.initial_residual_magnitude(),
450 : subdomain_solve_has_converged.residual_magnitude());
451 : } else if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
452 : ::Verbosity::Debug)) {
453 : Parallel::printf(
454 : "%s %s(%zu): Subdomain solver converged in %zu iterations (%s): %e "
455 : "-> %e\n",
456 : element_id, pretty_type::name<OptionsGroup>(), iteration_id,
457 : subdomain_solve_has_converged.num_iterations(),
458 : subdomain_solve_has_converged.reason(),
459 : subdomain_solve_has_converged.initial_residual_magnitude(),
460 : subdomain_solve_has_converged.residual_magnitude());
461 : }
462 : }
463 : const std::optional<std::string> section_observation_key =
464 : observers::get_section_observation_key<ArraySectionIdTag>(box);
465 : if (section_observation_key.has_value()) {
466 : contribute_to_subdomain_stats_observation<OptionsGroup,
467 : ParallelComponent>(
468 : iteration_id + 1, subdomain_solve_has_converged.num_iterations(),
469 : cache, element_id, *section_observation_key,
470 : db::get<Tags::ObservePerCoreReductions<OptionsGroup>>(box));
471 : }
472 :
473 : // Record subdomain solution for volume data output
474 : if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
475 : db::mutate<Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>>(
476 : [&subdomain_solution](const auto volume_data) {
477 : volume_data->element_data = subdomain_solution.element_data;
478 : },
479 : make_not_null(&box));
480 : }
481 :
482 : // Apply weighting
483 : if (LIKELY(max_overlap != 0)) {
484 : subdomain_solution.element_data *=
485 : get(db::get<Tags::Weight<OptionsGroup>>(box));
486 : }
487 :
488 : // Apply solution to central element
489 : db::mutate<fields_tag>(
490 : [&subdomain_solution](const auto fields) {
491 : *fields += subdomain_solution.element_data;
492 : },
493 : make_not_null(&box));
494 :
495 : // Send overlap solutions back to the neighbors that they are on
496 : if (LIKELY(max_overlap != 0)) {
497 : auto& receiver_proxy =
498 : Parallel::get_parallel_component<ParallelComponent>(cache);
499 : for (auto& [overlap_id, overlap_solution] :
500 : subdomain_solution.overlap_data) {
501 : const auto& direction = overlap_id.direction();
502 : const auto& neighbor_id = overlap_id.id();
503 : const auto& orientation =
504 : element.neighbors().at(direction).orientation(neighbor_id);
505 : const auto direction_from_neighbor = orientation(direction.opposite());
506 : Parallel::receive_data<overlap_solution_inbox_tag>(
507 : receiver_proxy[neighbor_id], iteration_id,
508 : std::make_pair(
509 : OverlapId<Dim>{direction_from_neighbor, element.id()},
510 : std::move(overlap_solution)));
511 : }
512 : }
513 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
514 : }
515 : };
516 :
517 : // Wait for the subdomain solutions on regions within this element that overlap
518 : // with neighboring element-centered subdomains. Combine the solutions as a
519 : // weighted sum.
520 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator>
521 : struct ReceiveOverlapSolution {
522 : private:
523 : using fields_tag = FieldsTag;
524 : using residual_tag =
525 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
526 : static constexpr size_t Dim = SubdomainOperator::volume_dim;
527 : using SubdomainData =
528 : ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
529 : using OverlapSolution = typename SubdomainData::OverlapData;
530 : using overlap_solution_inbox_tag =
531 : OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapSolution>;
532 :
533 : public:
534 : using const_global_cache_tags = tmpl::list<Tags::MaxOverlap<OptionsGroup>>;
535 : using inbox_tags = tmpl::list<overlap_solution_inbox_tag>;
536 :
537 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
538 : typename ActionList, typename ParallelComponent>
539 : static Parallel::iterable_action_return_t apply(
540 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
541 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
542 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
543 : const ParallelComponent* const /*meta*/) {
544 : const size_t iteration_id =
545 : get<Convergence::Tags::IterationId<OptionsGroup>>(box);
546 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
547 :
548 : // Nothing to do if overlap is empty
549 : if (UNLIKELY(db::get<Tags::MaxOverlap<OptionsGroup>>(box) == 0 or
550 : element.number_of_neighbors() == 0)) {
551 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
552 : }
553 :
554 : if (not dg::has_received_from_all_mortars<overlap_solution_inbox_tag>(
555 : iteration_id, element, inboxes)) {
556 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
557 : }
558 :
559 : // Do some logging
560 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
561 : ::Verbosity::Debug)) {
562 : Parallel::printf("%s %s(%zu): Receive overlap solution\n", element_id,
563 : pretty_type::name<OptionsGroup>(), iteration_id);
564 : }
565 :
566 : // Extract overlap solutions from inbox
567 : auto received_overlap_solutions =
568 : std::move(tuples::get<overlap_solution_inbox_tag>(inboxes)
569 : .extract(iteration_id)
570 : .mapped());
571 :
572 : // Record overlap solutions for volume data output
573 : if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
574 : db::mutate<Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>>(
575 : [&received_overlap_solutions](const auto volume_data) {
576 : volume_data->overlap_data = received_overlap_solutions;
577 : },
578 : make_not_null(&box));
579 : }
580 :
581 : // Add overlap solutions to this element's solution in a weighted sum
582 : db::mutate<fields_tag>(
583 : [&received_overlap_solutions](
584 : const auto fields, const Index<Dim>& full_extents,
585 : const std::array<size_t, Dim>& all_intruding_extents,
586 : const DirectionMap<Dim, Scalar<DataVector>>&
587 : all_intruding_overlap_weights) {
588 : for (auto& [overlap_id, overlap_solution] :
589 : received_overlap_solutions) {
590 : const auto& direction = overlap_id.direction();
591 : const auto& intruding_extents =
592 : gsl::at(all_intruding_extents, direction.dimension());
593 : const auto& overlap_weight =
594 : all_intruding_overlap_weights.at(direction);
595 : overlap_solution *= get(overlap_weight);
596 : LinearSolver::Schwarz::add_overlap_data(
597 : fields, overlap_solution, full_extents, intruding_extents,
598 : direction);
599 : }
600 : },
601 : make_not_null(&box), db::get<domain::Tags::Mesh<Dim>>(box).extents(),
602 : db::get<Tags::IntrudingExtents<Dim, OptionsGroup>>(box),
603 : db::get<domain::Tags::Faces<Dim, Tags::Weight<OptionsGroup>>>(box));
604 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
605 : }
606 : };
607 :
608 : } // namespace LinearSolver::Schwarz::detail
|
