diff --git a/src/DataStructures/Tensor/TypeAliases.hpp b/src/DataStructures/Tensor/TypeAliases.hpp index 647bd601c4ad..9b226e9bd877 100644 --- a/src/DataStructures/Tensor/TypeAliases.hpp +++ b/src/DataStructures/Tensor/TypeAliases.hpp @@ -359,6 +359,30 @@ using ijaa = Tensor, SpacetimeIndex, SpacetimeIndex>>; template +using iiaa = Tensor, + index_list, + SpatialIndex, + SpacetimeIndex, + SpacetimeIndex>>; +template +using iiAA = Tensor, + index_list, + SpatialIndex, + SpacetimeIndex, + SpacetimeIndex>>; +template +using iAbb = Tensor, + index_list, + SpacetimeIndex, + SpacetimeIndex, + SpacetimeIndex>>; +template +using iabb = Tensor, + index_list, + SpacetimeIndex, + SpacetimeIndex, + SpacetimeIndex>>; +template using Ijaa = Tensor, index_list, SpatialIndex, diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 35ea0c1e5727..7b34dc4a3ce3 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -8,6 +8,7 @@ add_spectre_library(${LIBRARY}) spectre_target_sources( ${LIBRARY} PRIVATE + KerrSchildDerivatives.cpp Tags.cpp PunctureField.cpp PunctureFieldOrder0.cpp @@ -20,6 +21,7 @@ spectre_target_headers( INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src HEADERS Inboxes.hpp + KerrSchildDerivatives.hpp SingletonChare.hpp Tags.hpp PunctureField.hpp diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp new file mode 100644 index 000000000000..ca242bb9907f --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.cpp @@ -0,0 +1,191 @@ + +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp" + +#include + +#include "DataStructures/Tensor/EagerMath/DotProduct.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/EagerMath/Trace.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/Gsl.hpp" + +namespace CurvedScalarWave::Worldtube { +tnsr::iAA spatial_derivative_inverse_ks_metric( + 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 = one_over_r_2 * one_over_r; + + tnsr::iAA di_imetric{}; + tnsr::ii delta_ll{0.}; + tnsr::Ij delta_ul{0.}; + tnsr::i pos_lower{}; + + for (size_t i = 0; i < 3; ++i) { + delta_ll.get(i, i) = 1.; + delta_ul.get(i, i) = 1.; + pos_lower.get(i) = pos.get(i); + } + + const auto d_imetric_ij = tenex::evaluate( + one_over_r_3 * + (6. * pos(ti::J) * pos(ti::K) * pos_lower(ti::i) * one_over_r_2 - + 2. * delta_ul(ti::J, ti::i) * pos(ti::K) - + 2. * delta_ul(ti::K, ti::i) * pos(ti::J))); + const auto d_imetric_i0 = tenex::evaluate( + one_over_r_2 * (-4. * pos_lower(ti::i) * pos(ti::J) * one_over_r_2 + + 2. * delta_ul(ti::J, ti::i))); + const auto d_imetric_00 = + tenex::evaluate(2. * pos_lower(ti::i) * one_over_r_3); + for (size_t i = 0; i < 3; ++i) { + di_imetric.get(i, 0, 0) = d_imetric_00.get(i); + for (size_t j = 0; j < 3; ++j) { + di_imetric.get(i, j + 1, 0) = d_imetric_i0.get(i, j); + for (size_t k = 0; k < 3; ++k) { + di_imetric.get(i, j + 1, k + 1) = d_imetric_ij.get(i, j, k); + } + } + } + return di_imetric; +} + +tnsr::iaa spatial_derivative_ks_metric( + const tnsr::aa& metric, + const tnsr::iAA& di_inverse_metric) { + tnsr::iaa di_metric{}; + tenex::evaluate( + make_not_null(&di_metric), -metric(ti::a, ti::c) * metric(ti::b, ti::d) * + di_inverse_metric(ti::i, ti::C, ti::D)); + return di_metric; +} + +tnsr::iiAA second_spatial_derivative_inverse_ks_metric( + 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 = one_over_r_2 * one_over_r; + const double one_over_r_4 = one_over_r_2 * one_over_r_2; + + tnsr::iiAA dij_imetric{}; + tnsr::ii delta_ll{0.}; + tnsr::Ij delta_ul{0.}; + tnsr::i pos_lower{}; + + for (size_t i = 0; i < 3; ++i) { + delta_ll.get(i, i) = 1.; + delta_ul.get(i, i) = 1.; + pos_lower.get(i) = pos.get(i); + } + + const auto d2_imetric_ij = tenex::evaluate( + one_over_r_3 * + (-2. * (delta_ul(ti::L, ti::i) * delta_ul(ti::K, ti::j) + + delta_ul(ti::K, ti::i) * delta_ul(ti::L, ti::j)) + + one_over_r_2 * + (6. * (delta_ll(ti::i, ti::j) * pos(ti::K) * pos(ti::L) + + delta_ul(ti::K, ti::i) * pos_lower(ti::j) * pos(ti::L) + + delta_ul(ti::K, ti::j) * pos_lower(ti::i) * pos(ti::L) + + delta_ul(ti::L, ti::i) * pos_lower(ti::j) * pos(ti::K) + + delta_ul(ti::L, ti::j) * pos_lower(ti::i) * pos(ti::K)) - + one_over_r_2 * 30. * pos_lower(ti::i) * pos_lower(ti::j) * + pos(ti::K) * pos(ti::L)))); + + const auto d2_imetric_i0 = tenex::evaluate( + one_over_r_4 * + (-4. * (delta_ll(ti::k, ti::j) * pos(ti::I) + + delta_ul(ti::I, ti::k) * pos_lower(ti::j) + + delta_ul(ti::I, ti::j) * pos_lower(ti::k)) + + one_over_r_2 * 16. * pos(ti::I) * pos_lower(ti::j) * pos_lower(ti::k))); + const auto d2_imetric_00 = tenex::evaluate( + one_over_r_3 * (2. * delta_ll(ti::i, ti::j) - + one_over_r_2 * 6. * pos_lower(ti::i) * pos_lower(ti::j))); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 3; ++j) { + dij_imetric.get(i, j, 0, 0) = d2_imetric_00.get(i, j); + for (size_t k = 0; k < 3; ++k) { + dij_imetric.get(i, j, k + 1, 0) = d2_imetric_i0.get(i, j, k); + for (size_t l = 0; l < 3; ++l) { + dij_imetric.get(i, j, k + 1, l + 1) = d2_imetric_ij.get(i, j, k, l); + } + } + } + } + return dij_imetric; +} + +tnsr::iiaa second_spatial_derivative_metric( + const tnsr::aa& metric, const tnsr::iaa& di_metric, + const tnsr::iAA& di_inverse_metric, + const tnsr::iiAA& dij_inverse_metric) { + tnsr::iiaa dij_metric{}; + tenex::evaluate( + make_not_null(&dij_metric), + -metric(ti::a, ti::c) * metric(ti::b, ti::d) * + dij_inverse_metric(ti::j, ti::i, ti::C, ti::D) - + 2. * metric(ti::a, ti::c) * di_metric(ti::j, ti::b, ti::d) * + di_inverse_metric(ti::i, ti::C, ti::D)); + return dij_metric; +} + +tnsr::iAbb spatial_derivative_christoffel( + const tnsr::iaa& di_metric, + const tnsr::iiaa& dij_metric, + const tnsr::AA& inverse_metric, + const tnsr::iAA& di_inverse_metric) { + tnsr::iAbb di_christoffel{}; + tnsr::abb d_metric{}; + tnsr::iabb di_d_metric{}; + for (size_t a = 0; a <= 3; ++a) { + for (size_t b = 0; b <= 3; ++b) { + d_metric.get(0, a, b) = 0.; + for (size_t i = 0; i < 3; ++i) { + d_metric.get(i + 1, a, b) = di_metric.get(i, a, b); + di_d_metric.get(i, 0, a, b) = 0.; + for (size_t j = 0; j < 3; ++j) { + di_d_metric.get(i, j + 1, a, b) = dij_metric.get(i, j, a, b); + } + } + } + } + tenex::evaluate( + make_not_null(&di_christoffel), + 0.5 * di_inverse_metric(ti::i, ti::A, ti::D) * + (d_metric(ti::b, ti::c, ti::d) + d_metric(ti::c, ti::b, ti::d) - + d_metric(ti::d, ti::b, ti::c)) + + 0.5 * inverse_metric(ti::A, ti::D) * + (di_d_metric(ti::i, ti::b, ti::c, ti::d) + + di_d_metric(ti::i, ti::c, ti::b, ti::d) - + di_d_metric(ti::i, ti::d, ti::b, ti::c))); + 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 new file mode 100644 index 000000000000..322477429ee7 --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp @@ -0,0 +1,64 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/Gsl.hpp" + +namespace CurvedScalarWave::Worldtube { +/*! + * \brief The spatial derivative of the zero spin inverse Kerr Schild metric, + * $\partial_i g^{\mu \nu}$, assuming a black hole at the coordinate center with + * mass M = 1. + */ +tnsr::iAA spatial_derivative_inverse_ks_metric( + const tnsr::I& pos); + +/*! + * \brief The spatial derivative of the spacetime metric, + * $\partial_i g_{\mu \nu}$. + */ +tnsr::iaa spatial_derivative_ks_metric( + const tnsr::aa& metric, + const tnsr::iAA& di_inverse_metric); + +/*! + * \brief The second spatial derivative of the zero spin inverse Kerr Schild + * metric, $\partial_i \partial_j g^{\mu \nu}$, assuming a black hole at the + * coordinate center with mass M = 1. + */ +tnsr::iiAA second_spatial_derivative_inverse_ks_metric( + const tnsr::I& pos); + +/*! + * \brief The spatial derivative of the spacetime metric, + * $\partial_i \partial_j g_{\mu \nu}$. + */ +tnsr::iiaa second_spatial_derivative_metric( + const tnsr::aa& metric, const tnsr::iaa& di_metric, + const tnsr::iAA& di_inverse_metric, + const tnsr::iiAA& dij_inverse_metric); + +/*! + * \brief The spatial derivative of the Christoffel + * symbols, $\partial_i \Gamma^\rho_{\mu \nu}$. + */ +tnsr::iAbb spatial_derivative_christoffel( + const tnsr::iaa& di_metric, + 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/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 64cecd9b5328..aeefc2156d13 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -4,6 +4,7 @@ set(LIBRARY "Test_ScalarWaveWorldtube") set(LIBRARY_SOURCES + Test_KerrSchildDerivatives.cpp Test_PunctureField.cpp Test_SelfForce.cpp Test_Tags.cpp diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp new file mode 100644 index 000000000000..9bc154b17d1d --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_KerrSchildDerivatives.cpp @@ -0,0 +1,359 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include + +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/EagerMath/Trace.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Evolution/Systems/CurvedScalarWave/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/GeneralRelativity/InverseSpacetimeMetric.hpp" +#include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "Utilities/Gsl.hpp" + +namespace CurvedScalarWave::Worldtube { +namespace { + +void test_derivative_inverse_metric() { + 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_imetric = spatial_derivative_inverse_ks_metric(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 b = 0; b <= a; ++b) { + for (size_t j = 0; j < 3; ++j) { + auto inverse_metric_helper = [&kerr_schild, a, + b](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::Shift, + gr::Tags::InverseSpatialMetric>{}); + const auto inverse_metric_local = gr::inverse_spacetime_metric( + get>(kerr_schild_quantities_local), + get>(kerr_schild_quantities_local), + get>( + kerr_schild_quantities_local)); + return inverse_metric_local.get(a, b); + }; + const auto numerical_deriv_imetric_j = + numerical_derivative(inverse_metric_helper, array_point, j, dx); + CHECK(di_imetric.get(j, a, b) == + local_approx(numerical_deriv_imetric_j)); + } + } + } +} + +void test_derivative_metric() { + 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 kerr_schild_quantities = kerr_schild.variables( + test_point, 0., + tmpl::list, gr::Tags::Shift, + gr::Tags::SpatialMetric, + gr::Tags::SpacetimeChristoffelSecondKind>{}); + const auto metric = gr::spacetime_metric( + get>(kerr_schild_quantities), + get>(kerr_schild_quantities), + get>(kerr_schild_quantities)); + + const auto& di_imetric = spatial_derivative_inverse_ks_metric(test_point); + const auto di_metric = spatial_derivative_ks_metric(metric, di_imetric); + + 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 b = 0; b <= a; ++b) { + for (size_t j = 0; j < 3; ++j) { + const auto metric_helper = [&kerr_schild, a, + b](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::Shift, + gr::Tags::SpatialMetric>{}); + const auto metric_local = gr::spacetime_metric( + get>(kerr_schild_quantities_local), + get>(kerr_schild_quantities_local), + get>( + kerr_schild_quantities_local)); + return metric_local.get(a, b); + }; + const auto numerical_deriv_metric_j = + numerical_derivative(metric_helper, array_point, j, dx); + CHECK(di_metric.get(j, a, b) == local_approx(numerical_deriv_metric_j)); + } + } + } +} + +void test_second_derivative_inverse_metric() { + 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& dij_imetric = + second_spatial_derivative_inverse_ks_metric(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 b = 0; b <= a; ++b) { + for (size_t j = 0; j < 3; ++j) { + for (size_t k = 0; k <= j; ++k) { + const auto di_imetric_helper = + [a, b, j](const std::array& point) { + const tnsr::I tensor_point(point); + const auto di_imetric_local = + spatial_derivative_inverse_ks_metric(tensor_point); + return di_imetric_local.get(j, a, b); + }; + const auto second_numerical_deriv_imetric_k = + numerical_derivative(di_imetric_helper, array_point, k, dx); + CHECK(dij_imetric.get(k, j, a, b) == + local_approx(second_numerical_deriv_imetric_k)); + } + } + } + } +} + +void test_second_derivative_metric() { + 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 kerr_schild_quantities = kerr_schild.variables( + test_point, 0., + tmpl::list, gr::Tags::Shift, + gr::Tags::SpatialMetric, + gr::Tags::SpacetimeChristoffelSecondKind>{}); + const auto metric = gr::spacetime_metric( + get>(kerr_schild_quantities), + get>(kerr_schild_quantities), + get>(kerr_schild_quantities)); + + const auto& di_imetric = spatial_derivative_inverse_ks_metric(test_point); + const auto& dij_imetric = + second_spatial_derivative_inverse_ks_metric(test_point); + const auto di_metric = spatial_derivative_ks_metric(metric, di_imetric); + const auto& dij_metric = second_spatial_derivative_metric( + metric, di_metric, di_imetric, dij_imetric); + + 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 b = 0; b <= a; ++b) { + for (size_t j = 0; j < 3; ++j) { + for (size_t k = 0; k <= j; ++k) { + const auto di_metric_helper = [&kerr_schild, j, a, + b](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::SpacetimeChristoffelSecondKind>{}); + const auto metric_local = gr::spacetime_metric( + get>(kerr_schild_quantities_local), + get>(kerr_schild_quantities_local), + get>( + kerr_schild_quantities_local)); + const auto di_imetric_local = + spatial_derivative_inverse_ks_metric(tensor_point); + const auto di_metric_local = + spatial_derivative_ks_metric(metric_local, di_imetric_local); + return di_metric_local.get(j, a, b); + }; + + const auto second_numerical_deriv_metric_k = + numerical_derivative(di_metric_helper, array_point, k, dx); + CHECK(dij_metric.get(k, j, a, b) == + local_approx(second_numerical_deriv_metric_k)); + } + } + } + } +} + +void test_derivative_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 kerr_schild_quantities = kerr_schild.variables( + test_point, 0., + tmpl::list, gr::Tags::Shift, + gr::Tags::SpatialMetric, + gr::Tags::InverseSpatialMetric, + gr::Tags::SpacetimeChristoffelSecondKind>{}); + const auto metric = gr::spacetime_metric( + get>(kerr_schild_quantities), + get>(kerr_schild_quantities), + get>(kerr_schild_quantities)); + const auto inverse_metric = gr::inverse_spacetime_metric( + get>(kerr_schild_quantities), + get>(kerr_schild_quantities), + get>(kerr_schild_quantities)); + + const auto& di_imetric = spatial_derivative_inverse_ks_metric(test_point); + const auto& dij_imetric = + second_spatial_derivative_inverse_ks_metric(test_point); + const auto di_metric = spatial_derivative_ks_metric(metric, di_imetric); + const auto& dij_metric = second_spatial_derivative_metric( + metric, di_metric, di_imetric, dij_imetric); + const auto di_christoffel = spatial_derivative_christoffel( + di_metric, dij_metric, inverse_metric, di_imetric); + + 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 b = 0; b < 4; ++b) { + for (size_t c = 0; c < 4; ++c) { + for (size_t j = 0; j < 3; ++j) { + const auto christoffel_helper = + [&kerr_schild, a, b, c](const std::array& point) { + const tnsr::I tensor_point(point); + const auto kerr_schild_quantities_local = kerr_schild.variables( + tensor_point, 0., + tmpl::list>{}); + const auto& christoffel_local = + get>( + kerr_schild_quantities_local); + return christoffel_local.get(a, b, c); + }; + + const auto numerical_deriv_christoffel_j = + numerical_derivative(christoffel_helper, array_point, j, dx); + CHECK(di_christoffel.get(j, a, b, c) == + local_approx(numerical_deriv_christoffel_j)); + } + } + } + } +} + +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( + "Unit.Evolution.Systems.CurvedScalarWave.KerrSchildDerivatives", + "[Unit][Evolution]") { + test_derivative_inverse_metric(); + test_derivative_metric(); + test_second_derivative_inverse_metric(); + test_second_derivative_metric(); + test_derivative_christoffel(); + test_derivative_contracted_christoffel(); +} +} // namespace +} // namespace CurvedScalarWave::Worldtube