Instantiate inertial->distorted for strahlkorper_in_different_frame

.clang-tidy

@@ -34,6 +34,8 @@ Checks: '*,
 -cppcoreguidelines-macro-usage,
 # public and protected member variables are fine,
 -cppcoreguidelines-non-private-member-variables-in-classes,
+# We do not use gsl::owner,
+-cppcoreguidelines-owning-memory,
 -fuchsia-*,
 # defaulting virtual functions in CoordinateMap,
 -google-default-arguments,

src/Domain/StrahlkorperTransformations.cpp

@@ -45,10 +45,11 @@ void coords_to_different_frame(
         std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
         functions_of_time,
     const double time) {
-  // If we ever want a source in a different frame than Grid, we
-  // need to add 'if constexpr's for that case.
-  static_assert(std::is_same_v<SrcFrame, ::Frame::Grid>,
-                "Source frame must currently be Grid frame");
+  // If ever additional cases besides Grid->Inertial, Grid->Distorted,
+  // and Inertial->Distorted are needed, add if constexprs below
+  static_assert(std::is_same_v<SrcFrame, ::Frame::Grid> or
+                    std::is_same_v<SrcFrame, ::Frame::Inertial>,
+                "Source frame must currently be Grid frame or Inertial frame");
   const auto block_logical_coords = block_logical_coordinates(
       domain, src_cartesian_coords, time, functions_of_time);
 
@@ -67,18 +68,36 @@ void coords_to_different_frame(
     const auto& block_id_and_coords = block_logical_coords[s].value();
     const auto& block = domain.blocks()[block_id_and_coords.id.get_index()];
 
-    if constexpr (std::is_same_v<DestFrame, ::Frame::Distorted>) {
+    if constexpr (std::is_same_v<DestFrame, ::Frame::Distorted> and
+                  std::is_same_v<SrcFrame, ::Frame::Grid>) {
       if (not block.has_distorted_frame()) {
         throw std::runtime_error(
             "Strahlkorper lies outside of distorted-frame region");
       }
       const auto& grid_to_distorted_map =
           block.moving_mesh_grid_to_distorted_map();
       x_dest = grid_to_distorted_map(x_src, time, functions_of_time);
+    } else if constexpr (std::is_same_v<DestFrame, ::Frame::Distorted> and
+                         std::is_same_v<SrcFrame, ::Frame::Inertial>) {
+      if (not block.has_distorted_frame()) {
+        throw std::runtime_error(
+            "Strahlkorper lies outside of distorted-frame region");
+      }
+      const auto& distorted_to_inertial_map =
+          block.moving_mesh_distorted_to_inertial_map();
+      const auto& inv_point =
+          distorted_to_inertial_map.inverse(x_src, time, functions_of_time);
+      if (inv_point.has_value()) {
+        x_dest = *inv_point;
+      } else {
+        throw std::runtime_error(
+            "Map from Frame::Distorted to Frame::Inertial is not invertible");
+      }
     } else {
-      static_assert(
-          std::is_same_v<DestFrame, ::Frame::Inertial>,
-          "Destination frame must currently be Distorted or Inertial");
+      static_assert(std::is_same_v<DestFrame, ::Frame::Inertial> and
+                        std::is_same_v<SrcFrame, ::Frame::Grid>,
+                    "Source frame -> destination frame must be Grid -> "
+                    "Distorted, Grid -> Inertial, or Inertial -> Distorted");
       const auto& grid_to_inertial_map =
           block.moving_mesh_grid_to_inertial_map();
       x_dest = grid_to_inertial_map(x_src, time, functions_of_time);
@@ -183,13 +202,22 @@ void strahlkorper_in_different_frame(
     const auto& block_id_and_coords = block_logical_coords[0].value();
     const auto& block = domain.blocks()[block_id_and_coords.id.get_index()];
 
-    static_assert(std::is_same_v<SrcFrame, ::Frame::Grid>,
-                  "Source frame must currently be Grid frame");
-    static_assert(std::is_same_v<DestFrame, ::Frame::Inertial>,
-                  "Destination frame must currently be Inertial frame");
-    const auto& x_src =
-        block.moving_mesh_logical_to_grid_map()(block_id_and_coords.data);
-
+    static_assert(
+        std::is_same_v<SrcFrame, ::Frame::Grid> or
+            std::is_same_v<SrcFrame, ::Frame::Inertial>,
+        "Source frame must currently be Grid frame or Inertial frame");
+    static_assert(std::is_same_v<DestFrame, ::Frame::Inertial> or
+                      std::is_same_v<DestFrame, ::Frame::Distorted>,
+                  "Destination frame must currently be Inertial frame or "
+                  "Distorted frame");
+    tnsr::I<double, 3, SrcFrame> x_src{};
+    if constexpr (std::is_same_v<SrcFrame, ::Frame::Grid>) {
+      x_src = block.moving_mesh_logical_to_grid_map()(block_id_and_coords.data);
+    } else {
+      x_src = block.moving_mesh_grid_to_inertial_map()(
+          block.moving_mesh_logical_to_grid_map()(block_id_and_coords.data),
+          time, functions_of_time);
+    }
     const double r_src =
         std::hypot(get<0>(x_src) - center_src[0], get<1>(x_src) - center_src[1],
                    get<2>(x_src) - center_src[2]);
@@ -406,6 +434,20 @@ void strahlkorper_coords_in_different_frame(
 GENERATE_INSTANTIATIONS(INSTANTIATEGENERAL, (::Frame::Grid),
                         (::Frame::Inertial))
 
+// Generate a specific instantiation: inertial -> distorted
+// (needed, e.g., for initializing a binary-black-hole ringdown).
+// Don't just add to GENERATE_INSTANTIATIONS above, to avoid also generating
+// unwanted instantiations, such as inertial -> inertial
+template void strahlkorper_in_different_frame(
+    const gsl::not_null<ylm::Strahlkorper<::Frame::Distorted>*>
+        dest_strahlkorper,
+    const ylm::Strahlkorper<::Frame::Inertial>& src_strahlkorper,
+    const Domain<3>& domain,
+    const std::unordered_map<
+        std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
+        functions_of_time,
+    const double time);
+
 #define INSTANTIATEALIGNED(_, data)                                  \
   template void strahlkorper_in_different_frame_aligned(             \
       const gsl::not_null<ylm::Strahlkorper<DESTFRAME(data)>*>       \

src/Domain/StrahlkorperTransformations.hpp

@@ -37,10 +37,11 @@ class not_null;
 /// requires adding member functions to Block.
 ///
 /// \note strahlkorper_in_different_frame is currently instantiated
-/// only for SrcFrame=Grid and DestFrame=Inertial.  In particular, we
-/// intentionally do not instantiate strahlkorper_in_different_frame
-/// for SrcFrame=Grid and DestFrame=Distorted, because our
-/// Grid->Distorted maps preserve centers and angles, so for
+/// only for i) SrcFrame=Grid and DestFrame=Inertial and ii) SrcFrame=Inertial
+/// and DestFrame=Distorted (necessary for initializing a binary-black-hole
+/// ringdown).  In particular, we intentionally do not instantiate
+/// strahlkorper_in_different_frame for SrcFrame=Grid and DestFrame=Distorted,
+/// because our Grid->Distorted maps preserve centers and angles, so for
 /// Grid->Distorted one should use
 /// strahlkorper_in_different_frame_aligned because it is more
 /// efficient.

tests/Unit/Domain/Test_StrahlkorperTransformations.cpp

@@ -25,18 +25,18 @@
 
 namespace {
 
-template <bool Aligned, typename DestFrame>
+template <bool Aligned, typename SrcFrame, typename DestFrame>
 void test_strahlkorper_in_different_frame() {
   const size_t grid_points_each_dimension = 5;
 
   // Set up a Strahlkorper corresponding to a Schwarzschild hole of
-  // mass 1, in the grid frame.
+  // mass 1, in the source frame.
   // Center the Strahlkorper at (0.03,0.02,0.01) so that we test a
   // nonzero center.
-  const std::array<double, 3> strahlkorper_grid_center = {0.03, 0.02, 0.01};
+  const std::array<double, 3> strahlkorper_src_center = {0.03, 0.02, 0.01};
   const size_t l_max = 8;
-  const ylm::Strahlkorper<Frame::Grid> strahlkorper_grid(
-      l_max, 2.0, strahlkorper_grid_center);
+  const ylm::Strahlkorper<SrcFrame> strahlkorper_src(l_max, 2.0,
+                                                     strahlkorper_src_center);
 
   // Create a Domain.
   // We choose a spherical shell domain extending from radius 1.9M to
@@ -56,15 +56,16 @@ void test_strahlkorper_in_different_frame() {
         std::make_unique<domain::creators::time_dependence::Shape<
             domain::ObjectLabel::None>>(0.0, l_max, 1.0,
                                         std::array<double, 3>{{0.1, 0.2, 0.3}},
-                                        strahlkorper_grid_center, 2.0, 12.0)));
+                                        strahlkorper_src_center, 2.0, 12.0)));
   } else {
     domain_creator.reset(new domain::creators::Sphere(
         1.9, 2.9, domain::creators::Sphere::Excision{}, 1_st,
         grid_points_each_dimension, false, std::nullopt, radial_partitioning,
         radial_distribution, ShellWedges::All,
         std::make_unique<
             domain::creators::time_dependence::UniformTranslation<3>>(
-            0.0, std::array<double, 3>({{0.01, 0.02, 0.03}}))));
+            0.0, std::array<double, 3>({{0.0, 0.0, 0.0}}),
+            std::array<double, 3>({{0.01, 0.02, 0.03}}))));
   }
   Domain<3> domain = domain_creator->create_domain();
   const auto functions_of_time = domain_creator->functions_of_time();
@@ -74,13 +75,13 @@ void test_strahlkorper_in_different_frame() {
   ylm::Strahlkorper<DestFrame> strahlkorper_dest{};
   if constexpr (Aligned) {
     strahlkorper_in_different_frame_aligned(make_not_null(&strahlkorper_dest),
-                                            strahlkorper_grid, domain,
+                                            strahlkorper_src, domain,
                                             functions_of_time, time);
 
   } else {
     strahlkorper_in_different_frame(make_not_null(&strahlkorper_dest),
-                                    strahlkorper_grid, domain,
-                                    functions_of_time, time);
+                                    strahlkorper_src, domain, functions_of_time,
+                                    time);
   }
 
   // Now compare.
@@ -92,13 +93,18 @@ void test_strahlkorper_in_different_frame() {
             2.0, ylm.theta_phi_points(), 1.0,
             std::array<double, 3>{{0.1, 0.2, 0.3}}));
     strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(
-        l_max, l_max, new_radius, strahlkorper_grid_center));
+        l_max, l_max, new_radius, strahlkorper_src_center));
+  } else if constexpr (std::is_same_v<SrcFrame, ::Frame::Inertial> and
+                       std::is_same_v<DestFrame, ::Frame::Distorted>) {
+    strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(
+        l_max, 2.0,
+        {{strahlkorper_src_center[0] - 0.005, strahlkorper_src_center[1] - 0.01,
+          strahlkorper_src_center[2] - 0.015}}));
   } else {
     strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(
         l_max, 2.0,
-        {{strahlkorper_grid_center[0] + 0.005,
-          strahlkorper_grid_center[1] + 0.01,
-          strahlkorper_grid_center[2] + 0.015}}));
+        {{strahlkorper_src_center[0] + 0.005, strahlkorper_src_center[1] + 0.01,
+          strahlkorper_src_center[2] + 0.015}}));
   }
   CHECK_ITERABLE_APPROX(strahlkorper_expected->physical_center(),
                         strahlkorper_dest.physical_center());
@@ -190,9 +196,11 @@ SPECTRE_TEST_CASE("Unit.Domain.StrahlkorperTransformations", "[Unit]") {
   domain::creators::register_derived_with_charm();
   domain::creators::time_dependence::register_derived_with_charm();
   domain::FunctionsOfTime::register_derived_with_charm();
-  test_strahlkorper_in_different_frame<false, Frame::Inertial>();
-  test_strahlkorper_in_different_frame<true, Frame::Inertial>();
-  test_strahlkorper_in_different_frame<true, Frame::Distorted>();
+  test_strahlkorper_in_different_frame<false, Frame::Grid, Frame::Inertial>();
+  test_strahlkorper_in_different_frame<false, Frame::Inertial,
+                                       Frame::Distorted>();
+  test_strahlkorper_in_different_frame<true, Frame::Grid, Frame::Inertial>();
+  test_strahlkorper_in_different_frame<true, Frame::Grid, Frame::Distorted>();
   test_strahlkorper_coords_in_different_frame<true, Frame::Grid>();
   test_strahlkorper_coords_in_different_frame<true, Frame::Distorted>();
   test_strahlkorper_coords_in_different_frame<false, Frame::Distorted>();