Support domain deformations in elliptic solver

This enables horizon-conforming (ellipsoidal) excision surfaces for spinning BH initial data. It will also allow for surface-fitting coordinates for BNS initial data.
sxs-collaboration · Oct 19, 2023 · 1cc38e4 · 1cc38e4
1 parent 598360a
commit 1cc38e4
Show file tree

Hide file tree

Showing 18 changed files with 92 additions and 41 deletions.
diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp
@@ -17,6 +17,7 @@
 #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
 #include "Domain/Creators/Tags/InitialExtents.hpp"
 #include "Domain/FaceNormal.hpp"
+#include "Domain/FunctionsOfTime/Tags.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionMap.hpp"
 #include "Domain/Structure/Element.hpp"
@@ -117,7 +118,8 @@ struct InitializeFacesMortarsAndBackground {
       db::mutate_apply<typename InitializeFacesAndMortars::return_tags,
                        typename InitializeFacesAndMortars::argument_tags>(
           InitializeFacesAndMortars{}, make_not_null(&box),
-          db::get<domain::Tags::InitialExtents<Dim>>(box), background,
+          db::get<domain::Tags::InitialExtents<Dim>>(box),
+          db::get<domain::Tags::FunctionsOfTime>(box), background,
           background_classes{});
       // Initialize background fields
       db::mutate_apply<typename InitializeBackground::return_tags,
@@ -129,7 +131,8 @@ struct InitializeFacesMortarsAndBackground {
       db::mutate_apply<typename InitializeFacesAndMortars::return_tags,
                        typename InitializeFacesAndMortars::argument_tags>(
           InitializeFacesAndMortars{}, make_not_null(&box),
-          db::get<domain::Tags::InitialExtents<Dim>>(box));
+          db::get<domain::Tags::InitialExtents<Dim>>(box),
+          db::get<domain::Tags::FunctionsOfTime>(box));
     }
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }

diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp
@@ -10,6 +10,7 @@
 
 #include "DataStructures/DataBox/DataBox.hpp"
 #include "Domain/Creators/Tags/Domain.hpp"
+#include "Domain/Creators/Tags/FunctionsOfTime.hpp"
 #include "Domain/Creators/Tags/InitialExtents.hpp"
 #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
@@ -40,6 +41,8 @@ namespace elliptic::dg::Actions {
  *   - `domain::Tags::Domain<Dim, Frame::Inertial>`
  *   - `domain::Tags::InitialExtents<Dim>`
  *   - `domain::Tags::InitialRefinementLevels<Dim>`
+ *   - `domain::Tags::FunctionsOfTime` (Parameters for maps that distort the
+ *     domain. These are always evaluated at t=0.)
  * - Adds:
  *   - `domain::Tags::Mesh<Dim>`
  *   - `domain::Tags::Element<Dim>`
@@ -61,6 +64,8 @@ struct InitializeDomain {
                  domain::Tags::InitialRefinementLevels<Dim>>;
   using simple_tags = typename InitializeGeometry::return_tags;
   using compute_tags = tmpl::list<>;
+  using const_global_cache_tags =
+      tmpl::list<domain::Tags::FunctionsOfTimeInitialize>;
 
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
             typename ActionList, typename ParallelComponent>

diff --git a/src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt b/src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt
@@ -31,6 +31,7 @@ target_link_libraries(
   Domain
   DomainStructure
   ErrorHandling
+  FunctionsOfTime
   LinearOperators
   Spectral
   Utilities

diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp
@@ -50,29 +50,26 @@ void InitializeGeometry<Dim>::operator()(
     const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
     const std::vector<std::array<size_t, Dim>>& initial_extents,
     const std::vector<std::array<size_t, Dim>>& initial_refinement,
-    const Domain<Dim>& domain, const ElementId<Dim>& element_id) const {
+    const Domain<Dim>& domain,
+    const domain::FunctionsOfTimeMap& functions_of_time,
+    const ElementId<Dim>& element_id) const {
   // Mesh
   const auto quadrature = Spectral::Quadrature::GaussLobatto;
   *mesh = domain::Initialization::create_initial_mesh(initial_extents,
                                                       element_id, quadrature);
   // Element
   const auto& block = domain.blocks()[element_id.block_id()];
-  if (block.is_time_dependent()) {
-    ERROR_NO_TRACE(
-        "The InitializeDomain action is for elliptic systems which do not have "
-        "any time-dependence, but the domain creator has set up the domain to "
-        "have time-dependence.");
-  }
   *element = domain::Initialization::create_initial_element(element_id, block,
                                                             initial_refinement);
   // Element map
-  *element_map = ElementMap<Dim, Frame::Inertial>{
-      element_id, block.stationary_map().get_clone()};
+  *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
   // Coordinates
   *logical_coords = logical_coordinates(*mesh);
-  *inertial_coords = element_map->operator()(*logical_coords);
+  *inertial_coords =
+      element_map->operator()(*logical_coords, 0., functions_of_time);
   // Jacobian
-  *inv_jacobian = element_map->inv_jacobian(*logical_coords);
+  *inv_jacobian =
+      element_map->inv_jacobian(*logical_coords, 0., functions_of_time);
   *det_inv_jacobian = determinant(*inv_jacobian);
 }
 

diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp
@@ -21,6 +21,8 @@
 #include "Domain/Domain.hpp"
 #include "Domain/ElementMap.hpp"
 #include "Domain/FaceNormal.hpp"
+#include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
+#include "Domain/FunctionsOfTime/Tags.hpp"
 #include "Domain/InterfaceLogicalCoordinates.hpp"
 #include "Domain/Structure/CreateInitialMesh.hpp"
 #include "Domain/Structure/Direction.hpp"
@@ -57,9 +59,10 @@ struct InitializeGeometry {
       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
                                     Frame::Inertial>,
       domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>;
-  using argument_tags = tmpl::list<domain::Tags::InitialExtents<Dim>,
-                                   domain::Tags::InitialRefinementLevels<Dim>,
-                                   domain::Tags::Domain<Dim>>;
+  using argument_tags =
+      tmpl::list<domain::Tags::InitialExtents<Dim>,
+                 domain::Tags::InitialRefinementLevels<Dim>,
+                 domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
   void operator()(
       gsl::not_null<Mesh<Dim>*> mesh, gsl::not_null<Element<Dim>*> element,
       gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
@@ -72,7 +75,9 @@ struct InitializeGeometry {
       gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
       const std::vector<std::array<size_t, Dim>>& initial_extents,
       const std::vector<std::array<size_t, Dim>>& initial_refinement,
-      const Domain<Dim>& domain, const ElementId<Dim>& element_id) const;
+      const Domain<Dim>& domain,
+      const domain::FunctionsOfTimeMap& functions_of_time,
+      const ElementId<Dim>& element_id) const;
 };
 
 namespace detail {
@@ -157,6 +162,7 @@ struct InitializeFacesAndMortars {
       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
                             Frame::Inertial>& inv_jacobian,
       const std::vector<std::array<size_t, Dim>>& initial_extents,
+      const domain::FunctionsOfTimeMap& functions_of_time,
       const Background& background = std::nullptr_t{},
       tmpl::list<BackgroundClasses...> /*meta*/ = tmpl::list<>{}) const {
     static_assert(std::is_same_v<InvMetricTag, void> or
@@ -191,10 +197,14 @@ struct InitializeFacesAndMortars {
       const auto face_logical_coords =
           interface_logical_coordinates(face_mesh, direction);
       auto& face_inertial_coords = (*faces_inertial_coords)[direction];
-      face_inertial_coords = element_map.operator()(face_logical_coords);
+      face_inertial_coords =
+          element_map(face_logical_coords, 0., functions_of_time);
       auto& face_normal = (*face_normals)[direction];
       auto& face_normal_magnitude = (*face_normal_magnitudes)[direction];
-      face_normal = unnormalized_face_normal(face_mesh, element_map, direction);
+      unnormalized_face_normal(
+          make_not_null(&face_normal), face_mesh,
+          element_map.inv_jacobian(face_logical_coords, 0., functions_of_time),
+          direction);
       if constexpr (std::is_same_v<InvMetricTag, void>) {
         magnitude(make_not_null(&face_normal_magnitude), face_normal);
       } else {
@@ -206,7 +216,8 @@ struct InitializeFacesAndMortars {
         face_normal.get(d) /= get(face_normal_magnitude);
       }
       get((*face_jacobians)[direction]) =
-          get(determinant(element_map.jacobian(face_logical_coords))) *
+          get(determinant(element_map.jacobian(face_logical_coords, 0.,
+                                               functions_of_time))) *
           get(face_normal_magnitude);
     }
     // Compute the Jacobian derivative numerically, because our coordinate maps
@@ -237,16 +248,16 @@ struct InitializeFacesAndMortars {
           const auto mortar_logical_coords = detail::mortar_logical_coordinates(
               mortar_mesh, mortar_size, direction);
           auto& mortar_jacobian = (*mortar_jacobians)[mortar_id];
-          mortar_jacobian =
-              determinant(element_map.jacobian(mortar_logical_coords));
+          mortar_jacobian = determinant(element_map.jacobian(
+              mortar_logical_coords, 0., functions_of_time));
           // These factors of two account for the mortar size
           for (const auto& mortar_size_i : mortar_size) {
             if (mortar_size_i != Spectral::MortarSize::Full) {
               get(mortar_jacobian) *= 0.5;
             }
           }
-          const auto inv_jacobian_on_mortar =
-              element_map.inv_jacobian(mortar_logical_coords);
+          const auto inv_jacobian_on_mortar = element_map.inv_jacobian(
+              mortar_logical_coords, 0., functions_of_time);
           const auto unnormalized_mortar_normal = unnormalized_face_normal(
               mortar_mesh, inv_jacobian_on_mortar, direction);
           Scalar<DataVector> mortar_normal_magnitude{};
@@ -255,7 +266,7 @@ struct InitializeFacesAndMortars {
                       unnormalized_mortar_normal);
           } else {
             const auto mortar_inertial_coords =
-                element_map(mortar_logical_coords);
+                element_map(mortar_logical_coords, 0., functions_of_time);
             const auto inv_metric_on_mortar =
                 get_inv_metric(mortar_inertial_coords);
             magnitude(make_not_null(&mortar_normal_magnitude),

diff --git a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/CMakeLists.txt b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/CMakeLists.txt
@@ -26,6 +26,7 @@ target_link_libraries(
   DataStructures
   DomainStructure
   DiscontinuousGalerkin
+  FunctionsOfTime
   ParallelSchwarz
   Spectral
   Utilities

diff --git a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp
@@ -25,6 +25,7 @@
 #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
 #include "Domain/ElementMap.hpp"
 #include "Domain/FaceNormal.hpp"
+#include "Domain/FunctionsOfTime/Tags.hpp"
 #include "Domain/Structure/CreateInitialMesh.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionMap.hpp"
@@ -217,6 +218,7 @@ struct InitializeSubdomain {
               tmpl::list<>>(InitializeFacesAndMortars{}, make_not_null(&box),
                             overlap_id,
                             db::get<domain::Tags::InitialExtents<Dim>>(box),
+                            db::get<domain::Tags::FunctionsOfTime>(box),
                             background, background_classes{});
           // Background fields
           initialize_background_fields(make_not_null(&box), overlap_id);
@@ -230,7 +232,8 @@ struct InitializeSubdomain {
                   typename InitializeFacesAndMortars::argument_tags>,
               tmpl::list<>>(InitializeFacesAndMortars{}, make_not_null(&box),
                             overlap_id,
-                            db::get<domain::Tags::InitialExtents<Dim>>(box));
+                            db::get<domain::Tags::InitialExtents<Dim>>(box),
+                            db::get<domain::Tags::FunctionsOfTime>(box));
         }
         // Faces on the other side of the overlapped element's mortars
         initialize_remote_faces(make_not_null(&box), overlap_id);
@@ -317,6 +320,8 @@ struct InitializeSubdomain {
     const auto& element =
         db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
     const auto& domain = db::get<domain::Tags::Domain<Dim>>(*box);
+    const auto& functions_of_time =
+        db::get<domain::Tags::FunctionsOfTime>(*box);
     const auto& neighbor_meshes = db::get<overlaps_tag<
         elliptic::dg::subdomain_operator::Tags::NeighborMortars<
             domain::Tags::Mesh<Dim>, Dim>>>(*box)
@@ -329,19 +334,23 @@ struct InitializeSubdomain {
         const auto neighbor_face_mesh =
             neighbor_meshes.at(mortar_id).slice_away(
                 direction_from_neighbor.dimension());
+        const auto neighbor_face_logical_coords = interface_logical_coordinates(
+            neighbor_face_mesh, direction_from_neighbor);
         const auto& neighbor_block = domain.blocks()[neighbor_id.block_id()];
-        ElementMap<Dim, Frame::Inertial> neighbor_element_map{
-            neighbor_id, neighbor_block.stationary_map().get_clone()};
+        const ElementMap<Dim, Frame::Inertial> neighbor_element_map{
+            neighbor_id, neighbor_block};
         const auto neighbor_face_normal = unnormalized_face_normal(
-            neighbor_face_mesh, neighbor_element_map, direction_from_neighbor);
+            neighbor_face_mesh,
+            neighbor_element_map.inv_jacobian(neighbor_face_logical_coords, 0.,
+                                              functions_of_time),
+            direction_from_neighbor);
         using neighbor_face_normal_magnitudes_tag = overlaps_tag<
             elliptic::dg::subdomain_operator::Tags::NeighborMortars<
                 domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>, Dim>>;
         if constexpr (is_curved) {
           const auto& background = db::get<BackgroundTag>(*box);
-          const auto neighbor_face_inertial_coords =
-              neighbor_element_map(interface_logical_coordinates(
-                  neighbor_face_mesh, direction_from_neighbor));
+          const auto neighbor_face_inertial_coords = neighbor_element_map(
+              neighbor_face_logical_coords, 0., functions_of_time);
           const auto inv_metric_on_face = get<typename System::inv_metric_tag>(
               elliptic::util::get_analytic_data<
                   tmpl::list<typename System::inv_metric_tag>>(

diff --git a/src/Elliptic/Executables/Elasticity/CMakeLists.txt b/src/Elliptic/Executables/Elasticity/CMakeLists.txt
@@ -18,6 +18,7 @@ set(LIBS_TO_LINK
   EllipticSubdomainPreconditioners
   Events
   EventsAndTriggers
+  FunctionsOfTime
   Informer
   LinearOperators
   Observer

diff --git a/src/Elliptic/Executables/Elasticity/SolveElasticity.cpp b/src/Elliptic/Executables/Elasticity/SolveElasticity.cpp
@@ -6,6 +6,7 @@
 #include <vector>
 
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
 #include "Elliptic/SubdomainPreconditioners/RegisterDerived.hpp"
 #include "Parallel/CharmMain.tpp"
 #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
@@ -17,6 +18,7 @@ extern "C" void CkRegisterMainModule() {
   Parallel::charmxx::register_main_module<metavariables>();
   Parallel::charmxx::register_init_node_and_proc(
       {&domain::creators::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
        &register_derived_classes_with_charm<
            metavariables::schwarz_smoother::subdomain_solver>,
        &elliptic::subdomain_preconditioners::register_derived_with_charm,

diff --git a/src/Elliptic/Executables/Poisson/CMakeLists.txt b/src/Elliptic/Executables/Poisson/CMakeLists.txt
@@ -12,6 +12,7 @@ set(LIBS_TO_LINK
   EllipticDgSubdomainOperator
   Events
   EventsAndTriggers
+  FunctionsOfTime
   Informer
   LinearOperators
   MathFunctions

diff --git a/src/Elliptic/Executables/Poisson/SolvePoisson.cpp b/src/Elliptic/Executables/Poisson/SolvePoisson.cpp
@@ -6,6 +6,7 @@
 #include <vector>
 
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
 #include "Parallel/CharmMain.tpp"
 #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
 
@@ -16,6 +17,7 @@ extern "C" void CkRegisterMainModule() {
   Parallel::charmxx::register_main_module<metavariables>();
   Parallel::charmxx::register_init_node_and_proc(
       {&domain::creators::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
        &register_derived_classes_with_charm<
            metavariables::schwarz_smoother::subdomain_solver>,
        &register_factory_classes_with_charm<metavariables>},

diff --git a/src/Elliptic/Executables/Punctures/CMakeLists.txt b/src/Elliptic/Executables/Punctures/CMakeLists.txt
@@ -26,6 +26,7 @@ target_link_libraries(
   ErrorHandling
   Events
   EventsAndTriggers
+  FunctionsOfTime
   Informer
   Initialization
   LinearOperators

diff --git a/src/Elliptic/Executables/Punctures/SolvePunctures.cpp b/src/Elliptic/Executables/Punctures/SolvePunctures.cpp
@@ -6,6 +6,7 @@
 #include <vector>
 
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
 #include "Elliptic/SubdomainPreconditioners/RegisterDerived.hpp"
 #include "Parallel/CharmMain.tpp"
 #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
@@ -14,6 +15,7 @@ extern "C" void CkRegisterMainModule() {
   Parallel::charmxx::register_main_module<Metavariables>();
   Parallel::charmxx::register_init_node_and_proc(
       {&domain::creators::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
        &register_derived_classes_with_charm<
            Metavariables::solver::schwarz_smoother::subdomain_solver>,
        &elliptic::subdomain_preconditioners::register_derived_with_charm,

diff --git a/src/Elliptic/Executables/Xcts/CMakeLists.txt b/src/Elliptic/Executables/Xcts/CMakeLists.txt
@@ -26,6 +26,7 @@ target_link_libraries(
   ErrorHandling
   Events
   EventsAndTriggers
+  FunctionsOfTime
   Hydro
   Informer
   Initialization

diff --git a/src/Elliptic/Executables/Xcts/SolveXcts.cpp b/src/Elliptic/Executables/Xcts/SolveXcts.cpp
@@ -6,6 +6,7 @@
 #include <vector>
 
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
 #include "Elliptic/SubdomainPreconditioners/RegisterDerived.hpp"
 #include "Parallel/CharmMain.tpp"
 #include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp"
@@ -14,9 +15,8 @@
 extern "C" void CkRegisterMainModule() {
   Parallel::charmxx::register_main_module<Metavariables>();
   Parallel::charmxx::register_init_node_and_proc(
-      {&setup_error_handling, &setup_memory_allocation_failure_reporting,
-       &disable_openblas_multithreading,
-       &domain::creators::register_derived_with_charm,
+      {&domain::creators::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
        &register_derived_classes_with_charm<
            Metavariables::solver::schwarz_smoother::subdomain_solver>,
        &elliptic::subdomain_preconditioners::register_derived_with_charm,

diff --git a/src/Executables/FindHorizons/CMakeLists.txt b/src/Executables/FindHorizons/CMakeLists.txt
@@ -7,6 +7,7 @@ set(LIBS_TO_LINK
   CoordinateMaps
   DomainCreators
   EllipticDg
+  FunctionsOfTime
   GeneralRelativity
   GrSurfaces
   Informer

diff --git a/src/Executables/FindHorizons/FindHorizons.cpp b/src/Executables/FindHorizons/FindHorizons.cpp
@@ -6,6 +6,7 @@
 #include <vector>
 
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
 #include "Parallel/CharmMain.tpp"
 #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
 
@@ -16,6 +17,7 @@ extern "C" void CkRegisterMainModule() {
   Parallel::charmxx::register_main_module<metavariables>();
   Parallel::charmxx::register_init_node_and_proc(
       {&domain::creators::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
        &register_factory_classes_with_charm<metavariables>},
       {});
 }