SpECTRE Documentation Coverage Report
Current view: top level - ParallelAlgorithms/LinearSolver/Schwarz - ElementActions.hpp Hit Total Coverage
Commit: 3c2e9d3ed337bca2146eee9de07432e292a38c3a Lines: 0 1 0.0 %
Date: 2024-06-11 22:56:19
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14