diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/CMakeLists.txt index d9a1fa4cefb7..533ae0ca4597 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/CMakeLists.txt +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/CMakeLists.txt @@ -20,4 +20,5 @@ spectre_target_headers( ReceiveElementData.hpp SendToElements.hpp UpdateAcceleration.hpp + UpdateFunctionsOfTime.hpp ) diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp index 03dea963b43e..71c36f79a55b 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp @@ -33,6 +33,7 @@ struct InitializeEvolvedVariables { using simple_tags = tmpl::list>; using return_tags = simple_tags; @@ -41,7 +42,8 @@ struct InitializeEvolvedVariables { using mutable_global_cache_tags = tmpl::list<>; using simple_tags_from_options = tmpl::list; using argument_tags = tmpl::list<::Tags::TimeStepper, - Tags::InitialPositionAndVelocity>; + Tags::InitialPositionAndVelocity, + ::Tags::Time, Tags::ExcisionSphere>; static void apply( const gsl::not_null, Tags::EvolvedVelocity>>*> @@ -51,11 +53,20 @@ struct InitializeEvolvedVariables { ::Tags::dt>>>*> 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 std::array, 2>& initial_pos_and_vel, + 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; *time_stepper_history = diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/ReceiveElementData.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/ReceiveElementData.hpp index 5c228a52cced..922158b881dd 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 new file mode 100644 index 000000000000..928a4a69a57b --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateFunctionsOfTime.hpp @@ -0,0 +1,153 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include + +#include "ControlSystem/UpdateFunctionOfTime.hpp" +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "Domain/Creators/Tags/Domain.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/Inboxes.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/RadiusFunctions.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/Worldtube.hpp" +#include "Parallel/AlgorithmExecution.hpp" +#include "Parallel/GlobalCache.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/Tags.hpp" +#include "ParallelAlgorithms/Initialization/MutateAssign.hpp" +#include "Time/Tags/TimeStepId.hpp" +#include "Time/TimeStepId.hpp" +#include "Time/TimeSteppers/TimeStepper.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace control_system { +// forward declare +struct UpdateMultipleFunctionsOfTime; +} // namespace control_system + +namespace CurvedScalarWave::Worldtube::Actions { + +/*! + * \brief Updates the functions of time according to the motion of the + * worldtube. + * + * \details We demand that the scalar charge is always in the center of the + * worldtube and therefore deform the grid to the track the particle's motion. + * In addition, the worldtube and black hole excision sphere radii are adjusted + * according to smooth_broken_power_law. + */ +struct UpdateFunctionsOfTime { + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, + tuples::TaggedTuple& /*inboxes*/, + Parallel::GlobalCache& cache, + const ArrayIndex& /*array_index*/, ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + const double current_expiration_time = db::get(box); + const double time = db::get<::Tags::Time>(box); + const auto& particle_pos_vel = + db::get>(box); + const double& x = get<0>(particle_pos_vel[0]); + const double& y = get<1>(particle_pos_vel[0]); + const double& xdot = get<0>(particle_pos_vel[1]); + const double& ydot = get<1>(particle_pos_vel[1]); + const double r = hypot(x, y); + const double radial_vel = (xdot * x + ydot * y) / r; + const auto& excision_sphere = db::get>(box); + const double grid_radius_particle = + get(magnitude(excision_sphere.center())); + + 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); + const auto& bh_radius_params = + db::get(box); + const double wt_radius_grid = excision_sphere.radius(); + const double wt_radius_inertial = + smooth_broken_power_law(r, wt_radius_params[0], wt_radius_params[1], + wt_radius_params[2], wt_radius_params[3]); + const double wt_radius_derivative = smooth_broken_power_law_derivative( + r, wt_radius_params[0], wt_radius_params[1], wt_radius_params[2], + wt_radius_params[3]); + const double wt_radius_time_derivative = wt_radius_derivative * radial_vel; + + const auto& bh_excision_sphere = + db::get>(box).excision_spheres().at( + "ExcisionSphereB"); + const double bh_excision_radius_grid = bh_excision_sphere.radius(); + const double bh_excision_radius_inertial = + smooth_broken_power_law(r, bh_radius_params[0], bh_radius_params[1], + bh_radius_params[2], bh_radius_params[3]); + const double bh_excision_radius_derivative = + smooth_broken_power_law_derivative( + r, bh_radius_params[0], bh_radius_params[1], bh_radius_params[2], + bh_radius_params[3]); + const double bh_excision_radius_time_derivative = + bh_excision_radius_derivative * radial_vel; + + 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. + 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])}; + + const DataVector size_b_update{ + sqrt_4_pi * (bh_excision_radius_grid - + 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])}; + + // 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 + // decreases but large enough to avoid floating point precision issues. + const double new_expiration_time = + time + + 0.01 * (db::get<::Tags::Next<::Tags::TimeStepId>>(box).substep_time() - + time); + + db::mutate( + [&new_expiration_time](const auto expiration_time) { + *expiration_time = new_expiration_time; + }, + make_not_null(&box)); + db::mutate( + [&wt_radius_inertial](const auto wt_radius) { + *wt_radius = wt_radius_inertial; + }, + make_not_null(&box)); + std::unordered_map> + all_updates{}; + all_updates["Rotation"] = + std::make_pair(angular_update, new_expiration_time); + all_updates["Expansion"] = + 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(size_b_update, new_expiration_time); + + Parallel::mutate<::domain::Tags::FunctionsOfTime, + control_system::UpdateMultipleFunctionsOfTime>( + cache, current_expiration_time, all_updates); + + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; +} // namespace CurvedScalarWave::Worldtube::Actions diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp index 9eb200f9785d..7db0ee87208b 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp @@ -430,6 +430,21 @@ struct CurrentIteration : db::SimpleTag { using type = size_t; }; +/*! + * \brief The current expiration time of the functions of time which are + * controlled by the worldtube singleton. + */ +struct ExpirationTime : db::SimpleTag { + using type = double; +}; + +/*! + * \brief The current worldtube radius held by the singleton. + */ +struct WorldtubeRadius : db::SimpleTag { + using type = double; +}; + /*! * \brief The initial position and velocity of the scalar charge in inertial * coordinates. diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 242493c661ec..3b87f4f0030e 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -20,6 +20,7 @@ set(LIBRARY_SOURCES SingletonActions/Test_InitializeEvolvedVariables.cpp SingletonActions/Test_ObserveWorldtubeSolution.cpp SingletonActions/Test_UpdateAcceleration.cpp + SingletonActions/Test_UpdateFunctionsOfTime.cpp ) add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") @@ -27,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 0c22153666e3..d462f4e03486 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 4632547b1224..28dc78bb6a6f 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 be04765a7379..0b6b02f53c39 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_InitializeEvolvedVariables.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_InitializeEvolvedVariables.cpp @@ -7,6 +7,7 @@ #include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp" #include "Framework/TestingFramework.hpp" #include "Time/Tags/HistoryEvolvedVariables.hpp" +#include "Time/Tags/Time.hpp" #include "Time/Tags/TimeStepper.hpp" #include "Time/TimeSteppers/AdamsBashforth.hpp" #include "Time/TimeSteppers/TimeStepper.hpp" @@ -19,21 +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::InitialPositionAndVelocity, Tags::CurrentIteration, + 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); + current_iteration, expiration_time, initial_wt_radius, initial_time, + excision_sphere); db::mutate_apply( make_not_null(&box)); @@ -47,6 +57,8 @@ SPECTRE_TEST_CASE( CHECK(db::get<::Tags::HistoryEvolvedVariables>(box) == 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 new file mode 100644 index 000000000000..e8b052ed4173 --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateFunctionsOfTime.cpp @@ -0,0 +1,215 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Block.hpp" +#include "Domain/Creators/Tags/Domain.hpp" +#include "Domain/Creators/Tags/FunctionsOfTime.hpp" +#include "Domain/Domain.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" +#include "Domain/FunctionsOfTime/Tags.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateFunctionsOfTime.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonChare.hpp" +#include "Framework/ActionTesting.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.hpp" +#include "Parallel/Phase.hpp" +#include "Parallel/PhaseDependentActionList.hpp" +#include "ParallelAlgorithms/Actions/FunctionsOfTimeAreReady.hpp" +#include "Time/Tags/TimeStepId.hpp" +#include "Time/Time.hpp" +#include "Time/TimeStepId.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace CurvedScalarWave::Worldtube { +namespace { + +template +struct MockWorldtubeSingleton { + using metavariables = Metavariables; + using mutable_global_cache_tags = + tmpl::list; + static constexpr size_t Dim = metavariables::volume_dim; + using chare_type = ActionTesting::MockSingletonChare; + using array_index = int; + using phase_dependent_action_list = tmpl::list< + Parallel::PhaseActions< + Parallel::Phase::Initialization, + tmpl::list, + ::Tags::Next<::Tags::TimeStepId>, + Tags::WorldtubeRadius, Tags::ExpirationTime>, + db::AddComputeTags<>>>>, + Parallel::PhaseActions< + Parallel::Phase::Testing, + tmpl::list>>>; + using component_being_mocked = WorldtubeSingleton; +}; + +template +struct MockMetavariables { + static constexpr size_t volume_dim = Dim; + using component_list = tmpl::list>; + using const_global_cache_tags = + tmpl::list, domain::Tags::Domain<3>, + Tags::WorldtubeRadiusParameters, + 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; + domain::creators::register_derived_with_charm(); + using worldtube_chare = MockWorldtubeSingleton; + const double initial_worldtube_radius = 0.456; + + const auto bco_domain_creator = + TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object< + true>(6., initial_worldtube_radius, 0.1); + domain::FunctionsOfTime::register_derived_with_charm(); + + const auto domain = bco_domain_creator->create_domain(); + const auto wt_excision_sphere = + domain.excision_spheres().at("ExcisionSphereA"); + 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}}; + + ActionTesting::MockRuntimeSystem runner{ + { + wt_excision_sphere, + bco_domain_creator->create_domain(), + wt_radius_options, + bh_radius_options, + }, + bco_domain_creator->functions_of_time()}; + + const double current_time = 1.; + const Time next_time{{current_time, 2.}, {1, 2}}; + const TimeStepId next_time_step_id{true, 123, next_time}; + std::array, 2> particle_pos_vel{}; + MAKE_GENERATOR(gen); + std::uniform_real_distribution pos_dist(-10., 10.); + std::uniform_real_distribution vel_dist(-0.2, 0.2); + particle_pos_vel.at(0) = + tnsr::I{{pos_dist(gen), pos_dist(gen), 0.}}; + particle_pos_vel.at(1) = + tnsr::I{{vel_dist(gen), vel_dist(gen), 0.}}; + const double initial_expiration_time = 1e-10; + ActionTesting::emplace_singleton_component_and_initialize( + &runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}, + {1., particle_pos_vel, next_time_step_id, initial_worldtube_radius, + initial_expiration_time}); + const auto& global_cache = ActionTesting::cache(runner, 0); + + const auto& fots = get(global_cache); + CHECK(fots.at("Rotation")->time_bounds()[1] == + approx(initial_expiration_time)); + CHECK(fots.at("Expansion")->time_bounds()[1] == + approx(initial_expiration_time)); + CHECK(fots.at("SizeA")->time_bounds()[1] == approx(initial_expiration_time)); + CHECK(fots.at("SizeB")->time_bounds()[1] == approx(initial_expiration_time)); + ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); + // UpdateFunctionsOfTime + CHECK(ActionTesting::next_action_if_ready( + make_not_null(&runner), 0)); + // CheckFunctionsOfTimeAreReady + CHECK(ActionTesting::next_action_if_ready( + make_not_null(&runner), 0)); + + const double new_expiration_time_expected = 1.005; + + CHECK(fots.at("Rotation")->time_bounds()[1] == + approx(new_expiration_time_expected)); + CHECK(fots.at("Expansion")->time_bounds()[1] == + approx(new_expiration_time_expected)); + CHECK(fots.at("SizeA")->time_bounds()[1] == + approx(new_expiration_time_expected)); + CHECK(fots.at("SizeB")->time_bounds()[1] == + approx(new_expiration_time_expected)); + + const auto& expiration_time_tag = + ActionTesting::get_databox_tag( + runner, 0); + CHECK(expiration_time_tag == approx(new_expiration_time_expected)); + + const auto& wt_neighbor_block = domain.blocks()[1]; + // does not include shape map + const auto& grid_to_inertial_map_particle = + wt_excision_sphere.moving_mesh_grid_to_inertial_map(); + // 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); + CHECK_ITERABLE_APPROX(std::get<0>(mapped_particle), particle_pos_vel.at(0)); + 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 = + grid_to_inertial_map_wt_boundary.coords_frame_velocity_jacobians( + wt_boundary_point, 0.5, fots); + const auto& mapped_wt_boundary_point = std::get<0>(mapped_wt_boundary_data); + const double wt_radius = get(magnitude(tenex::evaluate( + mapped_wt_boundary_point(ti::I) - particle_pos_vel.at(0)(ti::I)))); + const double wt_radius_expected = smooth_broken_power_law( + orbit_radius, wt_radius_options[0], wt_radius_options[1], + wt_radius_options[2], wt_radius_options[3]); + CHECK(wt_radius == approx(wt_radius_expected)); + + const auto& wt_radius_tag = + ActionTesting::get_databox_tag( + runner, 0); + CHECK(wt_radius_tag == approx(wt_radius_expected)); + + const auto& bh_neighbor_block = domain.blocks()[12]; + // includes bh shape map + const auto& grid_to_inertial_map_bh_boundary = + bh_neighbor_block.moving_mesh_grid_to_inertial_map(); + auto bh_boundary_point = bh_excision_sphere.center(); + bh_boundary_point[0] += bh_excision_sphere.radius(); + const auto mapped_bh_boundary_point = std::get<0>( + grid_to_inertial_map_bh_boundary.coords_frame_velocity_jacobians( + bh_boundary_point, 0.5, fots)); + const double bh_radius = get(magnitude(mapped_bh_boundary_point)); + const double bh_radius_expected = smooth_broken_power_law( + orbit_radius, bh_radius_options[0], bh_radius_options[1], + bh_radius_options[2], bh_radius_options[3]); + CHECK(bh_radius == approx(bh_radius_expected)); +} + +SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.UpdateFunctionsOfTime", + "[Unit]") { + test_particle_motion(); +} +} // namespace +} // namespace CurvedScalarWave::Worldtube diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp index 27ff1f216463..be35e78f033e 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp @@ -71,8 +71,8 @@ void test_initial_position_velocity_tag() { const double orbital_radius = 7.; const double angular_vel = 0.1; const auto domain_creator = - TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object( - orbital_radius, 0.2, angular_vel); + TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object< + false>(orbital_radius, 0.2, angular_vel); const tnsr::I expected_pos{{orbital_radius, 0., 0.}}; const tnsr::I expected_vel{{0., angular_vel * orbital_radius, 0.}}; CHECK(Tags::InitialPositionAndVelocity::create_from_options( @@ -192,8 +192,8 @@ void test_compute_face_coordinates_grid() { void test_compute_face_coordinates() { static constexpr size_t Dim = 3; const auto domain_creator = - TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object( - 7., 0.2, pow(7, -1.5)); + TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object< + false>(7., 0.2, pow(7, -1.5)); const double initial_time = 0.; auto domain = domain_creator->create_domain(); const auto excision_sphere = domain.excision_spheres().at("ExcisionSphereA"); @@ -321,8 +321,8 @@ void test_particle_position_velocity_compute() { const double orbit_radius = 9.; const double angular_velocity = 1. / (sqrt(orbit_radius) * orbit_radius); const auto domain_creator = - TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object( - orbit_radius, 0.2, angular_velocity); + TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object< + false>(orbit_radius, 0.2, angular_velocity); const double initial_time = 0.; auto domain = domain_creator->create_domain(); const auto excision_sphere = domain.excision_spheres().at("ExcisionSphereA"); @@ -632,15 +632,15 @@ void test_puncture_field() { void test_check_input_file() { const auto bbh_correct = - TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object( - 7., 0.2, pow(7., -1.5)); + TestHelpers::CurvedScalarWave::Worldtube::worldtube_binary_compact_object< + false>(7., 0.2, pow(7., -1.5)); const gr::Solutions::KerrSchild kerr_schild_correct( 1., make_array(0., 0., 0.), make_array(0., 0., 0.)); CHECK(Tags::CheckInputFile<3, gr::Solutions::KerrSchild>::create_from_options( bbh_correct, "ExcisionSphereA", kerr_schild_correct)); { const auto bbh_incorrect = TestHelpers::CurvedScalarWave::Worldtube:: - worldtube_binary_compact_object(7., 0.2, 1.); + worldtube_binary_compact_object(7., 0.2, 1.); CHECK_THROWS_WITH( (Tags::CheckInputFile<3, gr::Solutions::KerrSchild>:: create_from_options(bbh_incorrect, "ExcisionSphereA", @@ -756,6 +756,10 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.Worldtube.Tags", TestHelpers::db::test_simple_tag< CurvedScalarWave::Worldtube::Tags::SelfForceTurnOnInterval>( "SelfForceTurnOnInterval"); + TestHelpers::db::test_simple_tag< + CurvedScalarWave::Worldtube::Tags::ExpirationTime>("ExpirationTime"); + TestHelpers::db::test_simple_tag< + CurvedScalarWave::Worldtube::Tags::WorldtubeRadius>("WorldtubeRadius"); TestHelpers::db::test_simple_tag< CurvedScalarWave::Worldtube::Tags::MaxIterations>("MaxIterations"); TestHelpers::db::test_simple_tag< diff --git a/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp b/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp index 34a33a83ffa3..d496ab2fa529 100644 --- a/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp +++ b/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp @@ -14,6 +14,8 @@ #include "Helpers/Domain/BoundaryConditions/BoundaryCondition.hpp" namespace TestHelpers::CurvedScalarWave::Worldtube { + +template std::unique_ptr> worldtube_binary_compact_object( const double orbit_radius, const double worldtube_radius, const double angular_velocity) { @@ -70,21 +72,25 @@ std::unique_ptr> worldtube_binary_compact_object( " InitialAngularVelocity: [0.0, 0.0," + angular_velocity_stream.str() + "]\n" - " TranslationMap:\n" - " InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.], [0.0, 0.0, " - "0.0]]\n" + " TranslationMap: None\n" " ShapeMapA:\n" - " LMax: 8\n" + " LMax: 2\n" " InitialValues: Spherical\n" " SizeInitialValues: [0.0, 0.0, 0.0]\n" " TransitionEndsAtCube: false\n" " ShapeMapB:\n" - " LMax: 8\n" + " LMax: 2\n" " InitialValues: Spherical\n" " SizeInitialValues: [0.0, 0.0, 0.0]\n" " TransitionEndsAtCube: false\n"; return ::TestHelpers::test_option_tag<::domain::OptionTags::DomainCreator<3>, - Metavariables>( + Metavariables>( binary_compact_object_options); } + +template std::unique_ptr> +worldtube_binary_compact_object(double, double, double); + +template std::unique_ptr> +worldtube_binary_compact_object(double, double, double); } // namespace TestHelpers::CurvedScalarWave::Worldtube diff --git a/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.hpp b/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.hpp index 6d7841d06fcb..0852252b245a 100644 --- a/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.hpp +++ b/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.hpp @@ -15,6 +15,7 @@ #include "Utilities/TMPL.hpp" namespace TestHelpers::CurvedScalarWave::Worldtube { +template struct Metavariables { using system = TestHelpers::domain::BoundaryConditions::SystemWithBoundaryConditions<3>; @@ -22,9 +23,9 @@ struct Metavariables { : tt::ConformsTo { // we set `UseWorldtube` to `false` here so the functions of time are valid // which simplifies testing. - using factory_classes = tmpl::map< - tmpl::pair, - tmpl::list<::domain::creators::BinaryCompactObject>>>; + using factory_classes = tmpl::map, + tmpl::list<::domain::creators::BinaryCompactObject>>>; }; }; @@ -32,6 +33,7 @@ struct Metavariables { // corresponding to the circular worldtube setup: a central excision sphere // (the central black hole) and a worldtube excision sphere in circular orbit // around it with angular velocity R^{-3/2}, where R is the orbital radius. +template std::unique_ptr> worldtube_binary_compact_object( const double orbit_radius, const double worldtube_radius, const double angular_velocity);