Add BoostVelocity option to ShapeMap

src/Domain/Creators/TimeDependence/Shape.cpp

@@ -153,7 +153,8 @@ Shape<Label>::functions_of_time(const std::unordered_map<std::string, double>&
   const DataVector radial_distortion =
       inner_radius_ -
       get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-          inner_radius_, ylm.theta_phi_points(), mass_, spin_));
+          inner_radius_, ylm.theta_phi_points(), mass_, spin_,
+          std::array<double, 3>({0.0, 0.0, 0.0})));
   const auto radial_distortion_coefs = ylm.phys_to_spec(radial_distortion);
   const DataVector zeros =
       make_with_value<DataVector>(radial_distortion_coefs, 0.0);

src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp

@@ -29,8 +29,9 @@
 namespace domain::creators::time_dependent_options {
 KerrSchildFromBoyerLindquist::KerrSchildFromBoyerLindquist() = default;
 KerrSchildFromBoyerLindquist::KerrSchildFromBoyerLindquist(
-    const double mass_in, const std::array<double, 3> spin_in)
-    : mass(mass_in), spin(spin_in) {}
+    const double mass_in, const std::array<double, 3> spin_in,
+    const std::array<double, 3> boost_velocity_in)
+    : mass(mass_in), spin(spin_in), boost_velocity(boost_velocity_in) {}
 
 YlmsFromFile::YlmsFromFile() = default;
 YlmsFromFile::YlmsFromFile(std::string h5_filename_in,
@@ -73,13 +74,13 @@ initial_shape_and_size_funcs(
     if (std::holds_alternative<KerrSchildFromBoyerLindquist>(
             shape_options.initial_values.value())) {
       const ylm::Spherepack ylm{shape_options.l_max, shape_options.l_max};
-      const auto& mass_and_spin = std::get<KerrSchildFromBoyerLindquist>(
+      const auto& ks_from_bl_options = std::get<KerrSchildFromBoyerLindquist>(
           shape_options.initial_values.value());
       const DataVector radial_distortion =
           inner_radius -
           get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-              inner_radius, ylm.theta_phi_points(), mass_and_spin.mass,
-              mass_and_spin.spin));
+              inner_radius, ylm.theta_phi_points(), ks_from_bl_options.mass,
+              ks_from_bl_options.spin, ks_from_bl_options.boost_velocity));
       shape_funcs[0] = ylm.phys_to_spec(radial_distortion);
       // Transform from SPHEREPACK to actual Ylm for size func
       size_funcs[0][0] = shape_funcs[0][0] * sqrt(0.5 * M_PI);

src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp

@@ -23,9 +23,9 @@
 
 namespace domain::creators::time_dependent_options {
 /*!
- * \brief Mass and spin necessary for calculating the \f$ Y_{lm} \f$
- * coefficients of a Kerr horizon of certain Boyer-Lindquist radius for the
- * shape map of the Sphere domain creator.
+ * \brief Mass, spin, and boost velocity necessary for calculating the \f$
+ * Y_{lm} \f$ coefficients of a Kerr horizon of certain Boyer-Lindquist radius
+ * for the shape map of the Sphere domain creator.
  */
 struct KerrSchildFromBoyerLindquist {
   /// \brief The mass of the Kerr black hole.
@@ -39,8 +39,14 @@ struct KerrSchildFromBoyerLindquist {
     static constexpr Options::String help = {
         "The dim'less spin of the Kerr BH."};
   };
+  struct BoostVelocity {
+    using type = std::array<double, 3>;
+    static constexpr Options::String help = {
+        "Boost velocity to transform horizon with in Cartesian Kerr-Schild "
+        "coordinates."};
+  };
 
-  using options = tmpl::list<Mass, Spin>;
+  using options = tmpl::list<Mass, Spin, BoostVelocity>;
 
   static constexpr Options::String help = {
       "Conform to an ellipsoid of constant Boyer-Lindquist radius in "
@@ -49,12 +55,17 @@ struct KerrSchildFromBoyerLindquist {
       "choose an 'InnerRadius' of r_+ = M + sqrt(M^2-a^2)."};
 
   KerrSchildFromBoyerLindquist();
-  KerrSchildFromBoyerLindquist(double mass_in, std::array<double, 3> spin_in);
+  KerrSchildFromBoyerLindquist(double mass_in, std::array<double, 3> spin_in,
+                               std::array<double, 3> boost_velocity_in);
 
   double mass{std::numeric_limits<double>::signaling_NaN()};
   std::array<double, 3> spin{std::numeric_limits<double>::signaling_NaN(),
                              std::numeric_limits<double>::signaling_NaN(),
                              std::numeric_limits<double>::signaling_NaN()};
+  std::array<double, 3> boost_velocity{
+      std::numeric_limits<double>::signaling_NaN(),
+      std::numeric_limits<double>::signaling_NaN(),
+      std::numeric_limits<double>::signaling_NaN()};
 };
 
 /// Label for shape map options

src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.cpp

@@ -17,9 +17,10 @@ template <typename DataType>
 Scalar<DataType> kerr_horizon_radius(
     const std::array<DataType, 2>& theta_phi, const double mass,
     const std::array<double, 3>& dimensionless_spin) {
+  const std::array<double, 3> boost_velocity = {0.0, 0.0, 0.0};
   return kerr_schild_radius_from_boyer_lindquist(
       mass * (1.0 + sqrt(1.0 - square(magnitude(dimensionless_spin)))),
-      theta_phi, mass, dimensionless_spin);
+      theta_phi, mass, dimensionless_spin, boost_velocity);
 }
 
 template Scalar<DataVector> kerr_horizon_radius(
@@ -34,7 +35,8 @@ template <typename DataType>
 Scalar<DataType> kerr_schild_radius_from_boyer_lindquist(
     const double boyer_lindquist_radius,
     const std::array<DataType, 2>& theta_phi, const double mass,
-    const std::array<double, 3>& dimensionless_spin) {
+    const std::array<double, 3>& dimensionless_spin,
+    const std::array<double, 3>& boost_velocity) {
   const double spin_magnitude_squared = square(magnitude(dimensionless_spin));
   const double mass_squared = square(mass);
 
@@ -48,20 +50,37 @@ Scalar<DataType> kerr_schild_radius_from_boyer_lindquist(
                                  dimensionless_spin[1] * sin_theta * sin_phi +
                                  dimensionless_spin[2] * cos_theta;
 
-  return Scalar<DataType>{boyer_lindquist_radius *
-                          sqrt(square(boyer_lindquist_radius) +
-                               mass_squared * spin_magnitude_squared) /
-                          sqrt(square(boyer_lindquist_radius) +
-                               mass_squared * square(spin_dot_unit))};
+  const double boost_velocity_x = boost_velocity[0];
+  const double boost_velocity_y = boost_velocity[1];
+  const double boost_velocity_z = boost_velocity[2];
+  const double lorentz_factor_squared =
+      1.0 / (1.0 - (square(boost_velocity_x) + square(boost_velocity_y) +
+                    square(boost_velocity_z)));
+
+  Scalar<DataType> unboosted_kerr_schild_radius =
+      Scalar<DataType>{boyer_lindquist_radius *
+                       sqrt(square(boyer_lindquist_radius) +
+                            mass_squared * spin_magnitude_squared) /
+                       sqrt(square(boyer_lindquist_radius) +
+                            mass_squared * square(spin_dot_unit))};
+
+  return Scalar<DataType>{sqrt(square(get(unboosted_kerr_schild_radius)) +
+                          lorentz_factor_squared *
+                              square(get(unboosted_kerr_schild_radius)) *
+                              square(cos_phi * sin_theta * boost_velocity_x +
+                                     sin_phi * sin_theta * boost_velocity_y +
+                                     cos_theta * boost_velocity_z))};
 }
 
 template Scalar<DataVector> kerr_schild_radius_from_boyer_lindquist(
     const double boyer_lindquist_radius,
     const std::array<DataVector, 2>& theta_phi, const double mass,
-    const std::array<double, 3>& dimensionless_spin);
+    const std::array<double, 3>& dimensionless_spin,
+    const std::array<double, 3>& boost_velocity);
 
 template Scalar<double> kerr_schild_radius_from_boyer_lindquist(
     const double boyer_lindquist_radius, const std::array<double, 2>& theta_phi,
-    const double mass, const std::array<double, 3>& dimensionless_spin);
+    const double mass, const std::array<double, 3>& dimensionless_spin,
+    const std::array<double, 3>& boost_velocity);
 
 }  // namespace gr::Solutions

src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp

@@ -43,12 +43,25 @@ namespace gr::Solutions {
  * \f]
  * where the angles are encoded in \f$\hat x\f$ and everything else on the
  * right-hand side is constant.
+ *
+ * The Kerr-Schild radius can be Lorentz boosted in an arbitrary direction
+ * \f$\vec{\beta}\f$. We apply a Lorentz transformation on the Cartesian
+ * Kerr-Schild coordinates and find that the boosted Kerr-Schild radius \f$r'\f$
+ * as a function of angles satisfies
+ * \f[
+ * r'^2 = r^2 (1 + \gamma^2 \left( \beta_x\cos\phi\sin\theta +
+ * \beta_y\sin\phi\sin\theta + \beta_z\cos\theta \right)^2)
+ * \f]
+ * where \f$\gamma = \frac{1}{\sqrt{1 - \beta^2}}\f$ is the Lorentz factor, and
+ * the components of \f$\vec{\beta}\f$ are given in Cartesian Kerr-Schild
+ * coordinates.
  */
 template <typename DataType>
 Scalar<DataType> kerr_schild_radius_from_boyer_lindquist(
     const double boyer_lindquist_radius,
     const std::array<DataType, 2>& theta_phi, double mass,
-    const std::array<double, 3>& dimensionless_spin);
+    const std::array<double, 3>& dimensionless_spin,
+    const std::array<double, 3>& boost_velocity);
 /*!
  * \brief The Kerr-Schild radius corresponding to a Kerr horizon.
  *

support/Pipelines/Bbh/InitialData.yaml

@@ -143,13 +143,15 @@ DomainCreator:
         InitialValues:
           Mass: *mass_right
           Spin: *spin_right
+          BoostVelocity: [0.0, 0.0, 0.0]
         SizeInitialValues: [0, 0, 0]
         TransitionEndsAtCube: False
       ShapeMapB:
         LMax: 20
         InitialValues:
           Mass: *mass_left
           Spin: *spin_left
+          BoostVelocity: [0.0, 0.0, 0.0]
         SizeInitialValues: [0, 0, 0]
         TransitionEndsAtCube: False
 

support/Pipelines/Bbh/Inspiral.yaml

@@ -122,6 +122,7 @@ DomainCreator:
             - {{ ExcisionAShapeSpin_x }}
             - {{ ExcisionAShapeSpin_y }}
             - {{ ExcisionAShapeSpin_z }}
+          BoostVelocity: [0.0, 0.0, 0.0]
 {% endif %}
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
@@ -136,6 +137,7 @@ DomainCreator:
             - {{ ExcisionBShapeSpin_x }}
             - {{ ExcisionBShapeSpin_y }}
             - {{ ExcisionBShapeSpin_z }}
+          BoostVelocity: [0.0, 0.0, 0.0]
 {% endif %}
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true

tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp

@@ -846,6 +846,8 @@ void test_kerr_horizon_conforming() {
   const double mass_B = 1.2;
   const std::array<double, 3> spin_A{{0.0, 0.0, 0.9}};
   const std::array<double, 3> spin_B{{0.0, 0.2, 0.4}};
+  const std::array<double, 3> boost_velocity_A{{0.0, 0.0, 0.0}};
+  const std::array<double, 3> boost_velocity_B{{0.0, 0.0, 0.0}};
   const double r_plus_A = mass_A * (1. + sqrt(1. - dot(spin_A, spin_A)));
   const double r_plus_B = mass_B * (1. + sqrt(1. - dot(spin_B, spin_B)));
   const double inner_radius_A = r_plus_A;
@@ -874,11 +876,11 @@ void test_kerr_horizon_conforming() {
           std::nullopt,
           {{32_st,
             domain::creators::time_dependent_options::
-                KerrSchildFromBoyerLindquist{mass_A, spin_A},
+                KerrSchildFromBoyerLindquist{mass_A, spin_A, boost_velocity_A},
             std::nullopt}},
           {{32_st,
             domain::creators::time_dependent_options::
-                KerrSchildFromBoyerLindquist{mass_B, spin_B},
+                KerrSchildFromBoyerLindquist{mass_B, spin_B, boost_velocity_B},
             std::nullopt}}}};
   const auto domain = domain_creator.create_domain();
   const auto functions_of_time = domain_creator.functions_of_time();
@@ -894,10 +896,10 @@ void test_kerr_horizon_conforming() {
            make_not_null(&gen), make_not_null(&dist_phi), num_points)}};
   const auto radius_A =
       get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-          inner_radius_A, theta_phi, mass_A, spin_A));
+          inner_radius_A, theta_phi, mass_A, spin_A, boost_velocity_A));
   const auto radius_B =
       get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-          inner_radius_B, theta_phi, mass_B, spin_B));
+          inner_radius_B, theta_phi, mass_B, spin_B, boost_velocity_B));
   tnsr::I<DataVector, 3> x_A{};
   tnsr::I<DataVector, 3> x_B{};
   get<0>(x_A) =

tests/Unit/Domain/Creators/Test_Sphere.cpp

@@ -792,12 +792,13 @@ void test_shape_distortion() {
   const double time = 0.7;
   const double mass = 0.8;
   const std::array<double, 3> spin{{0.0, 0.0, 0.9}};
+  const std::array<double, 3> boost_velocity{{0.0, 0.0, 0.0}};
   const double r_plus = mass * (1. + sqrt(1. - dot(spin, spin)));
   const double inner_radius = r_plus;
 
   const DataVector radius =
       get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-          inner_radius, theta_phi, mass, spin));
+          inner_radius, theta_phi, mass, spin, boost_velocity));
   // Set up coordinates on an ellipsoid of constant Boyer-Lindquist radius
   tnsr::I<DataVector, 3> x{};
   get<0>(x) = radius * sin(get<0>(theta_phi)) * cos(get<1>(theta_phi));
@@ -822,7 +823,7 @@ void test_shape_distortion() {
     time_dependent_options = domain::creators::sphere::TimeDependentMapOptions{
         time,
         {{l_max, domain::creators::time_dependent_options::
-                     KerrSchildFromBoyerLindquist{mass, spin}}},
+                     KerrSchildFromBoyerLindquist{mass, spin, boost_velocity}}},
         std::nullopt,
         std::nullopt,
         std::nullopt};

tests/Unit/Domain/Creators/TimeDependence/Test_Shape.cpp

@@ -117,7 +117,8 @@ void test_r_theta_phi(const tnsr::I<double, 3, SourceFrame>& input,
   CHECK(input_theta_phi[1] == approx(output_centered_spherical[2]));
   const double kerr_schild_radius =
       gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-          inner_radius, input_theta_phi, mass, spin)
+          inner_radius, input_theta_phi, mass, spin,
+          std::array<double, 3>{0.0, 0.0, 0.0})
           .get();
   const double radius = magnitude(input_centered).get();
   const double transition_factor =

tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp

@@ -40,9 +40,12 @@ void test_kerr_schild_boyer_lindquist() {
   const auto kerr_schild_boyer_lindquist = TestHelpers::test_creation<
       domain::creators::time_dependent_options::KerrSchildFromBoyerLindquist>(
       "Mass: 1.7\n"
-      "Spin: [0.45, 0.12, 0.34]");
+      "Spin: [0.45, 0.12, 0.34]\n"
+      "BoostVelocity: [0.0, 0.0, 0.0]");
   CHECK(kerr_schild_boyer_lindquist.mass == 1.7);
   CHECK(kerr_schild_boyer_lindquist.spin == std::array{0.45, 0.12, 0.34});
+  CHECK(kerr_schild_boyer_lindquist.boost_velocity ==
+        std::array{0.0, 0.0, 0.0});
 }
 
 void test_ylms_from_file() {
@@ -97,6 +100,7 @@ void test_shape_map_options() {
         "InitialValues:\n"
         "  Mass: 1.7\n"
         "  Spin: [0.45, 0.12, 0.34]\n"
+        "  BoostVelocity: [0.0, 0.0, 0.0]\n"
         "SizeInitialValues: Auto\n"
         "TransitionEndsAtCube: True");
     CHECK(shape_map_options.name() == "ShapeMapB");
@@ -174,6 +178,7 @@ void test_funcs(const gsl::not_null<Generator*> generator) {
         "InitialValues:\n"
         "  Mass: 1.0\n"
         "  Spin: [0.0, 0.0, 0.0]\n"
+        "  BoostVelocity: [0.0, 0.0, 0.0]\n"
         "SizeInitialValues: [0.5, 1.0, 2.4]");
 
     const auto [shape_funcs, size_funcs] =
@@ -195,6 +200,7 @@ void test_funcs(const gsl::not_null<Generator*> generator) {
         "InitialValues:\n"
         "  Mass: 1.0\n"
         "  Spin: [0.0, 0.0, 0.0]\n"
+        "  BoostVelocity: [0.0, 0.0, 0.0]\n"
         "SizeInitialValues: Auto");
 
     const auto [shape_funcs, size_funcs] =

tests/Unit/Domain/Creators/TimeDependentOptions/Test_Sphere.cpp

@@ -38,10 +38,12 @@ std::string create_option_string(const std::optional<bool> use_non_zero_shape) {
                 "  LMax: 8\n"
                 "  SizeInitialValues: Auto\n"
                 "  InitialValues:" +
-                    (use_non_zero_shape.value() ? "\n"
-                                                  "    Mass: 1.0\n"
-                                                  "    Spin: [0.0, 0.0, 0.0]\n"s
-                                                : " Spherical\n"s)
+                    (use_non_zero_shape.value()
+                         ? "\n"
+                           "    Mass: 1.0\n"
+                           "    Spin: [0.0, 0.0, 0.0]\n"
+                           "    BoostVelocity: [0.0, 0.0, 0.0]\n"s
+                         : " Spherical\n"s)
               : " None\n") +
          "RotationMap: None\n"
          "ExpansionMap: None\n"

tests/Unit/Domain/Test_StrahlkorperTransformations.cpp

@@ -91,7 +91,8 @@ void test_strahlkorper_in_different_frame() {
     const DataVector new_radius =
         get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
             2.0, ylm.theta_phi_points(), 1.0,
-            std::array<double, 3>{{0.1, 0.2, 0.3}}));
+            std::array<double, 3>{{0.1, 0.2, 0.3}},
+            std::array<double, 3>{{0.0, 0.0, 0.0}}));
     strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(
         l_max, l_max, new_radius, strahlkorper_src_center));
   } else if constexpr (std::is_same_v<SrcFrame, ::Frame::Inertial>) {

tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_KerrHorizon.cpp

@@ -11,6 +11,7 @@
 #include "Framework/CheckWithRandomValues.hpp"
 #include "Framework/SetupLocalPythonEnvironment.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
+#include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp"
 #include "Utilities/ConstantExpressions.hpp"
 
 namespace gr::Solutions {
@@ -104,9 +105,35 @@
     const DataVector phi{M_PI_2 * 0.4, M_PI * 0.55, M_PI_4, M_PI_4 * 1.2, M_PI};
     const std::array<DataVector, 2> theta_phi{
         {std::move(theta), std::move(phi)}};
-    CHECK_ITERABLE_APPROX(kerr_schild_radius_from_boyer_lindquist(
-                              r_plus, theta_phi, mass, dimless_spin),
-                          kerr_horizon_radius(theta_phi, mass, dimless_spin));
+    const std::array<double, 3> no_boost_velocity = {0.0, 0.0, 0.0};
+    CHECK_ITERABLE_APPROX(
+        kerr_schild_radius_from_boyer_lindquist(
+            r_plus, theta_phi, mass, dimless_spin, no_boost_velocity),
+        kerr_horizon_radius(theta_phi, mass, dimless_spin));
+
+    const Scalar<DataVector> unboosted_radius =
+        kerr_horizon_radius(theta_phi, mass, dimless_spin);
+    tnsr::I<DataVector, 3> unboosted_coords = {};
+    get<0>(unboosted_coords) = get(unboosted_radius) * cos(phi) * sin(theta);
+    get<1>(unboosted_coords) = get(unboosted_radius) * sin(phi) * sin(theta);
+    get<2>(unboosted_coords) = get(unboosted_radius) * cos(theta);
+    const std::array<double, 3> boost_velocity = {0.1, -0.05, 0.07};
+    tnsr::I<DataVector, 3> boosted_coords =
+        make_with_value<tnsr::I<DataVector, 3>>(
+            get<0>(unboosted_coords),
+            std::numeric_limits<double>::signaling_NaN());
+    sr::lorentz_boost(make_not_null(&boosted_coords), unboosted_coords, 0.0,
+                      boost_velocity);
+    Scalar<DataVector> boosted_radius = Scalar<DataVector>{
+        sqrt(square(get<0>(boosted_coords)) + square(get<1>(boosted_coords)) +
+             square(get<2>(boosted_coords)))};
+    // Test for boosted Kerr horizon.
+    // Boosted horizon should be equivalent to applying Lorentz transformation
+    // on unboosted coordinates.
+    CHECK_ITERABLE_APPROX(
+        kerr_schild_radius_from_boyer_lindquist(r_plus, theta_phi, mass,
+                                                dimless_spin, boost_velocity),
+        boosted_radius);
   }
 }
 }  // namespace gr::Solutions