Add boost to KerrSchild

Add more Lorentz boost functions

Add instantiation of time deriv of spatial metric

Add boosted Kerr initial data

Expand doc

Fix initialization of lorentz boost matrix

Remove unused reference

Add test of boosted spacetime metric

Throw error for unsupported tag
>
>
Co-authored-by: Nils Vu <nilsvu@caltech.edu>

src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/CMakeLists.txt

@@ -44,6 +44,7 @@ target_link_libraries(
   GeneralRelativity
   Interpolation
   Serialization
+  SpecialRelativity
   )
 
 add_subdirectory(Python)

src/PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpatialMetric.cpp

@@ -81,7 +81,7 @@ tnsr::ii<DataType, SpatialDim, Frame> time_derivative_of_spatial_metric(
           extrinsic_curvature);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector),
-                        (Frame::Grid, Frame::Inertial))
+                        (Frame::Grid, Frame::Distorted, Frame::Inertial))
 
 #undef DIM
 #undef DTYPE

src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.cpp

@@ -3,6 +3,7 @@
 
 #include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp"
 
+#include <array>
 #include <cmath>
 #include <cstddef>
 #include <limits>
@@ -58,6 +59,75 @@ void lorentz_boost_matrix(
     (*boost_matrix).get(i + 1, i + 1) += 1.0;
   }
 }
+
+template <size_t SpatialDim>
+tnsr::Ab<double, SpatialDim, Frame::NoFrame> lorentz_boost_matrix(
+    const std::array<double, SpatialDim>& velocity) {
+  tnsr::I<double, SpatialDim, Frame::NoFrame> velocity_as_tensor(velocity);
+  return lorentz_boost_matrix(velocity_as_tensor);
+}
+
+template <typename DataType, size_t SpatialDim, typename Frame>
+void lorentz_boost(
+    const gsl::not_null<tnsr::I<DataType, SpatialDim, Frame>*> result,
+    const tnsr::I<DataType, SpatialDim, Frame>& vector,
+    const double vector_component_0,
+    const std::array<double, SpatialDim>& velocity) {
+  if (velocity == make_array<SpatialDim>(0.)) {
+    *result = vector;
+    return;
+  }
+  // Inverse matrix with respect to the boost applied to one forms
+  const auto boost_matrix = lorentz_boost_matrix(-velocity);
+  for (size_t i = 0; i < SpatialDim; ++i) {
+    result->get(i) = boost_matrix.get(i + 1, 0) * vector_component_0;
+    for (size_t j = 0; j < SpatialDim; ++j) {
+      result->get(i) += boost_matrix.get(i + 1, j + 1) * vector.get(j);
+    }
+  }
+}
+
+template <typename DataType, size_t SpatialDim, typename Frame>
+void lorentz_boost(
+    const gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> result,
+    const tnsr::a<DataType, SpatialDim, Frame>& one_form,
+    const std::array<double, SpatialDim>& velocity) {
+  if (velocity == make_array<SpatialDim>(0.)) {
+    *result = one_form;
+    return;
+  }
+  const auto boost_matrix = lorentz_boost_matrix(velocity);
+  for (size_t i = 0; i < SpatialDim + 1; ++i) {
+    result->get(i) = 0.;
+    for (size_t j = 0; j < SpatialDim + 1; ++j) {
+      result->get(i) += boost_matrix.get(i, j) * one_form.get(j);
+    }
+  }
+}
+
+template <typename DataType, size_t SpatialDim, typename Frame>
+void lorentz_boost(
+    const gsl::not_null<tnsr::ab<DataType, SpatialDim, Frame>*> result,
+    const tnsr::ab<DataType, SpatialDim, Frame>& tensor,
+    const std::array<double, SpatialDim>& velocity) {
+  if (velocity == make_array<SpatialDim>(0.)) {
+    *result = tensor;
+    return;
+  }
+  const auto boost_matrix = lorentz_boost_matrix(velocity);
+  for (size_t i = 0; i < SpatialDim + 1; ++i) {
+    for (size_t k = 0; k < SpatialDim + 1; ++k) {
+      result->get(i, k) = 0.;
+      for (size_t j = 0; j < SpatialDim + 1; ++j) {
+        for (size_t l = 0; l < SpatialDim + 1; ++l) {
+          result->get(i, k) += boost_matrix.get(i, j) * boost_matrix.get(k, l) *
+                               tensor.get(j, l);
+        }
+      }
+    }
+  }
+}
+
 }  // namespace sr
 
 // Explicit Instantiations
@@ -70,10 +140,40 @@ void lorentz_boost_matrix(
   template void sr::lorentz_boost_matrix(                          \
       gsl::not_null<tnsr::Ab<double, DIM(data), Frame::NoFrame>*>  \
           boost_matrix,                                            \
-      const tnsr::I<double, DIM(data), Frame::NoFrame>& velocity);
+      const tnsr::I<double, DIM(data), Frame::NoFrame>& velocity); \
+  template tnsr::Ab<double, DIM(data), Frame::NoFrame>             \
+  sr::lorentz_boost_matrix(const std::array<double, DIM(data)>& velocity);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
 
+#undef DIM
+#undef INSTANTIATE
+
+#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data)
+
+#define INSTANTIATE(_, data)                                              \
+  template void sr::lorentz_boost(                                        \
+      const gsl::not_null<tnsr::I<DTYPE(data), DIM(data), FRAME(data)>*>  \
+          result,                                                         \
+      const tnsr::I<DTYPE(data), DIM(data), FRAME(data)>& one_form,       \
+      const double one_form_component_0,                                  \
+      const std::array<double, DIM(data)>& velocity);                     \
+  template void sr::lorentz_boost(                                        \
+      const gsl::not_null<tnsr::a<DTYPE(data), DIM(data), FRAME(data)>*>  \
+          result,                                                         \
+      const tnsr::a<DTYPE(data), DIM(data), FRAME(data)>& one_form,       \
+      const std::array<double, DIM(data)>& velocity);                     \
+  template void sr::lorentz_boost(                                        \
+      const gsl::not_null<tnsr::ab<DTYPE(data), DIM(data), FRAME(data)>*> \
+          result,                                                         \
+      const tnsr::ab<DTYPE(data), DIM(data), FRAME(data)>& tensor,        \
+      const std::array<double, DIM(data)>& velocity);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector),
+                        (Frame::Grid, Frame::Distorted, Frame::Inertial))
+
 #undef DIM
 #undef DTYPE
 #undef FRAME

src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp

@@ -3,9 +3,12 @@
 
 #pragma once
 
+#include <array>
 #include <cstddef>
 
-#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/MakeArray.hpp"
 
 /// \cond
 namespace gsl {
@@ -53,5 +56,43 @@ template <size_t SpatialDim>
 void lorentz_boost_matrix(
     gsl::not_null<tnsr::Ab<double, SpatialDim, Frame::NoFrame>*> boost_matrix,
     const tnsr::I<double, SpatialDim, Frame::NoFrame>& velocity);
+
+template <size_t SpatialDim>
+tnsr::Ab<double, SpatialDim, Frame::NoFrame> lorentz_boost_matrix(
+    const std::array<double, SpatialDim>& velocity);
+/// @}
+
+/// @{
+/*!
+ * \ingroup SpecialRelativityGroup
+ * \brief Apply a Lorentz boost to the spatial part of a vector.
+ * \details This requires passing the 0th component of the vector as an
+ * additional argument.
+ */
+template <typename DataType, size_t SpatialDim, typename Frame>
+void lorentz_boost(gsl::not_null<tnsr::I<DataType, SpatialDim, Frame>*> result,
+                   const tnsr::I<DataType, SpatialDim, Frame>& vector,
+                   double vector_component_0,
+                   const std::array<double, SpatialDim>& velocity);
+
+/*!
+ * \brief Apply a Lorentz boost to a one form.
+ */
+template <typename DataType, size_t SpatialDim, typename Frame>
+void lorentz_boost(
+    gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> result,
+    const tnsr::a<DataType, SpatialDim, Frame>& one_form,
+    const std::array<double, SpatialDim>& velocity);
+
+/*!
+ * \brief Apply a Lorentz boost to each component of a rank-2 tensor with
+ * lower or covariant indices.
+ * \note In the future we might want to write a single function capable to boost
+ * a tensor of arbitrary rank.
+ */
+template <typename DataType, size_t SpatialDim, typename Frame>
+void lorentz_boost(gsl::not_null<tnsr::ab<DataType, SpatialDim, Frame>*> result,
+                   const tnsr::ab<DataType, SpatialDim, Frame>& tensor,
+                   const std::array<double, SpatialDim>& velocity);
 /// @}
 }  // namespace sr

support/Pipelines/Bbh/InitialData.yaml

@@ -34,6 +34,7 @@ Background: &background
           - {{ DimensionlessSpinLeft_y }}
           - {{ DimensionlessSpinLeft_z }}
         Center: [0., 0., 0.]
+        Velocity: [0., 0., 0.]
     ObjectRight: &kerr_right
       KerrSchild:
         Mass: &mass_right {{ MassRight }}
@@ -42,6 +43,7 @@ Background: &background
           - {{ DimensionlessSpinRight_y }}
           - {{ DimensionlessSpinRight_z }}
         Center: [0., 0., 0.]
+        Velocity: [0., 0., 0.]
     AngularVelocity: {{ OrbitalAngularVelocity }}
     Expansion: {{ RadialExpansionVelocity }}
     LinearVelocity: [0., 0., 0.]

tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml

@@ -26,6 +26,7 @@ BackgroundSpacetime:
     Mass: 1.
     Center: [0., 0., 0.]
     Spin: [0., 0., 0.]
+    Velocity: [0., 0., 0.]
 
 ResourceInfo:
   AvoidGlobalProc0: false

tests/InputFiles/GeneralizedHarmonic/KerrSchild.yaml

@@ -74,6 +74,7 @@ InitialData: &InitialData
     Mass: 1.0
     Spin: [0.0, 0.0, 0.0]
     Center: &Center [0.0, 0.0, 0.0]
+    Velocity: [0.0, 0.0, 0.0]
 
 DomainCreator:
   Sphere:

tests/InputFiles/Xcts/BinaryBlackHole.yaml

@@ -84,11 +84,13 @@ Background: &background
         Mass: 0.4229
         Spin: [0., 0., 0.]
         Center: [0., 0., 0.]
+        Velocity: [0., 0., 0.]
     ObjectRight: &kerr_right
       KerrSchild:
         Mass: 0.4229
         Spin: [0., 0., 0.]
         Center: [0., 0., 0.]
+        Velocity: [0., 0., 0.]
     AngularVelocity: 0.0144
     Expansion: 0.
     LinearVelocity: [0., 0., 0.]

tests/InputFiles/Xcts/KerrSchild.yaml

@@ -71,6 +71,7 @@ Background: &solution
     Mass: &KerrMass 1.
     Spin: &KerrSpin [0., 0., 0.6]
     Center: &KerrCenter [0., 0., 0.]
+    Velocity: [0., 0., 0.]
 
 InitialGuess: Flatness
 

tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_ApparentHorizon.cpp

@@ -411,7 +411,8 @@ void test_consistency_with_kerr(const bool compute_expansion) {
   std::array<double, 3> rotation =
       -0.5 * dimensionless_spin / horizon_kerrschild_radius;
   CAPTURE(rotation);
-  const KerrSchild solution{mass, dimensionless_spin, {{0., 0., 0.}}};
+  const KerrSchild solution{
+      mass, dimensionless_spin, {{0., 0., 0.}}, {{0., 0., 0.}}};
   const ApparentHorizon<Xcts::Geometry::Curved> kerr_horizon{
       center, rotation, std::make_unique<KerrSchild>(solution),
       // Check with and without the negative-expansion condition. Either the
@@ -579,8 +580,8 @@ SPECTRE_TEST_CASE("Unit.Xcts.BoundaryConditions.ApparentHorizon",
                 "  Lapse: Auto\n"
                 "  NegativeExpansion: None\n");
   test_creation({{1., 2., 3.}}, {{0.1, 0.2, 0.3}},
-                {{2.3, {{0.4, 0.5, 0.6}}, {{0., 0., 0.}}}},
-                {{3.4, {{0.3, 0.2, 0.1}}, {{0., 0., 0.}}}},
+                {{2.3, {{0.4, 0.5, 0.6}}, {{0., 0., 0.}}, {{0., 0., 0.}}}},
+                {{3.4, {{0.3, 0.2, 0.1}}, {{0., 0., 0.}}, {{0., 0., 0.}}}},
                 "ApparentHorizon:\n"
                 "  Center: [1., 2., 3.]\n"
                 "  Rotation: [0.1, 0.2, 0.3]\n"
@@ -589,11 +590,13 @@ SPECTRE_TEST_CASE("Unit.Xcts.BoundaryConditions.ApparentHorizon",
                 "      Mass: 2.3\n"
                 "      Spin: [0.4, 0.5, 0.6]\n"
                 "      Center: [0., 0., 0.]\n"
+                "      Velocity: [0., 0., 0.]\n"
                 "  NegativeExpansion:\n"
                 "    KerrSchild:\n"
                 "      Mass: 3.4\n"
                 "      Spin: [0.3, 0.2, 0.1]\n"
-                "      Center: [0., 0., 0.]\n");
+                "      Center: [0., 0., 0.]\n"
+                "      Velocity: [0., 0., 0.]\n");
   test_with_random_values();
   test_consistency_with_kerr(false);
   test_consistency_with_kerr(true);

tests/Unit/Evolution/Systems/CurvedScalarWave/Test_BackgroundSpacetime.cpp

@@ -37,7 +37,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.BackgroundSpacetime",
   const auto ks_background =
       gr::Solutions::KerrSchild(0.5, {{0.1, 0.2, 0.3}}, {{1.0, 3.0, 2.0}});
   const std::string& ks_option_string =
-      "Mass: 0.5\nSpin: [0.1,0.2,0.3]\nCenter: [1.0,3.0,2.0]";
+      "Mass: 0.5\nSpin: [0.1,0.2,0.3]\nCenter: [1.0,3.0,2.0]\nVelocity: "
+      "[0.0,0.0,0.0]";
   test_option_parsing(ks_background, ks_option_string);
 
   const auto spherical_ks_background = gr::Solutions::SphericalKerrSchild(

tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp

@@ -319,7 +319,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.Gh.NumericInitialData",
       "GeneralizedHarmonic(KerrSchild):\n"
       "  Mass: 1.\n"
       "  Spin: [0, 0, 0]\n"
-      "  Center: [0, 0, 0]",
+      "  Center: [0, 0, 0]\n"
+      "  Velocity: [0, 0, 0]",
       false);
 }
 

tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_AnalyticChristoffel.cpp

@@ -122,7 +122,8 @@ void test_ks(const Mesh<3>& mesh) {
           "    GeneralizedHarmonic(KerrSchild):\n"
           "      Mass: 1.2\n"
           "      Spin: [0.1, 0.2, 0.3]\n"
-          "      Center: [-0.1, -0.2, -0.4]\n")
+          "      Center: [-0.1, -0.2, -0.4]\n"
+          "      Velocity: [0.0, 0.0, 0.0]\n")
           ->get_clone());
   CHECK_FALSE(gauge_condition->is_harmonic());
 

tests/Unit/Helpers/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/VerifyGrSolution.hpp

@@ -426,5 +426,55 @@ void verify_consistency(const Solution& solution, const double time,
                                               derivative_delta),
       derivative_approx);
 }
+
+/// Check the consistency of quantities dependent on other metric quantities
+/// returned by a solution.  This includes checking pointwise relations such as
+/// consistency of the metric and inverse and comparing returned and
+/// numerical derivatives.
+template <typename Solution, typename Frame>
+void verify_spatial_consistency(const Solution& solution, const double time,
+                                const tnsr::I<double, 3, Frame>& position,
+                                const double derivative_delta,
+                                const double derivative_tolerance) {
+  using Lapse = gr::Tags::Lapse<double>;
+  using Shift = gr::Tags::Shift<double, 3, Frame>;
+  using SpatialMetric = gr::Tags::SpatialMetric<double, 3, Frame>;
+  using SqrtDetSpatialMetric = gr::Tags::SqrtDetSpatialMetric<double>;
+  using InverseSpatialMetric = gr::Tags::InverseSpatialMetric<double, 3, Frame>;
+  using ExtrinsicCurvature = gr::Tags::ExtrinsicCurvature<double, 3, Frame>;
+  using tags =
+      tmpl::list<SpatialMetric, SqrtDetSpatialMetric, InverseSpatialMetric,
+                 ExtrinsicCurvature, Lapse, Shift, detail::deriv<Shift, Frame>,
+                 Tags::dt<SpatialMetric>, detail::deriv<SpatialMetric, Frame>,
+                 Tags::dt<Lapse>, Tags::dt<Shift>, detail::deriv<Lapse, Frame>>;
+
+  auto derivative_approx = approx.epsilon(derivative_tolerance);
+
+  const auto vars = solution.variables(position, time, tags{});
+
+  const auto numerical_metric_det_and_inverse =
+      determinant_and_inverse(get<SpatialMetric>(vars));
+  CHECK_ITERABLE_APPROX(get(get<SqrtDetSpatialMetric>(vars)),
+                        sqrt(get(numerical_metric_det_and_inverse.first)));
+  CHECK_ITERABLE_APPROX(get<InverseSpatialMetric>(vars),
+                        numerical_metric_det_and_inverse.second);
+
+  CHECK_ITERABLE_CUSTOM_APPROX(
+      SINGLE_ARG(get<detail::deriv<Lapse, Frame>>(vars)),
+      detail::space_derivative<Lapse>(solution, position, time,
+                                      derivative_delta),
+      derivative_approx);
+  CHECK_ITERABLE_CUSTOM_APPROX(
+      SINGLE_ARG(get<detail::deriv<Shift, Frame>>(vars)),
+      detail::space_derivative<Shift>(solution, position, time,
+                                      derivative_delta),
+      derivative_approx);
+  CHECK_ITERABLE_CUSTOM_APPROX(
+      SINGLE_ARG(get<detail::deriv<SpatialMetric, Frame>>(vars)),
+      detail::space_derivative<SpatialMetric>(solution, position, time,
+                                              derivative_delta),
+      derivative_approx);
+}
+
 }  // namespace VerifyGrSolution
 }  // namespace TestHelpers