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