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

Generated by: LCOV version 1.14