15 #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 #include "DataStructures/DataBox/Tag.hpp"
18 #include "Domain/InterfaceComputeTags.hpp"
19 #include "Domain/InterfaceHelpers.hpp"
22 #include "Domain/Structure/OrientationMapHelpers.hpp"
24 #include "IO/Observer/Actions/RegisterWithObservers.hpp"
25 #include "IO/Observer/ArrayComponentId.hpp"
26 #include "IO/Observer/ObservationId.hpp"
27 #include "IO/Observer/ReductionActions.hpp"
28 #include "IO/Observer/TypeOfObservation.hpp"
29 #include "Informer/Tags.hpp"
30 #include "Informer/Verbosity.hpp"
31 #include "NumericalAlgorithms/Convergence/Tags.hpp"
32 #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
37 #include "Parallel/InboxInserters.hpp"
38 #include "Parallel/Invoke.hpp"
39 #include "Parallel/ParallelComponentHelpers.hpp"
41 #include "Parallel/Reduction.hpp"
42 #include "ParallelAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
43 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
44 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
46 #include "ParallelAlgorithms/LinearSolver/Schwarz/ComputeTags.hpp"
47 #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
48 #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
49 #include "ParallelAlgorithms/LinearSolver/Schwarz/Protocols.hpp"
50 #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
54 #include "Utilities/ProtocolHelpers.hpp"
56 #include "Utilities/TaggedTuple.hpp"
58 namespace LinearSolver::Schwarz::detail {
60 using reduction_data = Parallel::ReductionData<
73 template <
typename OptionsGroup>
74 struct RegisterObservers {
75 template <
typename ParallelComponent,
typename DbTagsList,
78 register_info(
const db::DataBox<DbTagsList>& ,
79 const ArrayIndex& ) noexcept {
86 template <
typename FieldsTag,
typename OptionsGroup,
typename SourceTag>
87 using RegisterElement =
90 template <
typename OptionsGroup,
typename ParallelComponent,
91 typename Metavariables,
typename ArrayIndex>
92 void contribute_to_subdomain_stats_observation(
93 const size_t iteration_id,
const size_t subdomain_solve_num_iterations,
95 const ArrayIndex& array_index) noexcept {
96 auto& local_observer =
97 *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
100 Parallel::simple_action<observers::Actions::ContributeReductionData>(
104 pretty_type::get_name<OptionsGroup>() +
"SubdomainSolves"),
108 std::string{
"/" + Options::name<OptionsGroup>() +
"SubdomainSolves"},
110 "MinNumIterations",
"MaxNumIterations"},
111 reduction_data{iteration_id, 1, subdomain_solve_num_iterations,
112 subdomain_solve_num_iterations,
113 subdomain_solve_num_iterations});
116 template <
typename SubdomainDataType,
typename OptionsGroup>
119 return "SubdomainData(" + Options::name<OptionsGroup>() +
")";
121 using type = SubdomainDataType;
124 template <
typename FieldsTag,
typename OptionsGroup,
typename SubdomainOperator>
125 struct InitializeElement {
127 using fields_tag = FieldsTag;
130 static constexpr
size_t Dim = SubdomainOperator::volume_dim;
131 using SubdomainData =
132 ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
136 using subdomain_solver_tag =
137 Tags::SubdomainSolver<LinearSolver::Serial::Gmres<SubdomainData>,
141 using initialization_tags =
142 tmpl::list<domain::Tags::InitialExtents<Dim>, subdomain_solver_tag>;
143 using initialization_tags_to_keep = tmpl::list<subdomain_solver_tag>;
144 using const_global_cache_tags = tmpl::list<Tags::MaxOverlap<OptionsGroup>>;
147 tmpl::list<Tags::Overlaps<residual_tag, Dim, OptionsGroup>,
148 SubdomainDataBufferTag<SubdomainData, OptionsGroup>>;
149 using compute_tags = tmpl::list<
157 Tags::IntrudingExtentsCompute<Dim, OptionsGroup>,
158 Tags::IntrudingOverlapWidthsCompute<Dim, OptionsGroup>,
160 Tags::ElementWeightCompute<Dim, OptionsGroup>,
163 Tags::IntrudingOverlapWeightCompute<Dim, OptionsGroup>>>;
164 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
165 typename ActionList,
typename ParallelComponent>
166 static auto apply(db::DataBox<DbTagsList>& box,
171 const ParallelComponent*
const ) noexcept {
172 const auto& element_mesh = db::get<domain::Tags::Mesh<Dim>>(box);
173 const size_t element_num_points = element_mesh.number_of_grid_points();
175 tmpl::list<SubdomainDataBufferTag<SubdomainData, OptionsGroup>>>(
177 return std::make_tuple(std::move(box));
183 template <
typename FieldsTag,
typename OptionsGroup,
typename SubdomainOperator>
185 tmpl::list<db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>>,
190 template <
typename FieldsTag,
typename OptionsGroup,
typename SubdomainOperator>
192 SubdomainOperator::volume_dim,
193 tmpl::list<db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>>,
196 template <
size_t Dim,
typename OptionsGroup,
typename OverlapSolution>
197 struct OverlapSolutionInboxTag
199 OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapSolution>> {
200 using temporal_id = size_t;
208 template <
typename FieldsTag,
typename OptionsGroup,
typename SubdomainOperator>
209 struct SolveSubdomain {
211 using fields_tag = FieldsTag;
214 static constexpr
size_t Dim = SubdomainOperator::volume_dim;
215 using SubdomainData =
216 ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
217 using OverlapData =
typename SubdomainData::OverlapData;
218 using overlap_solution_inbox_tag =
219 OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapData>;
222 using const_global_cache_tags =
223 tmpl::list<Tags::MaxOverlap<OptionsGroup>,
226 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
227 typename ActionList,
typename ParallelComponent>
229 db::DataBox<DbTagsList>& box,
233 const ParallelComponent*
const ) noexcept {
234 const size_t iteration_id =
235 get<Convergence::Tags::IterationId<OptionsGroup>>(box);
239 ::Verbosity::Debug)) {
241 "%s " + Options::name<OptionsGroup>() +
"(%zu): Solve subdomain\n",
242 element_id, iteration_id);
245 const auto& element = db::get<domain::Tags::Element<Dim>>(box);
246 const size_t max_overlap = db::get<Tags::MaxOverlap<OptionsGroup>>(box);
250 db::mutate<SubdomainDataBufferTag<SubdomainData, OptionsGroup>,
251 Tags::Overlaps<residual_tag, Dim, OptionsGroup>>(
253 [max_overlap, &element](
255 const auto overlap_residuals,
const auto& residual) noexcept {
256 subdomain_data->element_data = residual;
258 if (
LIKELY(max_overlap > 0 and element.number_of_neighbors() > 0)) {
259 subdomain_data->overlap_data = std::move(*overlap_residuals);
262 db::get<residual_tag>(box));
263 const auto& subdomain_residual =
264 db::get<SubdomainDataBufferTag<SubdomainData, OptionsGroup>>(box);
267 const size_t num_points =
268 db::get<domain::Tags::Mesh<Dim>>(box).number_of_grid_points();
269 SubdomainOperator subdomain_operator{num_points};
272 const auto apply_subdomain_operator =
274 const SubdomainData& operand) noexcept {
278 db::apply<typename SubdomainOperator::element_operator>(
280 tmpl::for_each<tmpl::list<domain::Tags::InternalDirections<Dim>,
282 [&box, &operand, &result,
283 &subdomain_operator](
auto directions_v) noexcept {
284 using directions = tmpl::type_from<decltype(directions_v)>;
285 using face_operator =
286 typename SubdomainOperator::template face_operator<directions>;
287 interface_apply<directions, face_operator>(
293 const auto& subdomain_solver =
294 get<Tags::SubdomainSolverBase<OptionsGroup>>(box);
295 auto subdomain_solve_initial_guess_in_solution_out =
296 make_with_value<SubdomainData>(subdomain_residual, 0.);
297 const auto subdomain_solve_has_converged = subdomain_solver.solve(
298 make_not_null(&subdomain_solve_initial_guess_in_solution_out),
299 apply_subdomain_operator, subdomain_residual);
301 auto& subdomain_solution = subdomain_solve_initial_guess_in_solution_out;
305 ::Verbosity::Quiet)) {
306 if (not subdomain_solve_has_converged or
307 subdomain_solve_has_converged.reason() ==
310 "%s WARNING: Subdomain solver did not converge in %zu iterations: "
312 element_id, subdomain_solve_has_converged.num_iterations(),
313 subdomain_solve_has_converged.initial_residual_magnitude(),
314 subdomain_solve_has_converged.residual_magnitude());
316 ::Verbosity::Debug)) {
318 "%s Subdomain solver converged in %zu iterations (%s): %e -> %e\n",
319 element_id, subdomain_solve_has_converged.num_iterations(),
320 subdomain_solve_has_converged.reason(),
321 subdomain_solve_has_converged.initial_residual_magnitude(),
322 subdomain_solve_has_converged.residual_magnitude());
325 contribute_to_subdomain_stats_observation<OptionsGroup, ParallelComponent>(
326 iteration_id + 1, subdomain_solve_has_converged.num_iterations(), cache,
330 if (
LIKELY(max_overlap > 0)) {
331 subdomain_solution.element_data *=
332 get(
db::get<Tags::Weight<OptionsGroup>>(box));
337 [&subdomain_solution](
const auto fields) noexcept {
338 *fields += subdomain_solution.element_data;
342 if (
LIKELY(max_overlap > 0)) {
343 auto& receiver_proxy =
344 Parallel::get_parallel_component<ParallelComponent>(cache);
345 for (
auto& [overlap_id, overlap_solution] :
346 subdomain_solution.overlap_data) {
347 const auto& direction = overlap_id.first;
348 const auto& neighbor_id = overlap_id.second;
349 const auto& orientation =
350 element.neighbors().at(direction).orientation();
351 const auto direction_from_neighbor = orientation(direction.opposite());
352 Parallel::receive_data<overlap_solution_inbox_tag>(
353 receiver_proxy[neighbor_id], iteration_id,
355 OverlapId<Dim>{direction_from_neighbor, element.id()},
356 std::move(overlap_solution)));
359 return {std::move(box)};
366 template <
typename FieldsTag,
typename OptionsGroup,
typename SubdomainOperator>
367 struct ReceiveOverlapSolution {
369 using fields_tag = FieldsTag;
372 static constexpr
size_t Dim = SubdomainOperator::volume_dim;
373 using SubdomainData =
374 ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
375 using OverlapSolution =
typename SubdomainData::OverlapData;
376 using overlap_solution_inbox_tag =
377 OverlapSolutionInboxTag<Dim, OptionsGroup, OverlapSolution>;
380 using const_global_cache_tags = tmpl::list<Tags::MaxOverlap<OptionsGroup>>;
381 using inbox_tags = tmpl::list<overlap_solution_inbox_tag>;
383 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
385 static bool is_ready(
const db::DataBox<DbTagsList>& box,
392 return dg::has_received_from_all_mortars<overlap_solution_inbox_tag>(
397 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
398 typename ActionList,
typename ParallelComponent>
403 const ParallelComponent*
const ) noexcept {
404 const size_t iteration_id =
405 get<Convergence::Tags::IterationId<OptionsGroup>>(box);
406 const auto& element = db::get<domain::Tags::Element<Dim>>(box);
410 element.number_of_neighbors() == 0)) {
411 return {std::move(box)};
416 ::Verbosity::Debug)) {
418 "(%zu): Receive overlap solution\n",
419 element_id, iteration_id);
423 const auto received_overlap_solutions =
424 std::move(tuples::get<overlap_solution_inbox_tag>(inboxes)
425 .extract(iteration_id)
427 db::mutate<fields_tag>(
429 [&received_overlap_solutions](
430 const auto fields,
const Index<Dim>& full_extents,
433 all_intruding_overlap_weights) noexcept {
434 for (
const auto& [overlap_id, overlap_solution] :
435 received_overlap_solutions) {
436 const auto& direction = overlap_id.first;
437 const auto& intruding_extents =
438 gsl::at(all_intruding_extents, direction.dimension());
439 const auto& overlap_weight =
440 all_intruding_overlap_weights.at(direction);
442 fields, overlap_solution *
get(overlap_weight), full_extents,
443 intruding_extents, direction);
446 db::get<domain::Tags::Mesh<Dim>>(box).extents(),
447 db::get<Tags::IntrudingExtents<Dim, OptionsGroup>>(box),
449 Tags::Weight<OptionsGroup>>>(box));
450 return {std::move(box)};