Communicate volume mesh in boundary data

src/Evolution/DgSubcell/Actions/Initialize.hpp

@@ -30,6 +30,7 @@
 #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
 #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
 #include "Evolution/DgSubcell/Tags/ReconstructionOrder.hpp"
 #include "Evolution/DgSubcell/Tags/StepsSinceTciCall.hpp"
 #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
@@ -68,6 +69,7 @@ namespace evolution::dg::subcell::Actions {
  *   - `System::variables_tag`
  * - Adds:
  *   - `subcell::Tags::Mesh<Dim>`
+ *   - `subcell::Tags::MeshForGhostData<Dim>`
  *   - `subcell::Tags::ActiveGrid`
  *   - `subcell::Tags::DidRollback`
  *   - `subcell::Tags::TciGridHistory`
@@ -94,6 +96,7 @@ struct SetSubcellGrid {
   using simple_tags = tmpl::list<
       Tags::ActiveGrid, Tags::DidRollback, Tags::TciGridHistory,
       Tags::TciCallsSinceRollback, Tags::StepsSinceTciCall,
+      evolution::dg::subcell::Tags::MeshForGhostData<Dim>,
       Tags::GhostDataForReconstruction<Dim>, Tags::TciDecision,
       Tags::NeighborTciDecisions<Dim>, Tags::DataForRdmpTci,
       subcell::Tags::CellCenteredFlux<typename System::flux_variables, Dim>,

src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp

@@ -43,6 +43,7 @@
 #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
 #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
 #include "Evolution/DgSubcell/Tags/Reconstructor.hpp"
 #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
 #include "Evolution/DiscontinuousGalerkin/BoundaryData.hpp"
@@ -257,6 +258,7 @@ struct SendDataForReconstruction {
                                 rdmp_tci_data.min_variables_values.size())));
 
         evolution::dg::BoundaryData<Dim> data{
+            subcell_mesh,
             subcell_mesh,
             dg_mesh.slice_away(direction.dimension()),
             std::move(subcell_data_to_send),
@@ -356,6 +358,7 @@ struct ReceiveDataForReconstruction {
                evolution::dg::Tags::MortarData<Dim>,
                evolution::dg::Tags::MortarNextTemporalId<Dim>,
                domain::Tags::NeighborMesh<Dim>,
+               evolution::dg::subcell::Tags::MeshForGhostData<Dim>,
                evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
         [&element,
          ghost_zone_size =
@@ -372,11 +375,14 @@ struct ReceiveDataForReconstruction {
                 mortar_next_time_step_id,
             const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
                 neighbor_mesh,
+            const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
+                mesh_for_ghost_data,
             const auto neighbor_tci_decisions,
             const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
                 neighbor_dg_to_fd_interpolants) {
           // Remove neighbor meshes for neighbors that don't exist anymore
           domain::remove_nonexistent_neighbors(neighbor_mesh, element);
+          domain::remove_nonexistent_neighbors(mesh_for_ghost_data, element);
 
           // Get the next time step id, and also the fluxes data if the neighbor
           // is doing DG.
@@ -402,8 +408,10 @@ struct ReceiveDataForReconstruction {
             }
             // Set new neighbor mesh
             neighbor_mesh->insert_or_assign(
-                mortar_id,
-                received_mortar_data.second.volume_mesh_ghost_cell_data);
+                mortar_id, received_mortar_data.second.volume_mesh);
+            mesh_for_ghost_data->insert_or_assign(
+                mortar_id, received_mortar_data.second
+                               .volume_mesh_ghost_cell_data.value());
           }
 
           ASSERT(ghost_data_ptr->empty(),
@@ -437,7 +445,7 @@ struct ReceiveDataForReconstruction {
                   rdmp_tci_data_ptr, ghost_data_ptr,
                   *received_data[directional_element_id].ghost_cell_data,
                   number_of_rdmp_vars, directional_element_id,
-                  neighbor_mesh->at(directional_element_id), element,
+                  mesh_for_ghost_data->at(directional_element_id), element,
                   subcell_mesh, ghost_zone_size,
                   neighbor_dg_to_fd_interpolants);
               ASSERT(neighbor_tci_decisions->contains(directional_element_id),

src/Evolution/DgSubcell/Actions/TciAndRollback.hpp

@@ -22,7 +22,6 @@
 #include "Domain/Structure/Element.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Tags.hpp"
-#include "Domain/Tags/NeighborMesh.hpp"
 #include "Evolution/DgSubcell/Actions/Labels.hpp"
 #include "Evolution/DgSubcell/ActiveGrid.hpp"
 #include "Evolution/DgSubcell/GhostData.hpp"
@@ -38,6 +37,7 @@
 #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
 #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
 #include "Evolution/DgSubcell/Tags/Reconstructor.hpp"
 #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
 #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
@@ -261,7 +261,8 @@ struct TciAndRollback {
           // method, since we need to lift G+D instead of the ingredients
           // that go into G+D, which is what we would be projecting here.
         },
-        make_not_null(&box), db::get<domain::Tags::NeighborMesh<Dim>>(box),
+        make_not_null(&box),
+        db::get<evolution::dg::subcell::Tags::MeshForGhostData<Dim>>(box),
         db::get<evolution::dg::subcell::Tags::Reconstructor>(box)
             .ghost_zone_size(),
         db::get<

src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp

@@ -21,6 +21,7 @@
 #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
 #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
 #include "Evolution/DiscontinuousGalerkin/BoundaryData.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "Time/TimeStepId.hpp"
@@ -35,9 +36,11 @@ void neighbor_reconstructed_face_solution(
         TimeStepId,
         DirectionalIdMap<VolumeDim, evolution::dg::BoundaryData<VolumeDim>>>*>
         received_temporal_id_and_data) {
-  db::mutate<subcell::Tags::GhostDataForReconstruction<VolumeDim>,
+  db::mutate<subcell::Tags::MeshForGhostData<VolumeDim>,
+             subcell::Tags::GhostDataForReconstruction<VolumeDim>,
              subcell::Tags::DataForRdmpTci>(
-      [&received_temporal_id_and_data](const auto subcell_ghost_data_ptr,
+      [&received_temporal_id_and_data](const auto mesh_for_ghost_data_ptr,
+                                       const auto subcell_ghost_data_ptr,
                                        const auto rdmp_tci_data_ptr) {
         const size_t number_of_evolved_vars =
             rdmp_tci_data_ptr->max_variables_values.size();
@@ -48,6 +51,9 @@ void neighbor_reconstructed_face_solution(
                  "The subcell mortar data was not sent at TimeStepId "
                      << received_temporal_id_and_data->first
                      << " with mortar id " << mortar_id);
+          mesh_for_ghost_data_ptr->insert_or_assign(
+              mortar_id,
+              received_mortar_data.second.volume_mesh_ghost_cell_data.value());
           const DataVector& neighbor_ghost_and_subcell_data =
               received_mortar_data.second.ghost_cell_data.value();
           // Compute min and max over neighbors

src/Evolution/DgSubcell/PrepareNeighborData.hpp

@@ -5,6 +5,7 @@
 
 #include <array>
 #include <cstddef>
+#include <optional>
 #include <utility>
 
 #include "DataStructures/DataBox/DataBox.hpp"
@@ -62,7 +63,7 @@ template <typename Metavariables, typename DbTagsList, size_t Dim>
 void prepare_neighbor_data(
     const gsl::not_null<DirectionMap<Dim, DataVector>*>
         all_neighbor_data_for_reconstruction,
-    const gsl::not_null<Mesh<Dim>*> ghost_data_mesh,
+    const gsl::not_null<std::optional<Mesh<Dim>>*> ghost_data_mesh,
     const gsl::not_null<db::DataBox<DbTagsList>*> box,
     [[maybe_unused]] const Variables<db::wrap_tags_in<
         ::Tags::Flux, typename Metavariables::system::flux_variables,

src/Evolution/DgSubcell/Tags/CMakeLists.txt

@@ -18,6 +18,7 @@ spectre_target_headers(
   Interpolators.hpp
   Jacobians.hpp
   Mesh.hpp
+  MeshForGhostData.hpp
   MethodOrder.hpp
   ObserverCoordinates.hpp
   ObserverMesh.hpp

src/Evolution/DgSubcell/Tags/MeshForGhostData.hpp

@@ -0,0 +1,29 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/DataBox/Tag.hpp"
+
+/// \cond
+template <size_t Dim, typename T>
+class DirectionalIdMap;
+template <size_t Dim>
+class Mesh;
+/// \endcond
+
+namespace evolution::dg::subcell::Tags {
+/*!
+ * \brief Holds the volume Mesh used by each neighbor for communicating subcell
+ * ghost data.
+ *
+ * \details This will be the subcell Mesh unless an Element doing DG has only
+ * neighbors also doing DG
+ */
+template <size_t Dim>
+struct MeshForGhostData : db::SimpleTag {
+  using type = DirectionalIdMap<Dim, ::Mesh<Dim>>;
+};
+}  // namespace evolution::dg::subcell::Tags

src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp

@@ -123,14 +123,14 @@ void retrieve_boundary_data_spsc(
         auto& current_inbox = (*boundary_data_ptr)[time_step_id];
         if (auto it = current_inbox.find(directional_element_id);
             it != current_inbox.end()) {
-          auto& [volume_mesh_of_ghost_cell_data, face_mesh, ghost_cell_data,
-                 boundary_data, boundary_data_validity_range,
+          auto& [volume_mesh, volume_mesh_of_ghost_cell_data, face_mesh,
+                 ghost_cell_data, boundary_data, boundary_data_validity_range,
                  boundary_tci_status, boundary_integration_order] = data;
           (void)ghost_cell_data;
-          auto& [current_volume_mesh_of_ghost_cell_data, current_face_mesh,
-                 current_ghost_cell_data, current_boundary_data,
-                 current_boundary_data_validity_range, current_tci_status,
-                 current_integration_order] = it->second;
+          auto& [current_volume_mesh, current_volume_mesh_of_ghost_cell_data,
+                 current_face_mesh, current_ghost_cell_data,
+                 current_boundary_data, current_boundary_data_validity_range,
+                 current_tci_status, current_integration_order] = it->second;
           // Need to use when optimizing subcell
           (void)current_volume_mesh_of_ghost_cell_data;
           // We have already received some data at this time. Receiving
@@ -326,8 +326,7 @@ bool receive_boundary_data_global_time_stepping(
              received_temporal_id_and_data.second) {
           const auto& mortar_id = received_mortar_data.first;
           neighbor_mesh->insert_or_assign(
-              mortar_id,
-              received_mortar_data.second.volume_mesh_ghost_cell_data);
+              mortar_id, received_mortar_data.second.volume_mesh);
           mortar_next_time_step_id->at(mortar_id) =
               received_mortar_data.second.validity_range;
           ASSERT(using_subcell_v<Metavariables> or
@@ -466,8 +465,7 @@ bool receive_boundary_data_local_time_stepping(
                        << mortar_id
                        << "\nTimeStepId: " << mortar_next_time_step_id);
             neighbor_mesh->insert_or_assign(
-                mortar_id,
-                received_mortar_data->second.volume_mesh_ghost_cell_data);
+                mortar_id, received_mortar_data->second.volume_mesh);
             neighbor_mortar_data.face_mesh =
                 received_mortar_data->second.interface_mesh;
             neighbor_mortar_data.mortar_mesh = mortar_meshes.at(mortar_id);

src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp

@@ -675,9 +675,8 @@ void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
   std::optional<DirectionMap<Dim, DataVector>>
       all_neighbor_data_for_reconstruction = std::nullopt;
   int tci_decision = 0;
-  // Set ghost_cell_mesh to the DG mesh, then update it below if we did a
-  // projection.
-  Mesh<Dim> ghost_data_mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
+  const Mesh<Dim>& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
+  std::optional<Mesh<Dim>> ghost_data_mesh = std::nullopt;
   if constexpr (using_subcell_v<Metavariables>) {
     if (not all_neighbor_data_for_reconstruction.has_value()) {
       all_neighbor_data_for_reconstruction = DirectionMap<Dim, DataVector>{};
@@ -734,15 +733,17 @@ void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
       SendData data{};
 
       if (neighbor_count == total_neighbors) {
-        data = SendData{ghost_data_mesh,
+        data = SendData{volume_mesh,
+                        ghost_data_mesh,
                         face_mesh_for_neighbor,
                         std::move(ghost_and_subcell_data),
                         {std::move(neighbor_boundary_data_on_mortar)},
                         next_time_step_id,
                         tci_decision,
                         integration_order};
       } else {
-        data = SendData{ghost_data_mesh,
+        data = SendData{volume_mesh,
+                        ghost_data_mesh,
                         face_mesh_for_neighbor,
                         ghost_and_subcell_data,
                         {std::move(neighbor_boundary_data_on_mortar)},
@@ -804,7 +805,7 @@ void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
                evolution::dg::Tags::MortarDataHistory<
                    Dim, typename dt_variables_tag::type>>(
         [&element, integration_order, &time_step_id, using_gauss_points,
-         &volume_det_inv_jacobian](
+         &volume_det_inv_jacobian, &volume_mesh](
             const gsl::not_null<
                 DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
                 mortar_data,
@@ -814,7 +815,6 @@ void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
                                           evolution::dg::MortarData<Dim>,
                                           typename dt_variables_tag::type>>*>
                 boundary_data_history,
-            const Mesh<Dim>& volume_mesh,
             const DirectionMap<Dim,
                                std::optional<Variables<tmpl::list<
                                    evolution::dg::Tags::MagnitudeOfNormal,
@@ -878,7 +878,7 @@ void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
             }
           }
         },
-        box, db::get<domain::Tags::Mesh<Dim>>(*box),
+        box,
         db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box));
   }
 }

src/Evolution/DiscontinuousGalerkin/BoundaryData.cpp

@@ -13,6 +13,7 @@
 namespace evolution::dg {
 template <size_t Dim>
 void BoundaryData<Dim>::pup(PUP::er& p) {
+  p | volume_mesh;
   p | volume_mesh_ghost_cell_data;
   p | interface_mesh;
   p | ghost_cell_data;
@@ -24,7 +25,8 @@ void BoundaryData<Dim>::pup(PUP::er& p) {
 
 template <size_t Dim>
 bool operator==(const BoundaryData<Dim>& lhs, const BoundaryData<Dim>& rhs) {
-  return lhs.volume_mesh_ghost_cell_data == rhs.volume_mesh_ghost_cell_data and
+  return lhs.volume_mesh == rhs.volume_mesh and
+         lhs.volume_mesh_ghost_cell_data == rhs.volume_mesh_ghost_cell_data and
          lhs.interface_mesh == rhs.interface_mesh and
          lhs.ghost_cell_data == rhs.ghost_cell_data and
          lhs.boundary_correction_data == rhs.boundary_correction_data and
@@ -40,7 +42,8 @@ bool operator!=(const BoundaryData<Dim>& lhs, const BoundaryData<Dim>& rhs) {
 
 template <size_t Dim>
 std::ostream& operator<<(std::ostream& os, const BoundaryData<Dim>& value) {
-  return os << "Ghost mesh: " << value.volume_mesh_ghost_cell_data << '\n'
+  return os << "Volume mesh: " << value.volume_mesh << '\n'
+            << "Ghost mesh: " << value.volume_mesh_ghost_cell_data << '\n'
             << "Interface mesh: " << value.interface_mesh << '\n'
             << "Ghost cell data: " << value.ghost_cell_data << '\n'
             << "Boundary correction: " << value.boundary_correction_data << '\n'

src/Evolution/DiscontinuousGalerkin/BoundaryData.hpp

@@ -17,25 +17,27 @@ namespace evolution::dg {
  *
  * The stored data consists of the following:
  *
- * 1. the mesh of the ghost cell data we received. This allows eliding
+ * 1. the volume mesh of the element.
+ * 2. the volume mesh corresponding to the ghost cell data. This allows eliding
  *    projection when all neighboring elements are doing DG.
- * 2. the mesh of the neighboring element's face (not the mortar mesh!)
- * 3. the variables at the ghost zone cells for finite difference/volume
+ * 3. the mesh of the neighboring element's face (not the mortar mesh!)
+ * 4. the variables at the ghost zone cells for finite difference/volume
  *    reconstruction
- * 4. the data on the mortar needed for computing the boundary corrections (e.g.
+ * 5. the data on the mortar needed for computing the boundary corrections (e.g.
  *    fluxes, characteristic speeds, conserved variables)
- * 5. the TimeStepId beyond which the boundary terms are no longer valid, when
+ * 6. the TimeStepId beyond which the boundary terms are no longer valid, when
  *    using local time stepping.
- * 6. the troubled cell indicator status used for determining halos around
+ * 7. the troubled cell indicator status used for determining halos around
  *    troubled cells.
- * 7. the integration order of the time-stepper
+ * 8. the integration order of the time-stepper
  */
 template <size_t Dim>
 struct BoundaryData {
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p);
 
-  Mesh<Dim> volume_mesh_ghost_cell_data{};
+  Mesh<Dim> volume_mesh{};
+  std::optional<Mesh<Dim>> volume_mesh_ghost_cell_data{};
   Mesh<Dim - 1> interface_mesh{};
   std::optional<DataVector> ghost_cell_data{};
   std::optional<DataVector> boundary_correction_data{};