Merge pull request sxs-collaboration#5681 from nikwit/geodesic-accele…

…ration

Add geodesic coordinate acceleration function

src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt

@@ -30,6 +30,8 @@ target_link_libraries(
   PUBLIC
   DataStructures
   Domain
+  GeneralRelativity
+  GeneralRelativitySolutions
   Options
   Parallel
   Utilities

src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.cpp

@@ -18,6 +18,9 @@
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/PunctureField.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
+#include "PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/SetNumberOfGridPoints.hpp"
 
@@ -139,7 +142,23 @@ void ParticlePositionVelocityCompute<Dim>::function(
   (*position_and_velocity)[1] = std::move(std::get<3>(mapped_tuple));
 }
 
+template <size_t Dim>
+void GeodesicAccelerationCompute<Dim>::function(
+    gsl::not_null<tnsr::I<double, Dim, Frame::Inertial>*> acceleration,
+    const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>&
+        position_velocity,
+    const gr::Solutions::KerrSchild& background_spacetime) {
+  const auto christoffel = get<
+      gr::Tags::SpacetimeChristoffelSecondKind<double, Dim, Frame::Inertial>>(
+      background_spacetime.variables(
+          position_velocity.at(0), 0.,
+          tmpl::list<gr::Tags::SpacetimeChristoffelSecondKind<
+              double, Dim, Frame::Inertial>>{}));
+  gr::geodesic_acceleration(acceleration, position_velocity.at(1), christoffel);
+}
+
 template struct ParticlePositionVelocityCompute<3>;
+template struct GeodesicAccelerationCompute<3>;
 template struct PunctureFieldCompute<3>;
 
 template struct FaceCoordinatesCompute<3, Frame::Grid, true>;

src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp

@@ -217,6 +217,31 @@ struct ParticlePositionVelocityCompute : ParticlePositionVelocity<Dim>,
 };
 /// @}
 
+/// @{
+/*!
+ * \brief Computes the coordinate geodesic acceleration of the particle in the
+ * inertial frame in Kerr-Schild coordinates.
+ */
+template <size_t Dim>
+struct GeodesicAcceleration : db::SimpleTag {
+  using type = tnsr::I<double, Dim, Frame::Inertial>;
+};
+
+template <size_t Dim>
+struct GeodesicAccelerationCompute : GeodesicAcceleration<Dim>, db::ComputeTag {
+  using base = GeodesicAcceleration<Dim>;
+  using return_type = tnsr::I<double, Dim, Frame::Inertial>;
+  using argument_tags = tmpl::list<
+      ParticlePositionVelocity<Dim>,
+      CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>>;
+  static void function(
+      gsl::not_null<tnsr::I<double, Dim, Frame::Inertial>*> acceleration,
+      const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>&
+          position_velocity,
+      const gr::Solutions::KerrSchild& background_spacetime);
+};
+/// @}
+
 /// @{
 /*!
  * \brief An optional that holds the coordinates of an element face abutting the

src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt

@@ -12,6 +12,7 @@ spectre_target_sources(
   DerivativeSpatialMetric.cpp
   DerivativesOfSpacetimeMetric.cpp
   ExtrinsicCurvature.cpp
+  GeodesicAcceleration.cpp
   InterfaceNullNormal.cpp
   InverseSpacetimeMetric.cpp
   KerrSchildCoords.cpp
@@ -43,6 +44,7 @@ spectre_target_headers(
   DerivativesOfSpacetimeMetric.hpp
   DetAndInverseSpatialMetric.hpp
   ExtrinsicCurvature.hpp
+  GeodesicAcceleration.hpp
   InterfaceNullNormal.hpp
   InverseSpacetimeMetric.hpp
   KerrSchildCoords.hpp

src/PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.cpp

@@ -0,0 +1,61 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace gr {
+template <typename DataType, size_t Dim>
+void geodesic_acceleration(
+    gsl::not_null<tnsr::I<DataType, Dim>*> acceleration,
+    const tnsr::I<DataType, Dim>& velocity,
+    const tnsr::Abb<DataType, Dim>& christoffel_second_kind) {
+  for (size_t i = 0; i < Dim; ++i) {
+    acceleration->get(i) =
+        velocity.get(i) * christoffel_second_kind.get(0, 0, 0) -
+        christoffel_second_kind.get(i + 1, 0, 0);
+    for (size_t j = 0; j < Dim; ++j) {
+      acceleration->get(i) +=
+          2. * velocity.get(j) *
+          (velocity.get(i) * christoffel_second_kind.get(0, j + 1, 0) -
+           christoffel_second_kind.get(i + 1, j + 1, 0));
+      for (size_t k = 0; k < Dim; ++k) {
+        acceleration->get(i) +=
+            velocity.get(j) * velocity.get(k) *
+            (velocity.get(i) * christoffel_second_kind.get(0, j + 1, k + 1) -
+             christoffel_second_kind.get(i + 1, j + 1, k + 1));
+      }
+    }
+  }
+}
+
+template <typename DataType, size_t Dim>
+tnsr::I<DataType, Dim> geodesic_acceleration(
+    const tnsr::I<DataType, Dim>& velocity,
+    const tnsr::Abb<DataType, Dim>& christoffel_second_kind) {
+  tnsr::I<DataType, Dim> acceleration{get_size(get<0>(velocity))};
+  geodesic_acceleration(make_not_null(&acceleration), velocity,
+                        christoffel_second_kind);
+  return acceleration;
+}
+
+}  // namespace gr
+template void gr::geodesic_acceleration(
+    const gsl::not_null<tnsr::I<double, 3>*> acceleration,
+    const tnsr::I<double, 3>& velocity,
+    const tnsr::Abb<double, 3>& christoffel_second_kind);
+template tnsr::I<double, 3> gr::geodesic_acceleration(
+    const tnsr::I<double, 3>& velocity,
+    const tnsr::Abb<double, 3>& christoffel_second_kind);
+template void gr::geodesic_acceleration(
+    const gsl::not_null<tnsr::I<DataVector, 3>*> acceleration,
+    const tnsr::I<DataVector, 3>& velocity,
+    const tnsr::Abb<DataVector, 3>& christoffel_second_kind);
+template tnsr::I<DataVector, 3> gr::geodesic_acceleration(
+    const tnsr::I<DataVector, 3>& velocity,
+    const tnsr::Abb<DataVector, 3>& christoffel_second_kind);

src/PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp

@@ -0,0 +1,38 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace gr {
+/// @{
+
+/*!
+ * \brief Computes the coordinate geodesic acceleration in the inertial frame.
+ *
+ * \details The geodesic acceleration in coordinate form is given by
+ * \f{equation}
+ *  \frac{d^2 x^i}{d t^2} = (v^i \Gamma^0_{00} - \Gamma^i_{00} )  + 2 v^j (v^i
+ * \Gamma^0_{j0} - \Gamma^i_{j0} ) + v^j v^k (v^i \Gamma^0_{jk} - \Gamma^i_{jk}
+ * ), \f}
+ * where \f$v^i\f$ is the coordinate velocity, \f$\Gamma^\mu_{\nu \rho}\f$ are
+ * the spacetime Christoffel symbols of the second kind, and all latin indices
+ * are spatial.
+ */
+template <typename DataType, size_t Dim>
+void geodesic_acceleration(
+    gsl::not_null<tnsr::I<DataType, Dim>*> acceleration,
+    const tnsr::I<DataType, Dim>& velocity,
+    const tnsr::Abb<DataType, Dim>& christoffel_second_kind);
+
+template <typename DataType, size_t Dim>
+tnsr::I<DataType, Dim> geodesic_acceleration(
+    const tnsr::I<DataType, Dim>& velocity,
+    const tnsr::Abb<DataType, Dim>& christoffel_second_kind);
+/// @}
+
+}  // namespace gr

tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt

@@ -41,6 +41,7 @@ target_link_libraries(
   DomainBoundaryConditions
   DomainBoundaryConditionsHelpers
   DomainCreators
+  GeneralRelativity
   GeneralRelativityHelpers
   GeneralRelativitySolutions
   MathFunctions

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp

@@ -38,6 +38,8 @@
 #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
+#include "PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 #include "Time/Tags/Time.hpp"
 #include "Utilities/CartesianProduct.hpp"
 #include "Utilities/TMPL.hpp"
@@ -342,6 +344,35 @@ void test_particle_position_velocity_compute() {
   }
 }
 
+void test_geodesic_acceleration_compute() {
+  static constexpr size_t Dim = 3;
+  MAKE_GENERATOR(gen);
+  std::uniform_real_distribution<> dist(1., 100.);
+  const auto random_position = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const auto random_velocity = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const gr::Solutions::KerrSchild kerr_schild(0.1, make_array(0.5, 0.1, 0.2),
+                                              make_array(0.3, 0.1, 0.2));
+  auto box =
+      db::create<db::AddSimpleTags<Tags::ParticlePositionVelocity<Dim>,
+                                   CurvedScalarWave::Tags::BackgroundSpacetime<
+                                       gr::Solutions::KerrSchild>>,
+                 db::AddComputeTags<Tags::GeodesicAccelerationCompute<Dim>>>(
+          std::array<tnsr::I<double, Dim>, 2>{
+              {random_position, random_velocity}},
+          kerr_schild);
+  const auto christoffel = get<
+      gr::Tags::SpacetimeChristoffelSecondKind<double, Dim, Frame::Inertial>>(
+      kerr_schild.variables(random_position, 0.,
+                            tmpl::list<gr::Tags::SpacetimeChristoffelSecondKind<
+                                double, Dim, Frame::Inertial>>{}));
+  const auto expected_acceleration =
+      gr::geodesic_acceleration(random_velocity, christoffel);
+  CHECK_ITERABLE_APPROX(db::get<Tags::GeodesicAcceleration<Dim>>(box),
+                        expected_acceleration);
+}
+
 void test_puncture_field() {
   static constexpr size_t Dim = 3;
   ::TestHelpers::db::test_compute_tag<Tags::PunctureFieldCompute<Dim>>(
@@ -481,10 +512,13 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.Worldtube.Tags",
       Tags::CheckInputFile<3, gr::Solutions::KerrSchild>>("CheckInputFile");
   TestHelpers::db::test_simple_tag<Tags::ObserveCoefficientsTrigger>(
       "ObserveCoefficientsTrigger");
+  TestHelpers::db::test_simple_tag<Tags::GeodesicAcceleration<3>>(
+      "GeodesicAcceleration");
   test_excision_sphere_tag();
   test_compute_face_coordinates_grid();
   test_compute_face_coordinates();
   test_particle_position_velocity_compute();
+  test_geodesic_acceleration_compute();
   test_puncture_field();
   test_check_input_file();
 }

tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt

@@ -7,6 +7,7 @@ set(LIBRARY_SOURCES
   Test_Christoffel.cpp
   Test_ComputeGhQuantities.cpp
   Test_ComputeSpacetimeQuantities.cpp
+  Test_GeodesicAcceleration.cpp
   Test_InterfaceNullNormal.cpp
   Test_KerrSchildCoords.cpp
   Test_ProjectionOperators.cpp