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/Mesh.hpp"
40 : #include "NumericalAlgorithms/Spectral/Spectral.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) + ").";
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<std::unique_ptr<SubdomainSolver>, OptionsGroup>;
209 :
210 : public: // Iterable action
211 : using simple_tags_from_options = tmpl::list<subdomain_solver_tag>;
212 :
213 : using const_global_cache_tags = tmpl::list<Tags::MaxOverlap<OptionsGroup>>;
214 :
215 : using simple_tags =
216 : tmpl::list<Tags::IntrudingExtents<Dim, OptionsGroup>,
217 : Tags::Weight<OptionsGroup>,
218 : domain::Tags::Faces<Dim, Tags::Weight<OptionsGroup>>,
219 : SubdomainDataBufferTag<SubdomainData, OptionsGroup>>;
220 : using compute_tags = tmpl::list<>;
221 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
222 : typename ActionList, typename ParallelComponent>
223 : static Parallel::iterable_action_return_t apply(
224 : db::DataBox<DbTagsList>& box,
225 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
226 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
227 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
228 : const ParallelComponent* const /*meta*/) {
229 : db::mutate_apply<InitializeElement>(make_not_null(&box));
230 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
231 : }
232 :
233 : public: // amr::protocols::Projector
234 : using return_tags = tmpl::append<simple_tags, simple_tags_from_options>;
235 : using argument_tags =
236 : tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>,
237 : domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
238 : Tags::MaxOverlap<OptionsGroup>>;
239 :
240 : template <typename... AmrData>
241 : static void apply(
242 : const gsl::not_null<std::array<size_t, Dim>*> intruding_extents,
243 : const gsl::not_null<Scalar<DataVector>*> element_weight,
244 : const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
245 : intruding_overlap_weights,
246 : const gsl::not_null<SubdomainData*> subdomain_data,
247 : [[maybe_unused]] const gsl::not_null<std::unique_ptr<SubdomainSolver>*>
248 : subdomain_solver,
249 : const Element<Dim>& element, const Mesh<Dim>& mesh,
250 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>& logical_coords,
251 : const size_t max_overlap, const AmrData&... amr_data) {
252 : const size_t num_points = mesh.number_of_grid_points();
253 :
254 : // Intruding overlaps
255 : std::array<double, Dim> intruding_overlap_widths{};
256 : for (size_t d = 0; d < Dim; ++d) {
257 : gsl::at(*intruding_extents, d) =
258 : LinearSolver::Schwarz::overlap_extent(mesh.extents(d), max_overlap);
259 : const auto& collocation_points =
260 : Spectral::collocation_points(mesh.slice_through(d));
261 : gsl::at(intruding_overlap_widths, d) =
262 : LinearSolver::Schwarz::overlap_width(gsl::at(*intruding_extents, d),
263 : collocation_points);
264 : }
265 :
266 : // Element weight
267 : // For max_overlap > 0 all overlaps will have non-zero extents on an LGL
268 : // mesh (because it has at least 2 points per dimension), so we don't need
269 : // to check their extents are non-zero individually
270 : if (LIKELY(max_overlap > 0)) {
271 : LinearSolver::Schwarz::element_weight(element_weight, logical_coords,
272 : intruding_overlap_widths,
273 : element.external_boundaries());
274 : } else {
275 : set_number_of_grid_points(element_weight, num_points);
276 : get(*element_weight) = 1.;
277 : }
278 :
279 : // Intruding overlap weights
280 : intruding_overlap_weights->clear();
281 : for (const auto& direction : element.internal_boundaries()) {
282 : const size_t dim = direction.dimension();
283 : if (gsl::at(*intruding_extents, dim) > 0) {
284 : const auto intruding_logical_coords =
285 : LinearSolver::Schwarz::data_on_overlap(
286 : logical_coords, mesh.extents(),
287 : gsl::at(*intruding_extents, dim), direction);
288 : (*intruding_overlap_weights)[direction] =
289 : LinearSolver::Schwarz::intruding_weight(
290 : intruding_logical_coords, direction, intruding_overlap_widths,
291 : element.neighbors().at(direction).size(),
292 : element.external_boundaries());
293 : }
294 : }
295 :
296 : // Subdomain data buffer
297 : *subdomain_data = SubdomainData{num_points};
298 :
299 : // Subdomain solver
300 : // The subdomain solver initially gets created from options on each element.
301 : // Then we have to copy it around during AMR.
302 : if constexpr (sizeof...(AmrData) == 1) {
303 : if constexpr (tt::is_a_v<tuples::TaggedTuple, AmrData...>) {
304 : // h-refinement: copy from the parent
305 : *subdomain_solver = get<subdomain_solver_tag>(amr_data...)->get_clone();
306 : } else if constexpr (tt::is_a_v<std::unordered_map, AmrData...>) {
307 : // h-coarsening, copy from one of the children (doesn't matter which)
308 : *subdomain_solver =
309 : get<subdomain_solver_tag>(amr_data.begin()->second...)->get_clone();
310 : }
311 : (*subdomain_solver)->reset();
312 : }
313 : }
314 : };
315 :
316 : // Restrict the residual to neighboring subdomains that overlap with this
317 : // element and send the data to those elements
318 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator>
319 : using SendOverlapData = LinearSolver::Schwarz::Actions::SendOverlapFields<
320 : tmpl::list<db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>>,
321 : OptionsGroup, true>;
322 :
323 : template <size_t Dim, typename OptionsGroup, typename OverlapSolution>
324 : struct OverlapSolutionInboxTag
325 : : public Parallel::InboxInserters::Map<
326 : OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapSolution>> {
327 : using temporal_id = size_t;
328 : using type = std::map<temporal_id, OverlapMap<Dim, OverlapSolution>>;
329 : };
330 :
331 : // Wait for the residual data on regions of this element's subdomain that
332 : // overlap with other elements. Once the residual data is available on all
333 : // overlaps, solve the restricted problem for this element-centered subdomain.
334 : // Apply the weighted solution on this element directly and send the solution on
335 : // overlap regions to the neighbors that they overlap with.
336 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
337 : typename ArraySectionIdTag>
338 : struct SolveSubdomain {
339 : private:
340 : using fields_tag = FieldsTag;
341 : using residual_tag =
342 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
343 : static constexpr size_t Dim = SubdomainOperator::volume_dim;
344 : using overlap_residuals_inbox_tag =
345 : Actions::detail::OverlapFieldsTag<Dim, tmpl::list<residual_tag>,
346 : OptionsGroup>;
347 : using SubdomainData =
348 : ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
349 : using OverlapData = typename SubdomainData::OverlapData;
350 : using overlap_solution_inbox_tag =
351 : OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapData>;
352 :
353 : public:
354 : using const_global_cache_tags =
355 : tmpl::list<Tags::MaxOverlap<OptionsGroup>,
356 : logging::Tags::Verbosity<OptionsGroup>,
357 : Tags::ObservePerCoreReductions<OptionsGroup>>;
358 : using inbox_tags = tmpl::list<overlap_residuals_inbox_tag>;
359 :
360 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
361 : typename ActionList, typename ParallelComponent>
362 : static Parallel::iterable_action_return_t apply(
363 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
364 : Parallel::GlobalCache<Metavariables>& cache,
365 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
366 : const ParallelComponent* const /*meta*/) {
367 : const size_t iteration_id =
368 : get<Convergence::Tags::IterationId<OptionsGroup>>(box);
369 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
370 : const size_t max_overlap = db::get<Tags::MaxOverlap<OptionsGroup>>(box);
371 :
372 : // Wait for communicated overlap data
373 : const bool has_overlap_data =
374 : max_overlap > 0 and element.number_of_neighbors() > 0;
375 : if (LIKELY(has_overlap_data) and
376 : not dg::has_received_from_all_mortars<overlap_residuals_inbox_tag>(
377 : iteration_id, element, inboxes)) {
378 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
379 : }
380 :
381 : // Do some logging
382 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
383 : ::Verbosity::Debug)) {
384 : Parallel::printf("%s %s(%zu): Solve subdomain\n", element_id,
385 : pretty_type::name<OptionsGroup>(), iteration_id);
386 : }
387 :
388 : // Assemble the subdomain data from the data on the element and the
389 : // communicated overlap data
390 : db::mutate<SubdomainDataBufferTag<SubdomainData, OptionsGroup>>(
391 : [&inboxes, &iteration_id, &has_overlap_data](
392 : const gsl::not_null<SubdomainData*> subdomain_data,
393 : const auto& residual) {
394 : subdomain_data->element_data = residual;
395 : // Nothing was communicated if the overlaps are empty
396 : if (LIKELY(has_overlap_data)) {
397 : subdomain_data->overlap_data =
398 : std::move(tuples::get<overlap_residuals_inbox_tag>(inboxes)
399 : .extract(iteration_id)
400 : .mapped());
401 : }
402 : },
403 : make_not_null(&box), db::get<residual_tag>(box));
404 : const auto& subdomain_residual =
405 : db::get<SubdomainDataBufferTag<SubdomainData, OptionsGroup>>(box);
406 :
407 : // Allocate workspace memory for repeatedly applying the subdomain operator
408 : const SubdomainOperator subdomain_operator{};
409 :
410 : // Solve the subdomain problem
411 : const auto& subdomain_solver =
412 : get<Tags::SubdomainSolverBase<OptionsGroup>>(box);
413 : auto subdomain_solve_initial_guess_in_solution_out =
414 : make_with_value<SubdomainData>(subdomain_residual, 0.);
415 : const auto subdomain_solve_has_converged = subdomain_solver.solve(
416 : make_not_null(&subdomain_solve_initial_guess_in_solution_out),
417 : subdomain_operator, subdomain_residual, std::forward_as_tuple(box));
418 : // Re-naming the solution buffer for the code below
419 : auto& subdomain_solution = subdomain_solve_initial_guess_in_solution_out;
420 :
421 : // Do some logging and observing
422 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
423 : ::Verbosity::Quiet)) {
424 : if (not subdomain_solve_has_converged or
425 : subdomain_solve_has_converged.reason() ==
426 : Convergence::Reason::MaxIterations) {
427 : Parallel::printf(
428 : "%s %s(%zu): WARNING: Subdomain solver did not converge in %zu "
429 : "iterations: %e -> %e\n",
430 : element_id, pretty_type::name<OptionsGroup>(), iteration_id,
431 : subdomain_solve_has_converged.num_iterations(),
432 : subdomain_solve_has_converged.initial_residual_magnitude(),
433 : subdomain_solve_has_converged.residual_magnitude());
434 : } else if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
435 : ::Verbosity::Debug)) {
436 : Parallel::printf(
437 : "%s %s(%zu): Subdomain solver converged in %zu iterations (%s): %e "
438 : "-> %e\n",
439 : element_id, pretty_type::name<OptionsGroup>(), iteration_id,
440 : subdomain_solve_has_converged.num_iterations(),
441 : subdomain_solve_has_converged.reason(),
442 : subdomain_solve_has_converged.initial_residual_magnitude(),
443 : subdomain_solve_has_converged.residual_magnitude());
444 : }
445 : }
446 : const std::optional<std::string> section_observation_key =
447 : observers::get_section_observation_key<ArraySectionIdTag>(box);
448 : if (section_observation_key.has_value()) {
449 : contribute_to_subdomain_stats_observation<OptionsGroup,
450 : ParallelComponent>(
451 : iteration_id + 1, subdomain_solve_has_converged.num_iterations(),
452 : cache, element_id, *section_observation_key,
453 : db::get<Tags::ObservePerCoreReductions<OptionsGroup>>(box));
454 : }
455 :
456 : // Apply weighting
457 : if (LIKELY(max_overlap > 0)) {
458 : subdomain_solution.element_data *=
459 : get(db::get<Tags::Weight<OptionsGroup>>(box));
460 : }
461 :
462 : // Apply solution to central element
463 : db::mutate<fields_tag>(
464 : [&subdomain_solution](const auto fields) {
465 : *fields += subdomain_solution.element_data;
466 : },
467 : make_not_null(&box));
468 :
469 : // Send overlap solutions back to the neighbors that they are on
470 : if (LIKELY(max_overlap > 0)) {
471 : auto& receiver_proxy =
472 : Parallel::get_parallel_component<ParallelComponent>(cache);
473 : for (auto& [overlap_id, overlap_solution] :
474 : subdomain_solution.overlap_data) {
475 : const auto& direction = overlap_id.direction();
476 : const auto& neighbor_id = overlap_id.id();
477 : const auto& orientation =
478 : element.neighbors().at(direction).orientation();
479 : const auto direction_from_neighbor = orientation(direction.opposite());
480 : Parallel::receive_data<overlap_solution_inbox_tag>(
481 : receiver_proxy[neighbor_id], iteration_id,
482 : std::make_pair(
483 : OverlapId<Dim>{direction_from_neighbor, element.id()},
484 : std::move(overlap_solution)));
485 : }
486 : }
487 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
488 : }
489 : };
490 :
491 : // Wait for the subdomain solutions on regions within this element that overlap
492 : // with neighboring element-centered subdomains. Combine the solutions as a
493 : // weighted sum.
494 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator>
495 : struct ReceiveOverlapSolution {
496 : private:
497 : using fields_tag = FieldsTag;
498 : using residual_tag =
499 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
500 : static constexpr size_t Dim = SubdomainOperator::volume_dim;
501 : using SubdomainData =
502 : ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
503 : using OverlapSolution = typename SubdomainData::OverlapData;
504 : using overlap_solution_inbox_tag =
505 : OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapSolution>;
506 :
507 : public:
508 : using const_global_cache_tags = tmpl::list<Tags::MaxOverlap<OptionsGroup>>;
509 : using inbox_tags = tmpl::list<overlap_solution_inbox_tag>;
510 :
511 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
512 : typename ActionList, typename ParallelComponent>
513 : static Parallel::iterable_action_return_t apply(
514 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
515 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
516 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
517 : const ParallelComponent* const /*meta*/) {
518 : const size_t iteration_id =
519 : get<Convergence::Tags::IterationId<OptionsGroup>>(box);
520 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
521 :
522 : // Nothing to do if overlap is empty
523 : if (UNLIKELY(db::get<Tags::MaxOverlap<OptionsGroup>>(box) == 0 or
524 : element.number_of_neighbors() == 0)) {
525 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
526 : }
527 :
528 : if (not dg::has_received_from_all_mortars<overlap_solution_inbox_tag>(
529 : iteration_id, element, inboxes)) {
530 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
531 : }
532 :
533 : // Do some logging
534 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
535 : ::Verbosity::Debug)) {
536 : Parallel::printf("%s %s(%zu): Receive overlap solution\n", element_id,
537 : pretty_type::name<OptionsGroup>(), iteration_id);
538 : }
539 :
540 : // Add solutions on overlaps to this element's solution in a weighted sum
541 : const auto received_overlap_solutions =
542 : std::move(tuples::get<overlap_solution_inbox_tag>(inboxes)
543 : .extract(iteration_id)
544 : .mapped());
545 : db::mutate<fields_tag>(
546 : [&received_overlap_solutions](
547 : const auto fields, const Index<Dim>& full_extents,
548 : const std::array<size_t, Dim>& all_intruding_extents,
549 : const DirectionMap<Dim, Scalar<DataVector>>&
550 : all_intruding_overlap_weights) {
551 : for (const auto& [overlap_id, overlap_solution] :
552 : received_overlap_solutions) {
553 : const auto& direction = overlap_id.direction();
554 : const auto& intruding_extents =
555 : gsl::at(all_intruding_extents, direction.dimension());
556 : const auto& overlap_weight =
557 : all_intruding_overlap_weights.at(direction);
558 : LinearSolver::Schwarz::add_overlap_data(
559 : fields, overlap_solution * get(overlap_weight), full_extents,
560 : intruding_extents, direction);
561 : }
562 : },
563 : make_not_null(&box), db::get<domain::Tags::Mesh<Dim>>(box).extents(),
564 : db::get<Tags::IntrudingExtents<Dim, OptionsGroup>>(box),
565 : db::get<domain::Tags::Faces<Dim, Tags::Weight<OptionsGroup>>>(box));
566 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
567 : }
568 : };
569 :
570 : } // namespace LinearSolver::Schwarz::detail
|
