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 spatial_derivative_christoffel( return di_christoffel; } +tnsr::iA spatial_derivative_ks_contracted_christoffel( + const tnsr::I& 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 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 spatial_derivative_christoffel( const tnsr::iiaa& dij_metric, const tnsr::AA& inverse_metric, const tnsr::iAA& 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 spatial_derivative_ks_contracted_christoffel( + const tnsr::I& 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 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>( + 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 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& point) { + const tnsr::I tensor_point(point); + const auto kerr_schild_quantities_local = kerr_schild.variables( + tensor_point, 0., + tmpl::list< + gr::Tags::Lapse, gr::Tags::Shift, + gr::Tags::SpatialMetric, + gr::Tags::InverseSpatialMetric, + gr::Tags::SpacetimeChristoffelSecondKind>{}); + const auto imetric_local = gr::inverse_spacetime_metric( + get>(kerr_schild_quantities_local), + get>(kerr_schild_quantities_local), + get>( + kerr_schild_quantities_local)); + const auto& christoffel_local = + get>( + 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