Add ability to p-project MortarData

Use it in ProjectMortars to project the BoundaryHistory for LTS.

src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp

@@ -72,6 +72,96 @@ 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");
+
+  std::unordered_set<DirectionalId<Dim>> old_mortar_ids;
+  for (const auto& [mortar_id, _] : *mortar_mesh) {
+    old_mortar_ids.emplace(mortar_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 (old_mortar_ids.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_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_boundary_data);
+            remote_history.for_each(project_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 +282,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>

src/Evolution/DiscontinuousGalerkin/MortarData.cpp

@@ -12,9 +12,11 @@
 #include <utility>
 #include <vector>
 
+#include "DataStructures/ApplyMatrices.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Projection.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Serialization/PupStlCpp17.hpp"
@@ -31,6 +33,54 @@ void MortarData<Dim>::pup(PUP::er& p) {
   p | volume_mesh;
 }
 
+template <size_t Dim>
+void p_project(
+    const gsl::not_null<::evolution::dg::MortarData<Dim>*> mortar_data,
+    const Mesh<Dim - 1>& new_mortar_mesh, const Mesh<Dim - 1>& new_face_mesh,
+    const Mesh<Dim>& new_volume_mesh) {
+  // nothing needs to be done in 1D as mortars/faces are a single point...
+  if constexpr (Dim > 1) {
+    if (mortar_data->mortar_data.has_value()) {
+      const auto& old_mortar_mesh = mortar_data->mortar_mesh.value();
+      if (old_mortar_mesh != new_mortar_mesh) {
+        const auto mortar_projection_matrices =
+            Spectral::p_projection_matrices(old_mortar_mesh, new_mortar_mesh);
+        DataVector& vars = mortar_data->mortar_data.value();
+        vars = apply_matrices(mortar_projection_matrices, vars,
+                              old_mortar_mesh.extents());
+        mortar_data->mortar_mesh = new_mortar_mesh;
+      }
+    }
+    if (mortar_data->face_normal_magnitude.has_value()) {
+      const auto& old_face_mesh = mortar_data->face_mesh.value();
+      if (old_face_mesh != new_face_mesh) {
+        const auto face_projection_matrices =
+            Spectral::p_projection_matrices(old_face_mesh, new_face_mesh);
+        DataVector& n = get(mortar_data->face_normal_magnitude.value());
+        n = apply_matrices(face_projection_matrices, n,
+                           old_face_mesh.extents());
+        if (mortar_data->face_det_jacobian.has_value()) {
+          DataVector& det_j = get(mortar_data->face_det_jacobian.value());
+          det_j = apply_matrices(face_projection_matrices, det_j,
+                                 old_face_mesh.extents());
+        }
+        mortar_data->face_mesh = new_face_mesh;
+      }
+    }
+  }
+  if (mortar_data->volume_det_inv_jacobian.has_value()) {
+    const auto& old_volume_mesh = mortar_data->volume_mesh.value();
+    if (old_volume_mesh != new_volume_mesh) {
+      const auto volume_projection_matrices =
+          Spectral::p_projection_matrices(old_volume_mesh, new_volume_mesh);
+      DataVector& det_inv_j = get(mortar_data->volume_det_inv_jacobian.value());
+      det_inv_j = apply_matrices(volume_projection_matrices, det_inv_j,
+                                 old_volume_mesh.extents());
+      mortar_data->volume_mesh = new_volume_mesh;
+    }
+  }
+}
+
 template <size_t Dim>
 bool operator==(const MortarData<Dim>& lhs, const MortarData<Dim>& rhs) {
   return lhs.mortar_data == rhs.mortar_data and
@@ -60,13 +110,19 @@ std::ostream& operator<<(std::ostream& os, const MortarData<Dim>& mortar_data) {
 
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
-#define INSTANTIATION(r, data)                                \
-  template class MortarData<DIM(data)>;                       \
-  template bool operator==(const MortarData<DIM(data)>& lhs,  \
-                           const MortarData<DIM(data)>& rhs); \
-  template bool operator!=(const MortarData<DIM(data)>& lhs,  \
-                           const MortarData<DIM(data)>& rhs); \
-  template std::ostream& operator<<(std::ostream& os,         \
+#define INSTANTIATION(r, data)                                     \
+  template class MortarData<DIM(data)>;                            \
+  template void p_project(                                         \
+      const gsl::not_null<::evolution::dg::MortarData<DIM(data)>*> \
+          mortar_data,                                             \
+      const Mesh<DIM(data) - 1>& new_mortar_mesh,                  \
+      const Mesh<DIM(data) - 1>& new_face_mesh,                    \
+      const Mesh<DIM(data)>& volume_mesh);                         \
+  template bool operator==(const MortarData<DIM(data)>& lhs,       \
+                           const MortarData<DIM(data)>& rhs);      \
+  template bool operator!=(const MortarData<DIM(data)>& lhs,       \
+                           const MortarData<DIM(data)>& rhs);      \
+  template std::ostream& operator<<(std::ostream& os,              \
                                     const MortarData<DIM(data)>& mortar_data);
 
 GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))

src/Evolution/DiscontinuousGalerkin/MortarData.hpp

@@ -72,6 +72,19 @@ struct MortarData {
   void pup(PUP::er& p);
 };
 
+/// \brief Projects the mortar data when p-refined
+///
+/// \details only updates the stored mesh if the corresponding data exists
+///
+/// \note The DG-subcell code stores the face mesh in the MortarData even when
+/// the geometric data are not used.  Currently projection of MortarData is only
+///  needed for local time-stepping which is not yet supported for DG-subcell.
+template <size_t Dim>
+void p_project(gsl::not_null<::evolution::dg::MortarData<Dim>*> mortar_data,
+               const Mesh<Dim - 1>& new_mortar_mesh,
+               const Mesh<Dim - 1>& new_face_mesh,
+               const Mesh<Dim>& new_volume_mesh);
+
 template <size_t Dim>
 bool operator==(const MortarData<Dim>& lhs, const MortarData<Dim>& rhs);
 

tests/Unit/Evolution/DiscontinuousGalerkin/CMakeLists.txt

@@ -29,6 +29,7 @@ add_test_library(${LIBRARY} "${LIBRARY_SOURCES}")
 target_link_libraries(
   ${LIBRARY}
   PRIVATE
+  Amr
   Boost::boost
   CoordinateMaps
   DataStructures