SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/DiscontinuousGalerkin/SubdomainOperator - SubdomainOperator.hpp Hit Total Coverage
Commit: 3f09028930c0450a2fb61ee918b22882f5d03d2b Lines: 2 33 6.1 %
Date: 2021-10-22 20:52:16
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 <optional>
       8             : #include <ostream>
       9             : #include <pup.h>
      10             : #include <tuple>
      11             : #include <utility>
      12             : 
      13             : #include "DataStructures/DataBox/DataBox.hpp"
      14             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      15             : #include "DataStructures/Variables.hpp"
      16             : #include "Domain/Domain.hpp"
      17             : #include "Domain/FaceNormal.hpp"
      18             : #include "Domain/Structure/Direction.hpp"
      19             : #include "Domain/Structure/DirectionMap.hpp"
      20             : #include "Domain/Structure/ElementId.hpp"
      21             : #include "Domain/Tags.hpp"
      22             : #include "Domain/Tags/FaceNormal.hpp"
      23             : #include "Domain/Tags/Faces.hpp"
      24             : #include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp"
      25             : #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
      26             : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
      27             : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
      28             : #include "Elliptic/Systems/GetSourcesComputer.hpp"
      29             : #include "Elliptic/Utilities/ApplyAt.hpp"
      30             : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
      31             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
      32             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      33             : #include "Parallel/PupStlCpp17.hpp"
      34             : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
      35             : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
      36             : #include "ParallelAlgorithms/LinearSolver/Schwarz/SubdomainOperator.hpp"
      37             : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
      38             : #include "Utilities/ErrorHandling/Assert.hpp"
      39             : #include "Utilities/Gsl.hpp"
      40             : #include "Utilities/TMPL.hpp"
      41             : 
      42             : /// Items related to the restriction of the DG operator to an element-centered
      43             : /// subdomain
      44             : namespace elliptic::dg::subdomain_operator {
      45             : 
      46             : namespace detail {
      47             : // Wrap the `Tag` in `LinearSolver::Schwarz::Tags::Overlaps`, except if it is
      48             : // included in `TakeFromCenterTags`.
      49             : template <typename Tag, typename Dim, typename OptionsGroup,
      50             :           typename TakeFromCenterTags>
      51             : struct make_overlap_tag_impl {
      52             :   using type = tmpl::conditional_t<
      53             :       tmpl::list_contains_v<TakeFromCenterTags, Tag>, Tag,
      54             :       LinearSolver::Schwarz::Tags::Overlaps<Tag, Dim::value, OptionsGroup>>;
      55             : };
      56             : 
      57             : // Wrap the `Tag` in `Tags::NeighborMortars`
      58             : template <typename Tag, typename Dim>
      59             : struct make_neighbor_mortars_tag_impl {
      60             :   using type = Tags::NeighborMortars<Tag, Dim::value>;
      61             : };
      62             : }  // namespace detail
      63             : 
      64             : /*!
      65             :  * \brief The elliptic DG operator on an element-centered subdomain
      66             :  *
      67             :  * This operator is a restriction of the full (linearized) DG-operator to an
      68             :  * element-centered subdomain with a few points overlap into neighboring
      69             :  * elements. It is a `LinearSolver::Schwarz::SubdomainOperator` to be used with
      70             :  * the Schwarz linear solver when it solves the elliptic DG operator.
      71             :  *
      72             :  * This operator requires the following tags are available on overlap regions
      73             :  * with neighboring elements:
      74             :  *
      75             :  * - Geometric quantities provided by
      76             :  *   `elliptic::dg::subdomain_operator::InitializeSubdomain`.
      77             :  * - All `System::fluxes_computer::argument_tags` and
      78             :  *   `System::sources_computer::argument_tags` (or
      79             :  *   `System::sources_computer_linearized::argument_tags` for nonlinear
      80             :  *   systems), except those listed in `ArgsTagsFromCenter`. The latter will be
      81             :  *   taken from the central element's DataBox, so they don't need to be made
      82             :  *   available on overlaps.
      83             :  * - The `System::fluxes_computer::argument_tags` on internal and external
      84             :  *   interfaces, except those listed in `System::fluxes_computer::volume_tags`.
      85             :  *
      86             :  * Some of these tags may require communication between elements. For example,
      87             :  * nonlinear system fields are constant background fields for the linearized DG
      88             :  * operator, but are updated in every nonlinear solver iteration. Therefore, the
      89             :  * updated nonlinear fields must be communicated across overlaps between
      90             :  * nonlinear solver iterations. To perform the communication you can use
      91             :  * `LinearSolver::Schwarz::Actions::SendOverlapFields` and
      92             :  * `LinearSolver::Schwarz::Actions::ReceiveOverlapFields`, setting
      93             :  * `RestrictToOverlap` to `false`. See
      94             :  * `LinearSolver::Schwarz::SubdomainOperator` for details.
      95             :  *
      96             :  * \par Overriding boundary conditions
      97             :  * Sometimes the subdomain operator should not use the boundary conditions that
      98             :  * have been selected when setting up the domain. For example, when the
      99             :  * subdomain operator is cached between non-linear solver iterations but the
     100             :  * boundary conditions depend on the non-linear fields, the preconditioning can
     101             :  * become ineffective (see
     102             :  * `LinearSolver::Schwarz::Actions::ResetSubdomainSolver`). Another example is
     103             :  * `elliptic::subdomain_preconditioners::MinusLaplacian`, where an auxiliary
     104             :  * Poisson system is used for preconditioning that doesn't have boundary
     105             :  * conditions set up in the domain. In these cases, the boundary conditions used
     106             :  * for the subdomain operator can be overridden with the
     107             :  * `override_boundary_conditions` argument in the constructor. Note that the
     108             :  * subdomain operator always applies the _linearized_ boundary conditions.
     109             :  *
     110             :  * \warning The subdomain operator hasn't been tested with periodic boundary
     111             :  * conditions so far.
     112             :  */
     113             : template <typename System, typename OptionsGroup,
     114             :           typename ArgsTagsFromCenter = tmpl::list<>>
     115           1 : struct SubdomainOperator
     116             :     : LinearSolver::Schwarz::SubdomainOperator<System::volume_dim> {
     117             :  private:
     118           0 :   static constexpr size_t Dim = System::volume_dim;
     119             : 
     120             :   // Operator applications happen sequentially so we don't have to keep track of
     121             :   // the temporal id
     122           0 :   static constexpr size_t temporal_id = 0;
     123             : 
     124             :   // The subdomain operator always applies the linearized DG operator
     125           0 :   static constexpr bool linearized = true;
     126             : 
     127           0 :   using BoundaryConditionsBase = typename System::boundary_conditions_base;
     128             : 
     129             :   // These are the arguments that we need to retrieve from the DataBox and pass
     130             :   // to the functions in `elliptic::dg`, both on the central element and on
     131             :   // neighbors
     132           0 :   using prepare_args_tags =
     133             :       tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>,
     134             :                  domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     135             :                                                Frame::Inertial>,
     136             :                  domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>,
     137             :                  domain::Tags::Faces<
     138             :                      Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>,
     139             :                  ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
     140             :                  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>;
     141           0 :   using apply_args_tags = tmpl::list<
     142             :       domain::Tags::Mesh<Dim>,
     143             :       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     144             :                                     Frame::Inertial>,
     145             :       domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
     146             :       domain::Tags::Faces<Dim,
     147             :                           domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>,
     148             :       ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
     149             :       ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
     150             :       elliptic::dg::Tags::PenaltyParameter, elliptic::dg::Tags::Massive>;
     151           0 :   using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
     152           0 :   using sources_args_tags =
     153             :       typename elliptic::get_sources_computer<System,
     154             :                                               linearized>::argument_tags;
     155             : 
     156             :   // We need the fluxes args also on interfaces (internal and external). The
     157             :   // volume tags are the subset that don't have to be taken from interfaces.
     158           0 :   using fluxes_args_volume_tags = typename System::fluxes_computer::volume_tags;
     159             : 
     160             :   // These tags can be taken directly from the central element's DataBox, even
     161             :   // when evaluating neighbors
     162           0 :   using args_tags_from_center = tmpl::remove_duplicates<
     163             :       tmpl::push_back<ArgsTagsFromCenter, elliptic::dg::Tags::PenaltyParameter,
     164             :                       elliptic::dg::Tags::Massive>>;
     165             : 
     166             :   // Data on neighbors is stored in the central element's DataBox in
     167             :   // `LinearSolver::Schwarz::Tags::Overlaps` maps, so we wrap the argument tags
     168             :   // with this prefix
     169           0 :   using make_overlap_tag =
     170             :       detail::make_overlap_tag_impl<tmpl::_1, tmpl::pin<tmpl::size_t<Dim>>,
     171             :                                     tmpl::pin<OptionsGroup>,
     172             :                                     tmpl::pin<args_tags_from_center>>;
     173           0 :   using prepare_args_tags_overlap =
     174             :       tmpl::transform<prepare_args_tags, make_overlap_tag>;
     175           0 :   using apply_args_tags_overlap =
     176             :       tmpl::transform<apply_args_tags, make_overlap_tag>;
     177           0 :   using fluxes_args_tags_overlap =
     178             :       tmpl::transform<fluxes_args_tags, make_overlap_tag>;
     179           0 :   using sources_args_tags_overlap =
     180             :       tmpl::transform<sources_args_tags, make_overlap_tag>;
     181           0 :   using fluxes_args_tags_overlap_faces = tmpl::transform<
     182             :       domain::make_faces_tags<Dim, fluxes_args_tags, fluxes_args_volume_tags>,
     183             :       make_overlap_tag>;
     184             : 
     185             :   // We also need some data on the remote side of all neighbors' mortars. Such
     186             :   // data is stored in the central element's DataBox in `Tags::NeighborMortars`
     187             :   // maps
     188           0 :   using make_neighbor_mortars_tag =
     189             :       detail::make_neighbor_mortars_tag_impl<tmpl::_1,
     190             :                                              tmpl::pin<tmpl::size_t<Dim>>>;
     191             : 
     192             :  public:
     193           0 :   SubdomainOperator(std::optional<std::unique_ptr<BoundaryConditionsBase>>
     194             :                         override_boundary_conditions = std::nullopt)
     195             :       : override_boundary_conditions_(std::move(override_boundary_conditions)) {
     196             :   }
     197             : 
     198             :   /// \warning This function is not thread-safe because it accesses mutable
     199             :   /// memory buffers.
     200             :   template <typename ResultTags, typename OperandTags, typename DbTagsList>
     201           1 :   void operator()(
     202             :       const gsl::not_null<
     203             :           LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, ResultTags>*>
     204             :           result,
     205             :       const LinearSolver::Schwarz::ElementCenteredSubdomainData<
     206             :           Dim, OperandTags>& operand,
     207             :       const db::DataBox<DbTagsList>& box) const {
     208             :     // Used to retrieve items out of the DataBox to forward to functions. This
     209             :     // replaces a long series of db::get calls.
     210             :     const auto get_items = [](const auto&... args) {
     211             :       return std::forward_as_tuple(args...);
     212             :     };
     213             : 
     214             :     // Retrieve data out of the DataBox
     215             :     using tags_to_retrieve = tmpl::flatten<tmpl::list<
     216             :         domain::Tags::Domain<Dim>, domain::Tags::Element<Dim>,
     217             :         ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
     218             :         // Data on overlaps with neighbors
     219             :         tmpl::transform<
     220             :             tmpl::flatten<tmpl::list<
     221             :                 Tags::ExtrudingExtent, domain::Tags::Element<Dim>,
     222             :                 domain::Tags::Mesh<Dim>,
     223             :                 domain::Tags::Faces<
     224             :                     Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>,
     225             :                 ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
     226             :                 ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
     227             :                 // Data on the remote side of the neighbor's mortars
     228             :                 tmpl::transform<
     229             :                     tmpl::list<
     230             :                         domain::Tags::Mesh<Dim>,
     231             :                         domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>,
     232             :                         domain::Tags::Mesh<Dim - 1>,
     233             :                         ::Tags::MortarSize<Dim - 1>>,
     234             :                     make_neighbor_mortars_tag>>>,
     235             :             make_overlap_tag>>>;
     236             :     const auto& [domain, central_element, central_mortar_meshes,
     237             :                  all_overlap_extents, all_neighbor_elements,
     238             :                  all_neighbor_meshes,
     239             :                  all_neighbor_face_normal_magnitudes,
     240             :                  all_neighbor_mortar_meshes, all_neighbor_mortar_sizes,
     241             :                  all_neighbors_neighbor_meshes,
     242             :                  all_neighbors_neighbor_face_normal_magnitudes,
     243             :                  all_neighbors_neighbor_mortar_meshes,
     244             :                  all_neighbors_neighbor_mortar_sizes] =
     245             :         db::apply<tags_to_retrieve>(get_items, box);
     246             :     const auto fluxes_args = db::apply<fluxes_args_tags>(get_items, box);
     247             :     const auto sources_args = db::apply<sources_args_tags>(get_items, box);
     248             :     using FluxesArgs = std::decay_t<decltype(fluxes_args)>;
     249             :     DirectionMap<Dim, FluxesArgs> fluxes_args_on_faces{};
     250             :     for (const auto& direction : Direction<Dim>::all_directions()) {
     251             :       fluxes_args_on_faces.emplace(
     252             :           direction, elliptic::util::apply_at<
     253             :                          domain::make_faces_tags<Dim, fluxes_args_tags,
     254             :                                                  fluxes_args_volume_tags>,
     255             :                          fluxes_args_volume_tags>(get_items, box, direction));
     256             :     }
     257             : 
     258             :     // Setup boundary conditions
     259             :     const auto apply_boundary_condition =
     260             :         [&box, &local_domain = domain,
     261             :          &override_boundary_conditions = override_boundary_conditions_](
     262             :             const ElementId<Dim>& local_element_id,
     263             :             const Direction<Dim>& local_direction, auto is_overlap,
     264             :             const auto& map_keys, const auto... fields_and_fluxes) {
     265             :           constexpr bool is_overlap_v =
     266             :               std::decay_t<decltype(is_overlap)>::value;
     267             :           // Get boundary conditions from domain, or use overridden boundary
     268             :           // conditions
     269             :           const auto& boundary_condition = [&local_domain, &local_element_id,
     270             :                                             &local_direction,
     271             :                                             &override_boundary_conditions]()
     272             :               -> const BoundaryConditionsBase& {
     273             :             if (override_boundary_conditions.has_value()) {
     274             :               return **override_boundary_conditions;
     275             :             }
     276             :             const auto& boundary_conditions =
     277             :                 local_domain.blocks()
     278             :                     .at(local_element_id.block_id())
     279             :                     .external_boundary_conditions();
     280             :             ASSERT(boundary_conditions.contains(local_direction),
     281             :                    "No boundary condition is available in block "
     282             :                        << local_element_id.block_id() << " in direction "
     283             :                        << local_direction
     284             :                        << ". Make sure you are setting up boundary conditions "
     285             :                           "when creating the domain.");
     286             :             ASSERT(
     287             :                 dynamic_cast<const BoundaryConditionsBase*>(
     288             :                     boundary_conditions.at(local_direction).get()) != nullptr,
     289             :                 "The boundary condition in block "
     290             :                     << local_element_id.block_id() << " in direction "
     291             :                     << local_direction
     292             :                     << " has an unexpected type. Make sure it derives off the "
     293             :                        "'boundary_conditions_base' class set in the system.");
     294             :             return dynamic_cast<const BoundaryConditionsBase&>(
     295             :                 *boundary_conditions.at(local_direction));
     296             :           }();
     297             :           elliptic::apply_boundary_condition<
     298             :               linearized,
     299             :               tmpl::conditional_t<is_overlap_v, make_overlap_tag, void>>(
     300             :               boundary_condition, box, map_keys, fields_and_fluxes...);
     301             :         };
     302             : 
     303             :     // The subdomain operator essentially does two sweeps over all elements in
     304             :     // the subdomain: In the first sweep it prepares the mortar data and stores
     305             :     // them on both sides of all mortars, and in the second sweep it consumes
     306             :     // the mortar data to apply the operator. This implementation is relatively
     307             :     // simple because it can re-use the implementation for the parallel DG
     308             :     // operator. However, it is also possible to apply the subdomain operator in
     309             :     // a single sweep over all elements, incrementally building up the mortar
     310             :     // data and applying boundary corrections immediately to both adjacent
     311             :     // elements once the data is available. That approach is possibly a
     312             :     // performance optimization but requires re-implementing a lot of logic for
     313             :     // the DG operator here. It should be considered once the subdomain operator
     314             :     // has been identified as the performance bottleneck. An alternative to
     315             :     // optimizing the subdomain operator performance is to precondition the
     316             :     // subdomain solve with a _much_ simpler subdomain operator, such as a
     317             :     // finite-difference Laplacian, so fewer applications of the more expensive
     318             :     // DG subdomain operator are necessary.
     319             : 
     320             :     // 1. Prepare mortar data on all elements in the subdomain and store them on
     321             :     //    mortars, reorienting if needed
     322             :     //
     323             :     // Prepare central element
     324             :     const auto apply_boundary_condition_center =
     325             :         [&apply_boundary_condition, &local_central_element = central_element](
     326             :             const Direction<Dim>& local_direction,
     327             :             const auto... fields_and_fluxes) {
     328             :           apply_boundary_condition(local_central_element.id(), local_direction,
     329             :                                    std::false_type{}, local_direction,
     330             :                                    fields_and_fluxes...);
     331             :         };
     332             :     db::apply<prepare_args_tags>(
     333             :         [this, &operand](const auto&... args) {
     334             :           elliptic::dg::prepare_mortar_data<System, linearized>(
     335             :               make_not_null(&central_auxiliary_vars_),
     336             :               make_not_null(&central_auxiliary_fluxes_),
     337             :               make_not_null(&central_primal_fluxes_),
     338             :               make_not_null(&central_mortar_data_), operand.element_data,
     339             :               args...);
     340             :         },
     341             :         box, temporal_id, apply_boundary_condition_center, fluxes_args,
     342             :         sources_args, fluxes_args_on_faces);
     343             :     // Prepare neighbors
     344             :     for (const auto& [direction, neighbors] : central_element.neighbors()) {
     345             :       const auto& orientation = neighbors.orientation();
     346             :       const auto direction_from_neighbor = orientation(direction.opposite());
     347             :       for (const auto& neighbor_id : neighbors) {
     348             :         const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
     349             :                                                                neighbor_id};
     350             :         const auto& overlap_extent = all_overlap_extents.at(overlap_id);
     351             :         const auto& neighbor = all_neighbor_elements.at(overlap_id);
     352             :         const auto& neighbor_mesh = all_neighbor_meshes.at(overlap_id);
     353             :         const auto& mortar_id = overlap_id;
     354             :         const auto& mortar_mesh = central_mortar_meshes.at(mortar_id);
     355             :         const ::dg::MortarId<Dim> mortar_id_from_neighbor{
     356             :             direction_from_neighbor, central_element.id()};
     357             : 
     358             :         // Intercept empty overlaps. In the unlikely case that overlaps have
     359             :         // zero extent, meaning no point of the neighbor is part of the
     360             :         // subdomain (which is fairly useless, except for testing), the
     361             :         // subdomain is identical to the central element and no communication
     362             :         // with neighbors is necessary. We can just handle the mortar between
     363             :         // central element and neighbor and continue.
     364             :         if (UNLIKELY(overlap_extent == 0)) {
     365             :           const auto& mortar_mesh_from_neighbor =
     366             :               all_neighbor_mortar_meshes.at(overlap_id)
     367             :                   .at(mortar_id_from_neighbor);
     368             :           const auto& mortar_size_from_neighbor =
     369             :               all_neighbor_mortar_sizes.at(overlap_id)
     370             :                   .at(mortar_id_from_neighbor);
     371             :           auto remote_boundary_data =
     372             :               elliptic::dg::zero_boundary_data_on_mortar<
     373             :                   typename System::primal_fields,
     374             :                   typename System::primal_fluxes>(
     375             :                   direction_from_neighbor, neighbor_mesh,
     376             :                   all_neighbor_face_normal_magnitudes.at(overlap_id)
     377             :                       .at(direction_from_neighbor),
     378             :                   mortar_mesh_from_neighbor, mortar_size_from_neighbor);
     379             :           if (not orientation.is_aligned()) {
     380             :             remote_boundary_data.orient_on_slice(
     381             :                 mortar_mesh_from_neighbor.extents(),
     382             :                 direction_from_neighbor.dimension(), orientation.inverse_map());
     383             :           }
     384             :           central_mortar_data_.at(mortar_id).remote_insert(
     385             :               temporal_id, std::move(remote_boundary_data));
     386             :           continue;
     387             :         }
     388             : 
     389             :         // Copy the central element's mortar data to the neighbor
     390             :         auto oriented_mortar_data =
     391             :             central_mortar_data_.at(mortar_id).local_data(temporal_id);
     392             :         if (not orientation.is_aligned()) {
     393             :           oriented_mortar_data.orient_on_slice(
     394             :               mortar_mesh.extents(), direction.dimension(), orientation);
     395             :         }
     396             :         neighbors_mortar_data_[overlap_id][::dg::MortarId<Dim>{
     397             :                                                direction_from_neighbor,
     398             :                                                central_element.id()}]
     399             :             .remote_insert(temporal_id, std::move(oriented_mortar_data));
     400             : 
     401             :         // Now we switch perspective to the neighbor. First, we extend the
     402             :         // overlap data to the full neighbor mesh by padding it with zeros. This
     403             :         // is necessary because spectral operators such as derivatives require
     404             :         // data on the full mesh.
     405             :         LinearSolver::Schwarz::extended_overlap_data(
     406             :             make_not_null(&extended_operand_vars_[overlap_id]),
     407             :             operand.overlap_data.at(overlap_id), neighbor_mesh.extents(),
     408             :             overlap_extent, direction_from_neighbor);
     409             : 
     410             :         const auto apply_boundary_condition_neighbor =
     411             :             [&apply_boundary_condition, &local_neighbor_id = neighbor_id,
     412             :              &overlap_id](const Direction<Dim>& local_direction,
     413             :                           const auto... fields_and_fluxes) {
     414             :               apply_boundary_condition(
     415             :                   local_neighbor_id, local_direction, std::true_type{},
     416             :                   std::forward_as_tuple(overlap_id, local_direction),
     417             :                   fields_and_fluxes...);
     418             :             };
     419             : 
     420             :         const auto fluxes_args_on_overlap =
     421             :             elliptic::util::apply_at<fluxes_args_tags_overlap,
     422             :                                      args_tags_from_center>(get_items, box,
     423             :                                                             overlap_id);
     424             :         const auto sources_args_on_overlap =
     425             :             elliptic::util::apply_at<sources_args_tags_overlap,
     426             :                                      args_tags_from_center>(get_items, box,
     427             :                                                             overlap_id);
     428             :         DirectionMap<Dim, FluxesArgs> fluxes_args_on_overlap_faces{};
     429             :         for (const auto& neighbor_direction :
     430             :              Direction<Dim>::all_directions()) {
     431             :           fluxes_args_on_overlap_faces.emplace(
     432             :               neighbor_direction,
     433             :               elliptic::util::apply_at<fluxes_args_tags_overlap_faces,
     434             :                                        args_tags_from_center>(
     435             :                   get_items, box,
     436             :                   std::forward_as_tuple(overlap_id, neighbor_direction)));
     437             :         }
     438             : 
     439             :         elliptic::util::apply_at<prepare_args_tags_overlap,
     440             :                                  args_tags_from_center>(
     441             :             [this, &overlap_id](const auto&... args) {
     442             :               elliptic::dg::prepare_mortar_data<System, linearized>(
     443             :                   make_not_null(&neighbors_auxiliary_vars_[overlap_id]),
     444             :                   make_not_null(&neighbors_auxiliary_fluxes_[overlap_id]),
     445             :                   make_not_null(&neighbors_primal_fluxes_[overlap_id]),
     446             :                   make_not_null(&neighbors_mortar_data_[overlap_id]),
     447             :                   extended_operand_vars_[overlap_id], args...);
     448             :             },
     449             :             box, overlap_id, temporal_id, apply_boundary_condition_neighbor,
     450             :             fluxes_args_on_overlap, sources_args_on_overlap,
     451             :             fluxes_args_on_overlap_faces);
     452             : 
     453             :         // Copy this neighbor's mortar data to the other side of the mortars. On
     454             :         // the other side we either have the central element, or another element
     455             :         // that may or may not be part of the subdomain.
     456             :         const auto& neighbor_mortar_meshes =
     457             :             all_neighbor_mortar_meshes.at(overlap_id);
     458             :         for (const auto& neighbor_mortar_id_and_data :
     459             :              neighbors_mortar_data_.at(overlap_id)) {
     460             :           // No structured bindings because capturing these in lambdas doesn't
     461             :           // work until C++20
     462             :           const auto& neighbor_mortar_id = neighbor_mortar_id_and_data.first;
     463             :           const auto& neighbor_mortar_data = neighbor_mortar_id_and_data.second;
     464             :           const auto& neighbor_direction = neighbor_mortar_id.first;
     465             :           const auto& neighbors_neighbor_id = neighbor_mortar_id.second;
     466             :           // No need to do anything on external boundaries
     467             :           if (neighbors_neighbor_id == ElementId<Dim>::external_boundary_id()) {
     468             :             continue;
     469             :           }
     470             :           const auto& neighbor_orientation =
     471             :               neighbor.neighbors().at(neighbor_direction).orientation();
     472             :           const auto neighbors_neighbor_direction =
     473             :               neighbor_orientation(neighbor_direction.opposite());
     474             :           const ::dg::MortarId<Dim> mortar_id_from_neighbors_neighbor{
     475             :               neighbors_neighbor_direction, neighbor_id};
     476             :           const auto send_mortar_data =
     477             :               [&neighbor_orientation, &neighbor_mortar_meshes,
     478             :                &neighbor_mortar_data, &neighbor_mortar_id,
     479             :                &neighbor_direction](auto& remote_mortar_data) {
     480             :                 const auto& neighbor_mortar_mesh =
     481             :                     neighbor_mortar_meshes.at(neighbor_mortar_id);
     482             :                 auto oriented_neighbor_mortar_data =
     483             :                     neighbor_mortar_data.local_data(temporal_id);
     484             :                 if (not neighbor_orientation.is_aligned()) {
     485             :                   oriented_neighbor_mortar_data.orient_on_slice(
     486             :                       neighbor_mortar_mesh.extents(),
     487             :                       neighbor_direction.dimension(), neighbor_orientation);
     488             :                 }
     489             :                 remote_mortar_data.remote_insert(
     490             :                     temporal_id, std::move(oriented_neighbor_mortar_data));
     491             :               };
     492             :           if (neighbors_neighbor_id == central_element.id() and
     493             :               mortar_id_from_neighbors_neighbor == mortar_id) {
     494             :             send_mortar_data(central_mortar_data_.at(mortar_id));
     495             :             continue;
     496             :           }
     497             :           // Determine whether the neighbor's neighbor overlaps with the
     498             :           // subdomain and find its overlap ID if it does.
     499             :           const auto neighbors_neighbor_overlap_id =
     500             :               [&local_all_neighbor_mortar_meshes = all_neighbor_mortar_meshes,
     501             :                &neighbors_neighbor_id, &mortar_id_from_neighbors_neighbor]()
     502             :               -> std::optional<LinearSolver::Schwarz::OverlapId<Dim>> {
     503             :             for (const auto& [local_overlap_id, local_mortar_meshes] :
     504             :                  local_all_neighbor_mortar_meshes) {
     505             :               if (local_overlap_id.second != neighbors_neighbor_id) {
     506             :                 continue;
     507             :               }
     508             :               for (const auto& local_mortar_id_and_mesh : local_mortar_meshes) {
     509             :                 if (local_mortar_id_and_mesh.first ==
     510             :                     mortar_id_from_neighbors_neighbor) {
     511             :                   return local_overlap_id;
     512             :                 }
     513             :               }
     514             :             }
     515             :             return std::nullopt;
     516             :           }();
     517             :           if (neighbors_neighbor_overlap_id.has_value()) {
     518             :             // The neighbor's neighbor is part of the subdomain so we copy the
     519             :             // mortar data over. Once the loop is complete we will also have
     520             :             // received mortar data back. At that point, both neighbors have a
     521             :             // copy of each other's mortar data, which is the subject of the
     522             :             // possible optimizations mentioned above. Note that the data may
     523             :             // differ by orientations.
     524             :             send_mortar_data(
     525             :                 neighbors_mortar_data_[*neighbors_neighbor_overlap_id]
     526             :                                       [mortar_id_from_neighbors_neighbor]);
     527             :           } else {
     528             :             // The neighbor's neighbor does not overlap with the subdomain, so
     529             :             // we don't copy mortar data and also don't expect to receive any.
     530             :             // Instead, we assume the data on it is zero and manufacture
     531             :             // appropriate remote boundary data.
     532             :             const auto& neighbors_neighbor_mortar_mesh =
     533             :                 all_neighbors_neighbor_mortar_meshes.at(overlap_id)
     534             :                     .at(neighbor_mortar_id);
     535             :             auto zero_mortar_data = elliptic::dg::zero_boundary_data_on_mortar<
     536             :                 typename System::primal_fields, typename System::primal_fluxes>(
     537             :                 neighbors_neighbor_direction,
     538             :                 all_neighbors_neighbor_meshes.at(overlap_id)
     539             :                     .at(neighbor_mortar_id),
     540             :                 all_neighbors_neighbor_face_normal_magnitudes.at(overlap_id)
     541             :                     .at(neighbor_mortar_id),
     542             :                 neighbors_neighbor_mortar_mesh,
     543             :                 all_neighbors_neighbor_mortar_sizes.at(overlap_id)
     544             :                     .at(neighbor_mortar_id));
     545             :             // The data is zero, but auxiliary quantities such as the face
     546             :             // normal magnitude may need re-orientation
     547             :             if (not neighbor_orientation.is_aligned()) {
     548             :               zero_mortar_data.orient_on_slice(
     549             :                   neighbors_neighbor_mortar_mesh.extents(),
     550             :                   neighbors_neighbor_direction.dimension(),
     551             :                   neighbor_orientation.inverse_map());
     552             :             }
     553             :             neighbors_mortar_data_.at(overlap_id)
     554             :                 .at(neighbor_mortar_id)
     555             :                 .remote_insert(temporal_id, std::move(zero_mortar_data));
     556             :           }
     557             :         }  // loop over neighbor's mortars
     558             :       }    // loop over neighbors
     559             :     }      // loop over directions
     560             : 
     561             :     // 2. Apply the operator on all elements in the subdomain
     562             :     //
     563             :     // Apply on central element
     564             :     db::apply<apply_args_tags>(
     565             :         [this, &result, &operand](const auto&... args) {
     566             :           elliptic::dg::apply_operator<System, linearized>(
     567             :               make_not_null(&result->element_data),
     568             :               make_not_null(&central_mortar_data_), operand.element_data,
     569             :               central_primal_fluxes_, args...);
     570             :         },
     571             :         box, temporal_id, sources_args);
     572             :     // Apply on neighbors
     573             :     for (const auto& [direction, neighbors] : central_element.neighbors()) {
     574             :       const auto& orientation = neighbors.orientation();
     575             :       const auto direction_from_neighbor = orientation(direction.opposite());
     576             :       for (const auto& neighbor_id : neighbors) {
     577             :         const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
     578             :                                                                neighbor_id};
     579             :         const auto& overlap_extent = all_overlap_extents.at(overlap_id);
     580             :         const auto& neighbor_mesh = all_neighbor_meshes.at(overlap_id);
     581             : 
     582             :         if (UNLIKELY(overlap_extent == 0)) {
     583             :           continue;
     584             :         }
     585             : 
     586             :         elliptic::util::apply_at<apply_args_tags_overlap,
     587             :                                  args_tags_from_center>(
     588             :             [this, &overlap_id](const auto&... args) {
     589             :               elliptic::dg::apply_operator<System, linearized>(
     590             :                   make_not_null(&extended_results_[overlap_id]),
     591             :                   make_not_null(&neighbors_mortar_data_.at(overlap_id)),
     592             :                   extended_operand_vars_.at(overlap_id),
     593             :                   neighbors_primal_fluxes_.at(overlap_id), args...);
     594             :             },
     595             :             box, overlap_id, temporal_id,
     596             :             elliptic::util::apply_at<sources_args_tags_overlap,
     597             :                                      args_tags_from_center>(get_items, box,
     598             :                                                             overlap_id));
     599             : 
     600             :         // Restrict the extended operator data back to the subdomain, assuming
     601             :         // we can discard any data outside the overlaps. WARNING: This
     602             :         // assumption may break with changes to the DG operator that affect its
     603             :         // sparsity. For example, multiplying the DG operator with the _full_
     604             :         // inverse mass-matrix ("massless" scheme with no "mass-lumping"
     605             :         // approximation) means that lifted boundary corrections bleed into the
     606             :         // volume.
     607             :         if (UNLIKELY(
     608             :                 result->overlap_data[overlap_id].number_of_grid_points() !=
     609             :                 operand.overlap_data.at(overlap_id).number_of_grid_points())) {
     610             :           result->overlap_data[overlap_id].initialize(
     611             :               operand.overlap_data.at(overlap_id).number_of_grid_points());
     612             :         }
     613             :         LinearSolver::Schwarz::data_on_overlap(
     614             :             make_not_null(&result->overlap_data[overlap_id]),
     615             :             extended_results_.at(overlap_id), neighbor_mesh.extents(),
     616             :             overlap_extent, direction_from_neighbor);
     617             :       }  // loop over neighbors
     618             :     }    // loop over directions
     619             :   }
     620             : 
     621             :   // NOLINTNEXTLINE(google-runtime-references)
     622           0 :   void pup(PUP::er& p) { p | override_boundary_conditions_; }
     623             : 
     624             :  private:
     625             :   std::optional<std::unique_ptr<BoundaryConditionsBase>>
     626           0 :       override_boundary_conditions_{};
     627             : 
     628             :   // Memory buffers for repeated operator applications
     629             :   mutable Variables<typename System::auxiliary_fields>
     630           0 :       central_auxiliary_vars_{};
     631           0 :   mutable Variables<typename System::primal_fluxes> central_primal_fluxes_{};
     632             :   mutable Variables<typename System::auxiliary_fluxes>
     633           0 :       central_auxiliary_fluxes_{};
     634             :   mutable LinearSolver::Schwarz::OverlapMap<
     635             :       Dim, Variables<typename System::auxiliary_fields>>
     636           0 :       neighbors_auxiliary_vars_{};
     637             :   mutable LinearSolver::Schwarz::OverlapMap<
     638             :       Dim, Variables<typename System::primal_fluxes>>
     639           0 :       neighbors_primal_fluxes_{};
     640             :   mutable LinearSolver::Schwarz::OverlapMap<
     641             :       Dim, Variables<typename System::auxiliary_fluxes>>
     642           0 :       neighbors_auxiliary_fluxes_{};
     643             :   mutable LinearSolver::Schwarz::OverlapMap<
     644             :       Dim, Variables<typename System::primal_fields>>
     645           0 :       extended_operand_vars_{};
     646             :   mutable LinearSolver::Schwarz::OverlapMap<
     647             :       Dim, Variables<typename System::primal_fields>>
     648           0 :       extended_results_{};
     649             :   mutable ::dg::MortarMap<
     650             :       Dim, elliptic::dg::MortarData<size_t, typename System::primal_fields,
     651             :                                     typename System::primal_fluxes>>
     652           0 :       central_mortar_data_{};
     653             :   mutable LinearSolver::Schwarz::OverlapMap<
     654             :       Dim, ::dg::MortarMap<Dim, elliptic::dg::MortarData<
     655             :                                     size_t, typename System::primal_fields,
     656             :                                     typename System::primal_fluxes>>>
     657           0 :       neighbors_mortar_data_{};
     658             : };
     659             : 
     660             : }  // namespace elliptic::dg::subdomain_operator

Generated by: LCOV version 1.14