diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt index 9e6be3f60d7f3..b1a05c165fbd4 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt @@ -16,6 +16,7 @@ set(LIBRARY_SOURCES Test_PsiSquared.cpp Test_Tags.cpp Test_TimeDerivative.cpp + Worldtube/Test_PunctureField.cpp Worldtube/Test_Tags.cpp Worldtube/SingletonActions/Test_InitializeElementFacesGridCoordinates.cpp ) @@ -39,6 +40,7 @@ target_link_libraries( GeneralRelativitySolutions MathFunctions GeneralRelativitySolutions + ScalarWaveWorldtube Utilities WaveEquationSolutions ) diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_PunctureField.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_PunctureField.cpp new file mode 100644 index 0000000000000..ad81c8748d370 --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_PunctureField.cpp @@ -0,0 +1,141 @@ +// 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/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Evolution/Systems/CurvedScalarWave/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/PunctureField.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" +#include "Utilities/Gsl.hpp" + +namespace CurvedScalarWave { +namespace { + +using deriv_psi_tag = + ::Tags::deriv, Frame::Inertial>; +using puncture_vars = + Variables, deriv_psi_tag>>; +// tests that the puncture field corresponds to a circular orbit with angular +// velocity = orbit_radius^(-3/2) +void test_circular_orbit() { + MAKE_GENERATOR(gen) + // sample 100 random points around the worldtube + const double orbit_radius = 7.; + const double orbit_speed = pow(orbit_radius, -1.5); + const size_t num_points = 100; + const tnsr::I wt_coords{{orbit_radius, 0., 0.}}; + std::uniform_real_distribution dist_around_wt(-3., 3.); + auto sample_points = + make_with_random_values>( + make_not_null(&gen), make_not_null(&dist_around_wt), + DataVector(num_points)); + sample_points.get(0) += orbit_radius; + const double time_0 = 0.; + const double time_1 = 20.; + for (size_t order = 0; order < 1; ++order) { + CAPTURE(order); + puncture_vars puncture_t0{num_points}; + Worldtube::compute_puncture_field(make_not_null(&puncture_t0), + sample_points, time_0, orbit_radius, 1., + order); + + // rotate the sample points and check that the values don't change + tnsr::I sample_points_rotated(num_points, + 0.); + sample_points_rotated.get(0) = + sample_points.get(0) * cos(orbit_speed * time_1) - + sample_points.get(1) * sin(orbit_speed * time_1); + sample_points_rotated.get(1) = + sample_points.get(0) * sin(orbit_speed * time_1) + + sample_points.get(1) * cos(orbit_speed * time_1); + sample_points_rotated.get(2) = sample_points.get(2); + puncture_vars puncture_t1{num_points}; + Worldtube::compute_puncture_field(make_not_null(&puncture_t1), + sample_points_rotated, time_1, + orbit_radius, 1., order); + CHECK_ITERABLE_APPROX(get(puncture_t0).get(), + get(puncture_t1).get()); + CHECK_ITERABLE_APPROX(get<::Tags::dt>(puncture_t0).get(), + get<::Tags::dt>(puncture_t1).get()); + + // check that the spatial derivative also gets rotated + const auto& di_psi = get(puncture_t0); + tnsr::i di_psi_rotated(num_points); + di_psi_rotated.get(0) = di_psi.get(0) * cos(orbit_speed * time_1) - + di_psi.get(1) * sin(orbit_speed * time_1); + di_psi_rotated.get(1) = di_psi.get(0) * sin(orbit_speed * time_1) + + di_psi.get(1) * cos(orbit_speed * time_1); + di_psi_rotated.get(2) = di_psi.get(2); + + CHECK_ITERABLE_APPROX(di_psi_rotated, get(puncture_t1)); + } +} + +// tests the derivative of the puncture field against a finite difference +// calculation +void test_derivative() { + MAKE_GENERATOR(gen); + const double time = 17.4; + const double orbit_radius = 10.; + for (size_t order = 0; order < 1; ++order) { + CAPTURE(order); + const auto helper_func = [&time, &orbit_radius, + &order](const std::array& point) { + tnsr::I tensor_point(size_t(1)); + tensor_point.get(0) = point.at(0); + tensor_point.get(1) = point.at(1); + tensor_point.get(2) = point.at(2); + puncture_vars puncture_field{1}; + Worldtube::compute_puncture_field(make_not_null(&puncture_field), + tensor_point, time, orbit_radius, 1., + order); + return get(puncture_field).get()[0]; + }; + Approx local_approx = Approx::custom().epsilon(1e-10).scale(1.0); + std::uniform_real_distribution theta_dist{0, M_PI}; + std::uniform_real_distribution phi_dist{0, 2 * M_PI}; + const double wt_radius = 0.1; + for (size_t i = 0; i < 1000; ++i) { + const auto theta = theta_dist(gen); + const auto phi = phi_dist(gen); + + const std::array test_point{ + orbit_radius + wt_radius * sin(theta) * cos(phi), + wt_radius * sin(theta) * sin(phi), wt_radius * cos(theta)}; + const double dx = 1e-4; + tnsr::I tensor_point(size_t(1)); + tensor_point.get(0) = orbit_radius + wt_radius * sin(theta) * cos(phi); + tensor_point.get(1) = wt_radius * sin(theta) * sin(phi); + tensor_point.get(2) = wt_radius * cos(theta); + puncture_vars puncture_field{1}; + Worldtube::compute_puncture_field(make_not_null(&puncture_field), + tensor_point, time, orbit_radius, 1., + order); + const auto& di_psi = get(puncture_field); + for (size_t j = 0; j < 3; ++j) { + const auto numerical_deriv_j = + numerical_derivative(helper_func, test_point, j, dx); + CHECK(di_psi.get(j)[0] == local_approx(numerical_deriv_j)); + } + } + } +} + +SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.PunctureField", + "[Unit][Evolution]") { + test_circular_orbit(); + test_derivative(); +} +} // namespace +} // namespace CurvedScalarWave