diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index fd839cd320326..99bebd7da7775 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -36,7 +36,6 @@ spectre_target_headers( target_link_libraries( ${LIBRARY} PUBLIC - ControlSystem DataStructures Domain GeneralRelativity diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp index 857ec07cc3d21..71c36f79a55b9 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp @@ -33,7 +33,7 @@ struct InitializeEvolvedVariables { using simple_tags = tmpl::list>; using return_tags = simple_tags; @@ -41,9 +41,9 @@ struct InitializeEvolvedVariables { using const_global_cache_tags = tmpl::list<>; using mutable_global_cache_tags = tmpl::list<>; using simple_tags_from_options = tmpl::list; - using argument_tags = - tmpl::list<::Tags::TimeStepper, - Tags::InitialPositionAndVelocity, ::Tags::Time>; + using argument_tags = tmpl::list<::Tags::TimeStepper, + Tags::InitialPositionAndVelocity, + ::Tags::Time, Tags::ExcisionSphere>; static void apply( const gsl::not_null, Tags::EvolvedVelocity>>*> @@ -54,13 +54,18 @@ struct InitializeEvolvedVariables { dt_evolved_vars, const gsl::not_null current_iteration, const gsl::not_null expiration_time, + const gsl::not_null worldtube_radius, const gsl::not_null<::Tags::HistoryEvolvedVariables::type*> time_stepper_history, const TimeStepper& time_stepper, const std::array, 2>& initial_pos_and_vel, - const double initial_time) { + const double initial_time, const ExcisionSphere& excision_sphere) { *current_iteration = 0; + + // the functions of time should be ready during the first step but not the + // second. We choose the arbitrary value of 1e-10 here to ensure this. *expiration_time = initial_time + 1e-10; + *worldtube_radius = excision_sphere.radius(); const size_t starting_order = time_stepper.number_of_past_steps() == 0 ? time_stepper.order() : 1; diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/ReceiveElementData.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/ReceiveElementData.hpp index 5c228a52cced7..922158b881dd9 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/ReceiveElementData.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/ReceiveElementData.hpp @@ -102,7 +102,7 @@ struct ReceiveElementData { for (const auto& [_, element_ylm_coefs] : inbox.at(time_step_id)) { external_ylm_coefs += element_ylm_coefs; } - const double wt_radius = db::get>(box).radius(); + const double wt_radius = db::get(box); external_ylm_coefs /= wt_radius * wt_radius; DataVector& psi_ylm_coefs = diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateFunctionsOfTime.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateFunctionsOfTime.hpp index 42d35c60760ba..928a4a69a57b7 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateFunctionsOfTime.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateFunctionsOfTime.hpp @@ -26,6 +26,11 @@ #include "Utilities/Gsl.hpp" #include "Utilities/TMPL.hpp" +namespace control_system { +// forward declare +struct UpdateMultipleFunctionsOfTime; +} // namespace control_system + namespace CurvedScalarWave::Worldtube::Actions { /*! @@ -61,15 +66,10 @@ struct UpdateFunctionsOfTime { const double grid_radius_particle = get(magnitude(excision_sphere.center())); - DataVector angular_update(2); - DataVector expansion_update(2); - DataVector size_a_update(2); - DataVector size_b_update(2); - - angular_update[0] = atan2(y, x); - angular_update[1] = (x * ydot - y * xdot) / square(r); - expansion_update[0] = r / grid_radius_particle; - expansion_update[1] = radial_vel / grid_radius_particle; + const DataVector angular_update{atan2(y, x), + (x * ydot - y * xdot) / square(r)}; + const DataVector expansion_update{r / grid_radius_particle, + radial_vel / grid_radius_particle}; const auto& wt_radius_params = db::get(box); @@ -101,21 +101,20 @@ struct UpdateFunctionsOfTime { const double sqrt_4_pi = sqrt(4. * M_PI); // The expansion map has already stretched the excision spheres and we need // to account for that. - size_a_update[0] = - sqrt_4_pi * (wt_radius_grid - wt_radius_inertial / expansion_update[0]); - size_a_update[1] = sqrt_4_pi * - (wt_radius_inertial * expansion_update[1] - - wt_radius_time_derivative * expansion_update[0]) / - square(expansion_update[0]); + const DataVector size_a_update{ + sqrt_4_pi * (wt_radius_grid - wt_radius_inertial / expansion_update[0]), + sqrt_4_pi * + (wt_radius_inertial * expansion_update[1] - + wt_radius_time_derivative * expansion_update[0]) / + square(expansion_update[0])}; - size_b_update[0] = + const DataVector size_b_update{ sqrt_4_pi * (bh_excision_radius_grid - - bh_excision_radius_inertial / expansion_update[0]); - size_b_update[1] = + bh_excision_radius_inertial / expansion_update[0]), sqrt_4_pi * - (bh_excision_radius_inertial * expansion_update[1] - - bh_excision_radius_time_derivative * expansion_update[0]) / - square(expansion_update[0]); + (bh_excision_radius_inertial * expansion_update[1] - + bh_excision_radius_time_derivative * expansion_update[0]) / + square(expansion_update[0])}; // we set the new expiration time to 1/100th of a time step next to the // current time. This is small enough that it can handle rapid time step @@ -138,12 +137,11 @@ struct UpdateFunctionsOfTime { std::unordered_map> all_updates{}; all_updates["Rotation"] = - std::make_pair(std::move(angular_update), new_expiration_time); + std::make_pair(angular_update, new_expiration_time); all_updates["Expansion"] = - std::make_pair(std::move(expansion_update), new_expiration_time); + std::make_pair(expansion_update, new_expiration_time); all_updates["SizeA"] = std::make_pair(size_a_update, new_expiration_time); - all_updates["SizeB"] = - std::make_pair(std::move(size_b_update), new_expiration_time); + all_updates["SizeB"] = std::make_pair(size_b_update, new_expiration_time); Parallel::mutate<::domain::Tags::FunctionsOfTime, control_system::UpdateMultipleFunctionsOfTime>( diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 9c8d3e7f077eb..3b87f4f0030ec 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -28,6 +28,7 @@ add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") target_link_libraries( ${LIBRARY} PRIVATE + ControlSystem CurvedScalarWave CurvedScalarWaveHelpers DataStructures diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp index 0c22153666e38..d462f4e03486a 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp @@ -120,8 +120,8 @@ struct MockMetavariables { using const_global_cache_tags = tmpl::list< domain::Tags::Domain, CurvedScalarWave::Tags::BackgroundSpacetime, - Tags::ExcisionSphere, Tags::ExpansionOrder, Tags::MaxIterations, - Tags::Charge, Tags::Mass>; + Tags::ExcisionSphere, Tags::WorldtubeRadius, Tags::ExpansionOrder, + Tags::MaxIterations, Tags::Charge, Tags::Mass>; }; void test_iterations(const size_t max_iterations) { @@ -166,10 +166,15 @@ void test_iterations(const size_t max_iterations) { tuples::TaggedTuple< domain::Tags::Domain, CurvedScalarWave::Tags::BackgroundSpacetime, - Tags::ExcisionSphere, Tags::ExpansionOrder, Tags::MaxIterations, - Tags::Charge, Tags::Mass> - tuple_of_opts{shell.create_domain(), kerr_schild, excision_sphere, - expansion_order, max_iterations, charge, + Tags::ExcisionSphere, Tags::WorldtubeRadius, Tags::ExpansionOrder, + Tags::MaxIterations, Tags::Charge, Tags::Mass> + tuple_of_opts{shell.create_domain(), + kerr_schild, + excision_sphere, + excision_sphere.radius(), + expansion_order, + max_iterations, + charge, std::make_optional(mass)}; ActionTesting::MockRuntimeSystem runner{std::move(tuple_of_opts)}; const auto element_ids = initial_element_ids(initial_refinements); diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp index 4632547b1224a..28dc78bb6a6ff 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp @@ -122,8 +122,8 @@ struct MockMetavariables { using dg_element_array = MockElementArray; using const_global_cache_tags = tmpl::list, Tags::ExcisionSphere, - Tags::ExpansionOrder, Tags::MaxIterations, Tags::Charge, - Tags::Mass>; + Tags::WorldtubeRadius, Tags::ExpansionOrder, + Tags::MaxIterations, Tags::Charge, Tags::Mass>; }; // This test checks that `SendToWorldtube` integrates the regular field on the @@ -165,10 +165,11 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") { const auto& initial_extents = shell.initial_extents(); // self force and therefore iterative scheme is turned off tuples::TaggedTuple, Tags::ExcisionSphere, - Tags::ExpansionOrder, Tags::MaxIterations, Tags::Charge, - Tags::Mass> + Tags::WorldtubeRadius, Tags::ExpansionOrder, + Tags::MaxIterations, Tags::Charge, Tags::Mass> tuple_of_opts{shell.create_domain(), excision_sphere, + excision_sphere.radius(), expansion_order, static_cast(0), 0.1, diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_InitializeEvolvedVariables.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_InitializeEvolvedVariables.cpp index a8736cc446648..0b6b02f53c392 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_InitializeEvolvedVariables.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_InitializeEvolvedVariables.cpp @@ -20,24 +20,30 @@ SPECTRE_TEST_CASE( tmpl::list, Tags::EvolvedVelocity<3>>>; using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>; const tnsr::I initial_pos{{1., 2., 3.}}; + const tnsr::I initial_excision_pos{{1., 2., 3.}}; const tnsr::I initial_vel{{4., 5., 6.}}; const size_t current_iteration = 77; const double expiration_time = 1234.; const double initial_time = 0.; + const double initial_wt_radius = 12345.; + const double wt_radius = 2.; + const ExcisionSphere<3> excision_sphere(wt_radius, initial_excision_pos, {}); auto box = db::create, ::Tags::ConcreteTimeStepper, Tags::InitialPositionAndVelocity, Tags::CurrentIteration, - Tags::ExpirationTime, ::Tags::Time>, + Tags::ExpirationTime, Tags::WorldtubeRadius, ::Tags::Time, + Tags::ExcisionSphere<3>>, time_stepper_ref_tags>( variables_tag::type{}, dt_variables_tag::type{}, TimeSteppers::History{}, static_cast>( std::make_unique(4)), std::array, 2>{{initial_pos, initial_vel}}, - current_iteration, expiration_time, initial_time); + current_iteration, expiration_time, initial_wt_radius, initial_time, + excision_sphere); db::mutate_apply( make_not_null(&box)); @@ -52,6 +58,7 @@ SPECTRE_TEST_CASE( TimeSteppers::History(1)); CHECK(get(box) == 0); CHECK(get(box) == initial_time + 1e-10); + CHECK(get(box) == wt_radius); } } // namespace } // namespace CurvedScalarWave::Worldtube diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateFunctionsOfTime.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateFunctionsOfTime.cpp index 436fff6b5e3ae..e8b052ed4173a 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateFunctionsOfTime.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateFunctionsOfTime.cpp @@ -70,6 +70,10 @@ struct MockMetavariables { Tags::BlackHoleRadiusParameters>; }; +// this test takes a random position and velocity of the scalar charge and then +// calls `UpdateFunctionsOfTime`. After this action, the functions of time +// should be set in such a way that the worldtube excision sphere has the same +// position, velocity as the particle with the appropriate worldtube radius. void test_particle_motion() { static constexpr size_t Dim = 3; using metavars = MockMetavariables; @@ -88,6 +92,7 @@ void test_particle_motion() { const auto bh_excision_sphere = domain.excision_spheres().at("ExcisionSphereB"); + // these will be overwritten const std::array wt_radius_options{{0.1, 5., 1e-2, 2.3}}; const std::array bh_radius_options{{0.1, 5., 1e-2, 0.123}}; @@ -156,6 +161,9 @@ void test_particle_motion() { // includes worldtube shape map const auto& grid_to_inertial_map_wt_boundary = wt_neighbor_block.moving_mesh_grid_to_inertial_map(); + + // checks the position and velocity of the worldtube excision sphere are now + // mapped according to the position and velocity of the particle const auto mapped_particle = grid_to_inertial_map_particle.coords_frame_velocity_jacobians( wt_excision_sphere.center(), 0.5, fots); @@ -163,6 +171,8 @@ void test_particle_motion() { CHECK_ITERABLE_APPROX(std::get<3>(mapped_particle), particle_pos_vel.at(1)); const double orbit_radius = get(magnitude(particle_pos_vel[0])); + // check that the grid points on the worldtube boundary are mapped according + // to the prescribed worldtube radius auto wt_boundary_point = wt_excision_sphere.center(); wt_boundary_point[0] += wt_excision_sphere.radius(); const auto mapped_wt_boundary_data =