From 327d83da6dec243e674cfda2868c6726b0b54615 Mon Sep 17 00:00:00 2001
From: Nikolas <niko.wittek@gmail.com>
Date: Mon, 18 Mar 2024 21:32:50 +0100
Subject: [PATCH] Add derivative of contracted KS Christoffel

---
 .../Worldtube/KerrSchildDerivatives.cpp       | 22 ++++++++
 .../Worldtube/KerrSchildDerivatives.hpp       | 10 ++++
 .../Worldtube/Test_KerrSchildDerivatives.cpp  | 51 +++++++++++++++++++
 3 files changed, 83 insertions(+)

diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp
index 5a4ea2664a31..ca242bb9907f 100644
--- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp
+++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp
@@ -166,4 +166,26 @@ tnsr::iAbb<double, 3> spatial_derivative_christoffel(
   return di_christoffel;
 }
 
+tnsr::iA<double, 3> spatial_derivative_ks_contracted_christoffel(
+    const tnsr::I<double, 3>& pos) {
+  const double r_sq = get(dot_product(pos, pos));
+  const double r = sqrt(r_sq);
+  const double one_over_r = 1. / r;
+  const double one_over_r_2 = 1. / r_sq;
+  const double one_over_r_3 = cube(one_over_r);
+  const double one_over_r_4 = square(one_over_r_2);
+  const double one_over_r_5 = one_over_r_4 * one_over_r;
+
+  tnsr::iA<double, 3> di_contracted_christoffel{};
+  for (size_t i = 0; i < 3; ++i) {
+    di_contracted_christoffel.get(i, 0) = 4. * pos.get(i) * one_over_r_4;
+    for (size_t j = 0; j < 3; ++j) {
+      di_contracted_christoffel.get(i, j + 1) =
+          -6. * pos.get(i) * pos.get(j) * one_over_r_5;
+    }
+    di_contracted_christoffel.get(i, i + 1) += 2. * one_over_r_3;
+  }
+  return di_contracted_christoffel;
+}
+
 }  // namespace CurvedScalarWave::Worldtube
diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp
index c4170ee43cd8..322477429ee7 100644
--- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp
+++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp
@@ -51,4 +51,14 @@ tnsr::iAbb<double, 3> spatial_derivative_christoffel(
     const tnsr::iiaa<double, 3>& dij_metric,
     const tnsr::AA<double, 3>& inverse_metric,
     const tnsr::iAA<double, 3>& di_inverse_metric);
+
+/*!
+ * \brief The spatial derivative of the zero spin Kerr Schild contracted
+ * Christoffel symbols,
+ * $\partial_i g^{\mu \nu} \Gamma^\rho_{\mu \nu}$, assuming a black hole at the
+ * coordinate center with mass M = 1.
+ */
+tnsr::iA<double, 3> spatial_derivative_ks_contracted_christoffel(
+    const tnsr::I<double, 3>& pos);
+
 }  // namespace CurvedScalarWave::Worldtube
diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp
index 496b228cdbfe..9bc154b17d1d 100644
--- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp
+++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp
@@ -293,6 +293,56 @@ void test_derivative_christoffel() {
   }
 }
 
+void test_derivative_contracted_christoffel() {
+  static constexpr size_t Dim = 3;
+  MAKE_GENERATOR(gen);
+  Approx local_approx = Approx::custom().epsilon(1e-10).scale(1.0);
+  std::uniform_real_distribution<double> pos_dist{2., 10.};
+
+  const gr::Solutions::KerrSchild kerr_schild(1., {{0., 0., 0.}},
+                                              {{0., 0., 0.}});
+
+  const auto test_point =
+      make_with_random_values<tnsr::I<double, 3, Frame::Inertial>>(
+          make_not_null(&gen), make_not_null(&pos_dist), 1);
+  const auto di_contracted_christoffel =
+      spatial_derivative_ks_contracted_christoffel(test_point);
+
+  const std::array<double, 3> array_point{test_point.get(0), test_point.get(1),
+                                          test_point.get(2)};
+  const double dx = 1e-4;
+  for (size_t a = 0; a < 4; ++a) {
+    for (size_t j = 0; j < 3; ++j) {
+      const auto contracted_christoffel_helper =
+          [&kerr_schild, a](const std::array<double, 3>& point) {
+            const tnsr::I<double, 3> tensor_point(point);
+            const auto kerr_schild_quantities_local = kerr_schild.variables(
+                tensor_point, 0.,
+                tmpl::list<
+                    gr::Tags::Lapse<double>, gr::Tags::Shift<double, Dim>,
+                    gr::Tags::SpatialMetric<double, Dim>,
+                    gr::Tags::InverseSpatialMetric<double, Dim>,
+                    gr::Tags::SpacetimeChristoffelSecondKind<double, Dim>>{});
+            const auto imetric_local = gr::inverse_spacetime_metric(
+                get<gr::Tags::Lapse<double>>(kerr_schild_quantities_local),
+                get<gr::Tags::Shift<double, Dim>>(kerr_schild_quantities_local),
+                get<gr::Tags::InverseSpatialMetric<double, Dim>>(
+                    kerr_schild_quantities_local));
+            const auto& christoffel_local =
+                get<gr::Tags::SpacetimeChristoffelSecondKind<double, Dim>>(
+                    kerr_schild_quantities_local);
+            const auto contracted_christoffel_local =
+                trace_last_indices(christoffel_local, imetric_local);
+            return contracted_christoffel_local.get(a);
+          };
+      const auto numerical_deriv_contracted_christoffel_j =
+          numerical_derivative(contracted_christoffel_helper, array_point, j,
+                               dx);
+      CHECK(di_contracted_christoffel.get(j, a) ==
+            local_approx(numerical_deriv_contracted_christoffel_j));
+    }
+  }
+}
 
 // All the derivatives are checked against finite difference calculations
 SPECTRE_TEST_CASE(
@@ -303,6 +353,7 @@ SPECTRE_TEST_CASE(
   test_second_derivative_inverse_metric();
   test_second_derivative_metric();
   test_derivative_christoffel();
+  test_derivative_contracted_christoffel();
 }
 }  // namespace
 }  // namespace CurvedScalarWave::Worldtube