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/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/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/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