diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt
index ebcf8ccb673c..54babcdabb0c 100644
--- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt
+++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt
@@ -30,6 +30,8 @@ target_link_libraries(
   PUBLIC
   DataStructures
   Domain
+  GeneralRelativity
+  GeneralRelativitySolutions
   Options
   Parallel
   Utilities
diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.cpp
index fefebad312c5..bd2511f0ceba 100644
--- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.cpp
+++ b/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>;
diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp
index 8d259924b23e..49d662020f9d 100644
--- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp
+++ b/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
diff --git a/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt b/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt
index 00abdb311c26..ddcd59775a6f 100644
--- a/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt
+++ b/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
diff --git a/src/PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.cpp b/src/PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.cpp
new file mode 100644
index 000000000000..9756b3bc8c17
--- /dev/null
+++ b/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);
diff --git a/src/PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp b/src/PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp
new file mode 100644
index 000000000000..ca225f99e545
--- /dev/null
+++ b/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
diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt
index 424cc706654d..aa0562a1322c 100644
--- a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt
+++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt
@@ -41,6 +41,7 @@ target_link_libraries(
   DomainBoundaryConditions
   DomainBoundaryConditionsHelpers
   DomainCreators
+  GeneralRelativity
   GeneralRelativityHelpers
   GeneralRelativitySolutions
   MathFunctions
diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp
index c2d050d43473..7a6182d9c09f 100644
--- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp
+++ b/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();
 }
diff --git a/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt b/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt
index e88223e37a5b..70b23dfc437b 100644
--- a/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt
+++ b/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
diff --git a/tests/Unit/PointwiseFunctions/GeneralRelativity/Test_GeodesicAcceleration.cpp b/tests/Unit/PointwiseFunctions/GeneralRelativity/Test_GeodesicAcceleration.cpp
new file mode 100644
index 000000000000..237b4a5594e9
--- /dev/null
+++ b/tests/Unit/PointwiseFunctions/GeneralRelativity/Test_GeodesicAcceleration.cpp
@@ -0,0 +1,187 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <cstddef>
+#include <random>
+
+#include "DataStructures/BoostMultiArray.hpp"
+#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Framework/CheckWithRandomValues.hpp"
+#include "Helpers/DataStructures/DataBox/TestHelpers.hpp"
+#include "Helpers/PointwiseFunctions/GeneralRelativity/TestHelpers.hpp"
+#include "NumericalAlgorithms/OdeIntegration/OdeIntegration.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
+#include "PointwiseFunctions/GeneralRelativity/GeodesicAcceleration.hpp"
+#include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "Utilities/EqualWithinRoundoff.hpp"
+#include "Utilities/StdArrayHelpers.hpp"
+
+namespace {
+
+template <typename DataType>
+void test_circular_orbit_kerr_schild() {
+  MAKE_GENERATOR(gen);
+  const double mass = 3.2;
+  const std::array<double, 3> spin{{0., 0., 0.}};
+  const std::array<double, 3> center{{0., 0., 0.}};
+  gr::Solutions::KerrSchild kerr_schild(mass, spin, center);
+  const auto used_for_size = Scalar<DataType>(static_cast<size_t>(100));
+  std::uniform_real_distribution<> radius_dist(7., 100.);
+  std::uniform_real_distribution<> angle_dist(0., 2. * M_PI);
+  const auto radius = make_with_random_values<DataType>(
+      make_not_null(&gen), radius_dist, used_for_size);
+  const auto angle = make_with_random_values<DataType>(
+      make_not_null(&gen), angle_dist, used_for_size);
+  const auto angular_velocity = sqrt(mass) / (radius * sqrt(radius));
+
+  const tnsr::I<DataType, 3> position{
+      {radius * cos(angle), radius * sin(angle),
+       make_with_value<DataType>(used_for_size, 0.)}};
+  const tnsr::I<DataType, 3> velocity{
+      {-radius * angular_velocity * sin(angle),
+       radius * angular_velocity * cos(angle),
+       make_with_value<DataType>(used_for_size, 0.)}};
+  const tnsr::I<DataType, 3> expected_acceleration{
+      {-radius * square(angular_velocity) * cos(angle),
+       -radius * square(angular_velocity) * sin(angle),
+       make_with_value<DataType>(used_for_size, 0.)}};
+
+  const auto christoffel = get<
+      gr::Tags::SpacetimeChristoffelSecondKind<DataType, 3, Frame::Inertial>>(
+      kerr_schild.variables(position, 0.,
+                            tmpl::list<gr::Tags::SpacetimeChristoffelSecondKind<
+                                DataType, 3, Frame::Inertial>>{}));
+
+  const auto geodesic_acc = gr::geodesic_acceleration(velocity, christoffel);
+  CHECK_ITERABLE_APPROX(geodesic_acc, expected_acceleration);
+}
+
+std::array<tnsr::I<double, 3>, 2> state_to_tnsr(
+    const std::array<double, 6>& state) {
+  tnsr::I<double, 3> pos_tensor{{state[0], state[1], state[2]}};
+  tnsr::I<double, 3> vel_tensor{{state[3], state[4], state[5]}};
+  return {std::move(pos_tensor), std::move(vel_tensor)};
+}
+
+// struct with interface conforming to a `System` as used in
+// boost::numeric::odeint
+struct BoostGeodesicIntegrator {
+  gr::Solutions::KerrSchild kerr_schild;
+  void operator()(const std::array<double, 6>& state,
+                  std::array<double, 6>& derivative_of_state,
+                  const double /*time*/) const {
+    const auto [position, velocity] = state_to_tnsr(state);
+    const auto christoffel = get<
+        gr::Tags::SpacetimeChristoffelSecondKind<double, 3, Frame::Inertial>>(
+        kerr_schild.variables(
+            position, 0.,
+            tmpl::list<gr::Tags::SpacetimeChristoffelSecondKind<
+                double, 3, Frame::Inertial>>{}));
+    const auto geodesic_acc = gr::geodesic_acceleration(velocity, christoffel);
+    for (size_t i = 0; i < 3; ++i) {
+      gsl::at(derivative_of_state, i) = state.at(i + 3);
+      gsl::at(derivative_of_state, i + 3) = geodesic_acc.get(i);
+    }
+  }
+};
+
+// struct with interface conforming to an `Observer` as used in
+// boost::numeric::odeint
+struct BoostObserver {
+  std::vector<std::array<double, 6>>& states;
+  std::vector<double>& times;
+  void operator()(const std::array<double, 6>& state, double t) {
+    states.push_back(state);
+    times.push_back(t);
+  }
+};
+
+tnsr::a<double, 3> lower_four_velocity(
+    const std::array<double, 6>& state,
+    const gr::Solutions::KerrSchild& kerr_schild) {
+  const auto [pos, vel] = state_to_tnsr(state);
+  const auto spacetime_vars = kerr_schild.variables(
+      pos, 0.,
+      tmpl::list<gr::Tags::Lapse<double>,
+                 gr::Tags::Shift<double, 3, Frame::Inertial>,
+                 gr::Tags::SpatialMetric<double, 3, Frame::Inertial>>{});
+  const auto spacetime_metric = gr::spacetime_metric(
+      get<gr::Tags::Lapse<double>>(spacetime_vars),
+      get<gr::Tags::Shift<double, 3, Frame::Inertial>>(spacetime_vars),
+      get<gr::Tags::SpatialMetric<double, 3, Frame::Inertial>>(spacetime_vars));
+  double temp = spacetime_metric.get(0, 0);
+  for (size_t i = 0; i < 3; ++i) {
+    temp += 2. * spacetime_metric.get(i + 1, 0) * vel.get(i);
+    for (size_t j = 0; j < 3; ++j) {
+      temp += spacetime_metric.get(i + 1, j + 1) * vel.get(i) * vel.get(j);
+    }
+  }
+  // the expression for u0 follows directly from u^i = u^0 * v^i and the
+  // normalization of the four velocity
+  const double u0 = sqrt(-1. / temp);
+  const tnsr::A<double, 3> u{
+      {u0, vel.get(0) * u0, vel.get(1) * u0, vel.get(2) * u0}};
+  return tenex::evaluate<ti::a>(spacetime_metric(ti::a, ti::b) * u(ti::B));
+}
+
+// checks that the conserved quantities due to the spacetime symmetries i.e.
+// Killing vectors are conserved during geodesic evolution
+void test_conserved_quantities_kerr_schild() {
+  const double mass = 1.3;
+  const std::array<double, 3> spin{{0., 0., 0.3}};
+  const std::array<double, 3> center{{0., 0., 0.}};
+  gr::Solutions::KerrSchild kerr_schild(mass, spin, center);
+
+  const double t_max = 1000.0;
+
+  // corresponds to a strongly perturbed circular orbit
+  std::array<double, 6> initial_state = {
+      10., 1.0, 0.8, 0.1, 1. / sqrt(10.) * 1.2, 0.3};
+
+  std::vector<std::array<double, 6>> states{};
+  std::vector<double> times{};
+  BoostObserver observer{states, times};
+  boost::numeric::odeint::integrate_adaptive(
+      boost::numeric::odeint::make_controlled(
+          1e-15, 1e-15,
+          boost::numeric::odeint::runge_kutta_cash_karp54<
+              std::array<double, 6>>()),
+      BoostGeodesicIntegrator{kerr_schild}, initial_state, 0.0, t_max, 1e-5,
+      observer);
+
+  const auto initial_four_velocity =
+      lower_four_velocity(initial_state, kerr_schild);
+
+  // the two Killing vectors in Kerr-Schild coordinates
+  const tnsr::A<double, 3> timelike_killing{{1., 0., 0., 0.}};
+  const tnsr::A<double, 3> azimuthal_killing_initial{
+      {0., initial_state[1], -initial_state[0], 0.}};
+
+  const double initial_energy =
+      get(dot_product(initial_four_velocity, timelike_killing));
+  const double initial_angular_momentum =
+      get(dot_product(initial_four_velocity, azimuthal_killing_initial));
+
+  for (const auto& state : states) {
+    const auto four_velocity = lower_four_velocity(state, kerr_schild);
+    const double energy = get(dot_product(four_velocity, timelike_killing));
+    const tnsr::A<double, 3> azimuthal_killing{{0., state[1], -state[0], 0.}};
+    const double angular_momentum =
+        get(dot_product(four_velocity, azimuthal_killing));
+    CHECK(initial_energy == approx(energy));
+    CHECK(initial_angular_momentum == approx(angular_momentum));
+  }
+}
+}  // namespace
+
+SPECTRE_TEST_CASE(
+    "Unit.PointwiseFunctions.GeneralRelativity.GeodesicAcceleration",
+    "[Unit][PointwiseFunctions]") {
+  test_circular_orbit_kerr_schild<double>();
+  test_circular_orbit_kerr_schild<DataVector>();
+  test_conserved_quantities_kerr_schild();
+}