Add AMR projector for elliptic domain

sxs-collaboration · Feb 2, 2024 · f8a2d9b · f8a2d9b
1 parent bb8d2cb
commit f8a2d9b
Show file tree

Hide file tree

Showing 2 changed files with 139 additions and 18 deletions.
diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp
@@ -36,6 +36,44 @@
 
 namespace elliptic::dg {
 
+namespace detail{
+template <size_t Dim>
+void initialize_coords_and_jacobians(
+    gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
+        logical_coords,
+    gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
+    gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                                  Frame::Inertial>*>
+        inv_jacobian,
+    gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
+    gsl::not_null<Scalar<DataVector>*> det_jacobian,
+    gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                                  Frame::Inertial>*>
+        det_times_inv_jacobian,
+        const Mesh<Dim>& mesh,
+    const ElementMap<Dim, Frame::Inertial>& element_map,
+    const domain::FunctionsOfTimeMap& functions_of_time) {
+  // Coordinates
+  *logical_coords = logical_coordinates(mesh);
+  *inertial_coords =
+      element_map(*logical_coords, 0., functions_of_time);
+  // Jacobian
+  // Note: we can try to use `::dg::metric_identity_jacobian_quantities` here.
+  // When I tried (NV, Dec 2023) the DG residuals diverged on a sphere domain
+  // with a large outer boundary (1e9). This was possibly because no
+  // metric-identity Jacobians were used on faces, though I also tried slicing
+  // the metric-identity Jacobian to the faces and that didn't help.
+  *inv_jacobian =
+      element_map.inv_jacobian(*logical_coords, 0., functions_of_time);
+  *det_inv_jacobian = determinant(*inv_jacobian);
+  get(*det_jacobian) = 1. / get(*det_inv_jacobian);
+  *det_times_inv_jacobian = *inv_jacobian;
+  for (auto& component : *det_times_inv_jacobian) {
+    component *= get(*det_jacobian);
+  }
+    }
+} // namespace detail
+
 template <size_t Dim>
 void InitializeGeometry<Dim>::operator()(
     const gsl::not_null<Mesh<Dim>*> mesh,
@@ -73,24 +111,11 @@ void InitializeGeometry<Dim>::operator()(
                                                             initial_refinement);
   // Element map
   *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
-  // Coordinates
-  *logical_coords = logical_coordinates(*mesh);
-  *inertial_coords =
-      element_map->operator()(*logical_coords, 0., functions_of_time);
-  // Jacobian
-  // Note: we can try to use `::dg::metric_identity_jacobian_quantities` here.
-  // When I tried (NV, Dec 2023) the DG residuals diverged on a sphere domain
-  // with a large outer boundary (1e9). This was possibly because no
-  // metric-identity Jacobians were used on faces, though I also tried slicing
-  // the metric-identity Jacobian to the faces and that didn't help.
-  *inv_jacobian =
-      element_map->inv_jacobian(*logical_coords, 0., functions_of_time);
-  *det_inv_jacobian = determinant(*inv_jacobian);
-  get(*det_jacobian) = 1. / get(*det_inv_jacobian);
-  *det_times_inv_jacobian = *inv_jacobian;
-  for (auto& component : *det_times_inv_jacobian) {
-    component *= get(*det_jacobian);
-  }
+  // Coordinates and Jacobians
+  detail::initialize_coords_and_jacobians(
+      logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
+      det_jacobian, det_times_inv_jacobian, *mesh, *element_map,
+      functions_of_time);
 }
 
 namespace detail {

diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp
@@ -43,14 +43,34 @@
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "NumericalAlgorithms/Spectral/Projection.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
 #include "Utilities/CallWithDynamicType.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
 
 namespace elliptic::dg {
 
+namespace detail {
+template <size_t Dim>
+void initialize_coords_and_jacobians(
+    gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
+        logical_coords,
+    gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
+    gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                                  Frame::Inertial>*>
+        inv_jacobian,
+    gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
+    gsl::not_null<Scalar<DataVector>*> det_jacobian,
+    gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                                  Frame::Inertial>*>
+        det_times_inv_jacobian,
+    const Mesh<Dim>& mesh, const ElementMap<Dim, Frame::Inertial>& element_map,
+    const domain::FunctionsOfTimeMap& functions_of_time);
+}  // namespace detail
+
 /*!
  * \brief Initialize the background-independent geometry for the elliptic DG
  * operator
@@ -102,6 +122,82 @@ struct InitializeGeometry {
       Spectral::Quadrature quadrature, const ElementId<Dim>& element_id) const;
 };
 
+template <size_t Dim>
+struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> {
+  using return_tags = tmpl::list<
+      domain::Tags::ElementMap<Dim>,
+      domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
+      domain::Tags::Coordinates<Dim, Frame::Inertial>,
+      domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
+                                    Frame::Inertial>,
+      domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
+                                        Frame::Inertial>>;
+  using argument_tags =
+      tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
+                 domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
+
+  // p-refinement
+  static void apply(
+      const gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
+      const gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
+          logical_coords,
+      const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
+          inertial_coords,
+      const gsl::not_null<InverseJacobian<
+          DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
+          inv_jacobian,
+      const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
+      const gsl::not_null<Scalar<DataVector>*> det_jacobian,
+      const gsl::not_null<InverseJacobian<
+          DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
+          det_times_inv_jacobian,
+      const Mesh<Dim>& mesh, const Element<Dim>& /*element*/,
+      const Domain<Dim>& /*domain*/,
+      const domain::FunctionsOfTimeMap& functions_of_time,
+      const std::pair<Mesh<Dim>, Element<Dim>>& old_mesh_and_element) {
+    if (mesh == old_mesh_and_element.first) {
+      // Neighbors may have changed, but this element hasn't. Nothing to do.
+      return;
+    }
+    // The element map doesn't change with p-refinement
+    detail::initialize_coords_and_jacobians(
+        logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
+        det_jacobian, det_times_inv_jacobian, mesh, *element_map,
+        functions_of_time);
+  }
+
+  // h-refinement
+  template <typename... ParentOrChildrenItemsType>
+  static void apply(
+      const gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
+      const gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
+          logical_coords,
+      const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
+          inertial_coords,
+      const gsl::not_null<InverseJacobian<
+          DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
+          inv_jacobian,
+      const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
+      const gsl::not_null<Scalar<DataVector>*> det_jacobian,
+      const gsl::not_null<InverseJacobian<
+          DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
+          det_times_inv_jacobian,
+      const Mesh<Dim>& mesh, const Element<Dim>& element,
+      const Domain<Dim>& domain,
+      const domain::FunctionsOfTimeMap& functions_of_time,
+      const ParentOrChildrenItemsType&... /*parent_or_children_items*/) {
+    const auto& element_id = element.id();
+    const auto& block = domain.blocks()[element_id.block_id()];
+    *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
+    detail::initialize_coords_and_jacobians(
+        logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
+        det_jacobian, det_times_inv_jacobian, mesh, *element_map,
+        functions_of_time);
+  }
+};
+
 namespace detail {
 // Compute derivative of the Jacobian numerically
 template <size_t Dim>