Add geodesic acceleration compute tag

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

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();
 }