Merge pull request #6230 from kidder/project_mortar_data

Project mortar data

.clang-tidy

@@ -13,6 +13,10 @@ Checks: '*,
 # complains about Catch macros having infinite loops,
 -bugprone-infinite-loop,
 -bugprone-macro-parentheses,
+# useless as it triggers on calls to value() which already do a check,
+# and is too aggressive as it requires the check to be immediately before the
+# access.
+-bugprone-unchecked-optional-access,
 # Official justification is not public, unofficial justification has nothing,
 # specific to iterators and is basically "never create mutable temporaries",,
 # which is against standard practice everywhere else.,

src/Domain/Structure/OrientationMapHelpers.cpp

@@ -315,6 +315,30 @@ void orient_each_component(
 }
 }  // namespace
 
+template <size_t VolumeDim>
+Mesh<VolumeDim - 1> orient_mesh_on_slice(
+    const Mesh<VolumeDim - 1>& mesh_on_slice, const size_t sliced_dim,
+    const OrientationMap<VolumeDim>& orientation_of_neighbor) {
+  if constexpr (VolumeDim < 3) {
+    return mesh_on_slice;
+  } else {
+    if (orientation_of_neighbor.is_aligned()) {
+      return mesh_on_slice;
+    }
+    const size_t first_dim_of_slice = sliced_dim == 0 ? 1 : 0;
+    const size_t second_dim_of_slice = sliced_dim == 2 ? 1 : 2;
+    if (orientation_of_neighbor(first_dim_of_slice) >
+        orientation_of_neighbor(second_dim_of_slice)) {
+      return Mesh<2>{
+          {{mesh_on_slice.extents(1), mesh_on_slice.extents(0)}},
+          {{mesh_on_slice.basis(1), mesh_on_slice.basis(0)}},
+          {{mesh_on_slice.quadrature(1), mesh_on_slice.quadrature(0)}}};
+    } else {
+      return mesh_on_slice;
+    }
+  }
+}
+
 template <size_t VolumeDim>
 void orient_variables(
     const gsl::not_null<DataVector*> result, const DataVector& variables,
@@ -428,6 +452,9 @@ DataVector orient_variables_on_slice(
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
 #define INSTANTIATION(r, data)                                               \
+  template Mesh<DIM(data) - 1> orient_mesh_on_slice(                         \
+      const Mesh<DIM(data) - 1>& mesh_on_slice, const size_t sliced_dim,     \
+      const OrientationMap<DIM(data)>& orientation_of_neighbor);             \
   template void orient_variables(                                            \
       const gsl::not_null<DataVector*> result, const DataVector& variables,  \
       const Index<DIM(data)>& extents,                                       \

src/Domain/Structure/OrientationMapHelpers.hpp

@@ -13,13 +13,23 @@
 
 /// \cond
 template <size_t VolumeDim>
+class Mesh;
+template <size_t VolumeDim>
 class OrientationMap;
 template <size_t>
 class Index;
 template <typename TagsList>
 class Variables;
 /// \endcond
 
+/// \ingroup ComputationalDomainGroup
+/// \brief Orient a sliced Mesh to the logical frame of  a neighbor element with
+/// the given orientation.
+template <size_t VolumeDim>
+Mesh<VolumeDim - 1> orient_mesh_on_slice(
+    const Mesh<VolumeDim - 1>& mesh_on_slice, size_t sliced_dim,
+    const OrientationMap<VolumeDim>& orientation_of_neighbor);
+
 /// @{
 /// \ingroup ComputationalDomainGroup
 /// \brief Orient variables to the data-storage order of a neighbor element with

src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp

@@ -337,6 +337,7 @@ struct ReceiveDataForReconstruction {
     inbox.erase(current_time_step_id);
 
     const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
+    const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(box);
 
     db::mutate<Tags::GhostDataForReconstruction<Dim>, Tags::DataForRdmpTci,
                evolution::dg::Tags::MortarData<Dim>,
@@ -347,7 +348,7 @@ struct ReceiveDataForReconstruction {
          ghost_zone_size =
              db::get<evolution::dg::subcell::Tags::Reconstructor>(box)
                  .ghost_zone_size(),
-         &received_data, &subcell_mesh](
+         &received_data, &subcell_mesh, &mortar_meshes](
             const gsl::not_null<DirectionalIdMap<Dim, GhostData>*>
                 ghost_data_ptr,
             const gsl::not_null<RdmpTciData*> rdmp_tci_data_ptr,
@@ -381,6 +382,8 @@ struct ReceiveDataForReconstruction {
                     .has_value()) {
               mortar_data->at(mortar_id).neighbor().face_mesh =
                   received_mortar_data.second.interface_mesh;
+              mortar_data->at(mortar_id).neighbor().mortar_mesh =
+                  mortar_meshes.at(mortar_id);
               mortar_data->at(mortar_id).neighbor().mortar_data = std::move(
                   *received_mortar_data.second.boundary_correction_data);
             }

src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp

@@ -306,10 +306,12 @@ bool receive_boundary_data_global_time_stepping(
         make_not_null(&db::as_access(*box)), received_temporal_id_and_data);
   }
 
+  const auto& mortar_meshes =
+      get<evolution::dg::Tags::MortarMesh<volume_dim>>(*box);
   db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
              evolution::dg::Tags::MortarNextTemporalId<volume_dim>,
              domain::Tags::NeighborMesh<volume_dim>>(
-      [&received_temporal_id_and_data](
+      [&received_temporal_id_and_data, &mortar_meshes](
           const gsl::not_null<DirectionalIdMap<
               volume_dim, evolution::dg::MortarDataHolder<volume_dim>>*>
               mortar_data,
@@ -338,6 +340,8 @@ bool receive_boundary_data_global_time_stepping(
                   .has_value()) {
             mortar_data->at(mortar_id).neighbor().face_mesh =
                 received_mortar_data.second.interface_mesh;
+            mortar_data->at(mortar_id).neighbor().mortar_mesh =
+                mortar_meshes.at(mortar_id);
             mortar_data->at(mortar_id).neighbor().mortar_data = std::move(
                 received_mortar_data.second.boundary_correction_data.value());
           }
@@ -411,13 +415,14 @@ bool receive_boundary_data_local_time_stepping(
   }
   ASSERT(inbox_ptr != nullptr, "The inbox pointer should not be null.");
   InboxMap& inbox = *inbox_ptr;
+  const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(*box);
 
   const bool have_all_intermediate_messages = db::mutate<
       evolution::dg::Tags::MortarDataHistory<Dim,
                                              typename dt_variables_tag::type>,
       evolution::dg::Tags::MortarNextTemporalId<Dim>,
       domain::Tags::NeighborMesh<Dim>>(
-      [&inbox, &needed_time](
+      [&inbox, &needed_time, &mortar_meshes](
           const gsl::not_null<DirectionalIdMap<
               Dim,
               TimeSteppers::BoundaryHistory<evolution::dg::MortarData<Dim>,
@@ -464,6 +469,7 @@ bool receive_boundary_data_local_time_stepping(
                 received_mortar_data->second.volume_mesh_ghost_cell_data);
             neighbor_mortar_data.face_mesh =
                 received_mortar_data->second.interface_mesh;
+            neighbor_mortar_data.mortar_mesh = mortar_meshes.at(mortar_id);
             neighbor_mortar_data.mortar_data = std::move(
                 received_mortar_data->second.boundary_correction_data.value());
             boundary_data_history->at(mortar_id).remote().insert(

src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp

@@ -844,6 +844,7 @@ void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
               auto& local_mortar_data = mortar_data->at(mortar_id).local();
               local_mortar_data.face_normal_magnitude = face_normal_magnitude;
               if (using_gauss_points) {
+                local_mortar_data.volume_mesh = volume_mesh;
                 local_mortar_data.volume_det_inv_jacobian =
                     volume_det_inv_jacobian;
                 local_mortar_data.face_det_jacobian = face_det_jacobian;

src/Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp

@@ -253,6 +253,7 @@ void internal_mortar_data_impl(
         // Can use the local_mortar_data
         auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
         local_mortar.face_mesh = face_mesh;
+        local_mortar.mortar_mesh = mortar_mesh;
         // If this is the first time, initialize the data. If we don't do this,
         // then the DataVector will be non-owning which we don't want
         if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
@@ -294,6 +295,7 @@ void internal_mortar_data_impl(
       if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
         auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
         local_mortar.face_mesh = face_mesh;
+        local_mortar.mortar_mesh = mortar_mesh;
         // If this is the first time, initialize the data. If we don't do this,
         // then the DataVector will be non-owning which we don't want
         if (UNLIKELY(not local_mortar.mortar_data.has_value())) {

src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp

@@ -72,6 +72,98 @@ mortars_apply_impl(const std::vector<std::array<size_t, Dim>>& initial_extents,
                    Spectral::Quadrature quadrature, const Element<Dim>& element,
                    const TimeStepId& next_temporal_id,
                    const Mesh<Dim>& volume_mesh);
+
+template <size_t Dim, typename MortarDataHistoryType>
+void p_project(
+    const gsl::not_null<
+        ::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
+    /* mortar_data */,
+    const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
+    const gsl::not_null<
+        ::dg::MortarMap<Dim, std::array<Spectral::MortarSize, Dim - 1>>*>
+    /* mortar_size */,
+    const gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*>
+    /* mortar_next_temporal_id */,
+    const gsl::not_null<
+        DirectionMap<Dim, std::optional<Variables<tmpl::list<
+                              evolution::dg::Tags::MagnitudeOfNormal,
+                              evolution::dg::Tags::NormalCovector<Dim>>>>>*>
+        normal_covector_and_magnitude,
+    const gsl::not_null<MortarDataHistoryType*> mortar_data_history,
+    const Mesh<Dim>& new_mesh, const Element<Dim>& new_element,
+    const std::unordered_map<ElementId<Dim>, amr::Info<Dim>>& neighbor_info,
+    const std::pair<Mesh<Dim>, Element<Dim>>& old_mesh_and_element) {
+  const auto& [old_mesh, old_element] = old_mesh_and_element;
+  ASSERT(old_element.id() == new_element.id(),
+         "p-refinement should not have changed the element id");
+
+  const bool mesh_changed = old_mesh != new_mesh;
+
+  for (const auto& [direction, neighbors] : new_element.neighbors()) {
+    const auto sliced_away_dimension = direction.dimension();
+    const auto old_face_mesh = old_mesh.slice_away(sliced_away_dimension);
+    const auto new_face_mesh = new_mesh.slice_away(sliced_away_dimension);
+    const bool face_mesh_changed = old_face_mesh != new_face_mesh;
+    if (face_mesh_changed) {
+      (*normal_covector_and_magnitude)[direction] = std::nullopt;
+    }
+    for (const auto& neighbor : neighbors) {
+      const DirectionalId<Dim> mortar_id{direction, neighbor};
+      if (mortar_mesh->contains(mortar_id)) {
+        const auto new_neighbor_mesh = neighbors.orientation().inverse_map()(
+            neighbor_info.at(neighbor).new_mesh);
+        const auto& old_mortar_mesh = mortar_mesh->at(mortar_id);
+        auto new_mortar_mesh = ::dg::mortar_mesh(
+            new_mesh.slice_away(direction.dimension()),
+            new_neighbor_mesh.slice_away(direction.dimension()));
+        const bool mortar_mesh_changed = old_mortar_mesh != new_mortar_mesh;
+
+        if (mortar_mesh_changed or mesh_changed) {
+          // mortar_data does not need projecting as it has already been used
+          // and will be resized automatically
+          // mortar_size does not change as the mortar has not changed
+          // next_temporal_id does not change as the mortar has not changed
+          if (not mortar_data_history->empty()) {
+            auto& boundary_history = mortar_data_history->at(mortar_id);
+            auto local_history = boundary_history.local();
+            auto remote_history = boundary_history.remote();
+            const auto project_local_boundary_data =
+                [&new_mortar_mesh, &new_face_mesh, &new_mesh](
+                    const TimeStepId& /* id */,
+                    const gsl::not_null<::evolution::dg::MortarData<Dim>*>
+                        mortar_data) {
+                  p_project(mortar_data, new_mortar_mesh, new_face_mesh,
+                            new_mesh);
+                };
+            local_history.for_each(project_local_boundary_data);
+            const auto project_remote_boundary_data =
+                [&new_mortar_mesh](
+                    const TimeStepId& /* id */,
+                    const gsl::not_null<::evolution::dg::MortarData<Dim>*>
+                        mortar_data) {
+                  p_project_only_mortar_data(mortar_data, new_mortar_mesh);
+                };
+            remote_history.for_each(project_remote_boundary_data);
+            boundary_history.clear_coupling_cache();
+          }
+          mortar_mesh->at(mortar_id) = std::move(new_mortar_mesh);
+        }
+      } else {
+        ERROR("h-refinement not implemented yet");
+      }
+    }
+  }
+
+  for (const auto& direction : new_element.external_boundaries()) {
+    const auto sliced_away_dimension = direction.dimension();
+    const auto old_face_mesh = old_mesh.slice_away(sliced_away_dimension);
+    const auto new_face_mesh = new_mesh.slice_away(sliced_away_dimension);
+    const bool face_mesh_changed = old_face_mesh != new_face_mesh;
+    if (face_mesh_changed) {
+      (*normal_covector_and_magnitude)[direction] = std::nullopt;
+    }
+  }
+}
 }  // namespace detail
 
 /*!
@@ -192,48 +284,21 @@ struct ProjectMortars : tt::ConformsTo<amr::protocols::Projector> {
       const gsl::not_null<
           ::dg::MortarMap<dim, std::array<Spectral::MortarSize, dim - 1>>*>
           mortar_size,
-      const gsl::not_null<
-          ::dg::MortarMap<dim, TimeStepId>*> /*mortar_next_temporal_id*/,
+      const gsl::not_null<::dg::MortarMap<dim, TimeStepId>*>
+          mortar_next_temporal_id,
       const gsl::not_null<
           DirectionMap<dim, std::optional<Variables<tmpl::list<
                                 evolution::dg::Tags::MagnitudeOfNormal,
                                 evolution::dg::Tags::NormalCovector<dim>>>>>*>
           normal_covector_and_magnitude,
-      const gsl::not_null<mortar_data_history_type*>
-      /*mortar_data_history*/,
+      const gsl::not_null<mortar_data_history_type*> mortar_data_history,
       const Mesh<dim>& new_mesh, const Element<dim>& new_element,
       const std::unordered_map<ElementId<dim>, amr::Info<dim>>& neighbor_info,
-      const std::pair<Mesh<dim>, Element<dim>>& /*old_mesh_and_element*/) {
-    if (Metavariables::local_time_stepping) {
-      ERROR("AMR with local time-stepping is not yet supported");
-    }
-
-    mortar_data->clear();
-    mortar_mesh->clear();
-    mortar_size->clear();
-    // mortar_next_temporal_id is not changed, but this will break when
-    // h-refinement is enabled and the neighbors are no longer the same
-    for (const auto& [direction, neighbors] : new_element.neighbors()) {
-      (*normal_covector_and_magnitude)[direction] = std::nullopt;
-      for (const auto& neighbor : neighbors) {
-        const DirectionalId<dim> mortar_id{direction, neighbor};
-        mortar_data->emplace(mortar_id, MortarDataHolder<dim>{});
-        const auto new_neighbor_mesh = neighbors.orientation().inverse_map()(
-            neighbor_info.at(neighbor).new_mesh);
-        mortar_mesh->emplace(
-            mortar_id,
-            ::dg::mortar_mesh(
-                new_mesh.slice_away(direction.dimension()),
-                new_neighbor_mesh.slice_away(direction.dimension())));
-        mortar_size->emplace(
-            mortar_id,
-            ::dg::mortar_size(new_element.id(), neighbor, direction.dimension(),
-                              neighbors.orientation()));
-      }
-    }
-    for (const auto& direction : new_element.external_boundaries()) {
-      (*normal_covector_and_magnitude)[direction] = std::nullopt;
-    }
+      const std::pair<Mesh<dim>, Element<dim>>& old_mesh_and_element) {
+    detail::p_project(mortar_data, mortar_mesh, mortar_size,
+                      mortar_next_temporal_id, normal_covector_and_magnitude,
+                      mortar_data_history, new_mesh, new_element, neighbor_info,
+                      old_mesh_and_element);
   }
 
   template <typename... Tags>