Merge pull request sxs-collaboration#5610 from guilara/PRs-KerrBoost

src/PointwiseFunctions/AnalyticData/GrMhd/SpecInitialData.hpp

@@ -61,6 +61,9 @@ class SpecInitialData : public evolution::initial_data::InitialData,
                  gr::Tags::Lapse<DataType>, gr::Tags::Shift<DataType, 3>>,
       hydro::grmhd_tags<DataType>>;
 
+  static std::string name() {
+    return "SpecInitialData" + std::to_string(ThermodynamicDim) + "dEos";
+  }
   struct DataDirectory {
     using type = std::string;
     static constexpr Options::String help = {

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/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml

@@ -415,7 +415,7 @@ ControlSystems:
     ControlError:
 
 InitialData:
-  # SpecInitialData:
+  # SpecInitialData3dEos:
   #   DataDirectory: "/path/to/EvID"
   #   EquationOfState: *eos
   #   DensityCutoff: *AtmosphereDensityCutoff

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