diff --git a/src/ControlSystem/Measurements/BNSCenterOfMass.hpp b/src/ControlSystem/Measurements/BNSCenterOfMass.hpp index f88c1a6e5a04..b287133282bd 100644 --- a/src/ControlSystem/Measurements/BNSCenterOfMass.hpp +++ b/src/ControlSystem/Measurements/BNSCenterOfMass.hpp @@ -14,6 +14,9 @@ #include "ControlSystem/RunCallbacks.hpp" #include "DataStructures/LinkedMessageId.hpp" #include "DataStructures/Tensor/TypeAliases.hpp" +#include "Domain/BlockLogicalCoordinates.hpp" +#include "Domain/Creators/Tags/Domain.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/Tags.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "IO/Logging/Verbosity.hpp" @@ -26,6 +29,7 @@ #include "ParallelAlgorithms/Actions/FunctionsOfTimeAreReady.hpp" #include "ParallelAlgorithms/Events/Tags.hpp" #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp" +#include "Utilities/Algorithm.hpp" #include "Utilities/GetOutput.hpp" #include "Utilities/ProtocolHelpers.hpp" #include "Utilities/Serialization/CharmPupable.hpp" @@ -273,18 +277,69 @@ struct PostReductionSendBNSStarCentersToControlSystem { center_databox, cache, measurement_id); if (Parallel::get(cache)) { - std::vector data_to_write{ + std::vector grid_data_to_write{ measurement_id.id, center_a[0], center_a[1], center_a[2], center_b[0], center_b[1], center_b[2]}; + // To convert grid coords to inertial coords, we must find the block that + // these coords are in and use that grid to inertial map + const Domain<3>& domain = Parallel::get>(cache); + const domain::FunctionsOfTimeMap& functions_of_time = + Parallel::get(cache); + tnsr::I grid_tnsr_center{}; + for (size_t i = 0; i < 3; i++) { + grid_tnsr_center.get(i) = + DataVector{gsl::at(center_a, i), gsl::at(center_b, i)}; + } + + const auto block_logical_coords = block_logical_coordinates( + domain, grid_tnsr_center, measurement_id.id, functions_of_time); + + ASSERT(alg::all_of(block_logical_coords, + [](const auto& logical_coord_holder) { + return logical_coord_holder.has_value(); + }), + "Grid centers of BNS (" + << center_a << ", " << center_b + << ") could not be mapped to the logical frame."); + + const auto& blocks = domain.blocks(); + std::vector inertial_data_to_write{measurement_id.id}; + + ASSERT(block_logical_coords.size() == 2, + "There should be exactly 2 block logical coordinates for the two " + "centers of the BNS, but instead there are " + << block_logical_coords.size()); + + for (size_t n = 0; n < block_logical_coords.size(); n++) { + const auto& logical_coord_holder = block_logical_coords[n]; + const auto& block_id = logical_coord_holder.value().id; + + const auto& block = blocks[block_id.get_index()]; + const auto& grid_to_inertial_map = + block.moving_mesh_grid_to_inertial_map(); + + const auto inertial_center = grid_to_inertial_map( + tnsr::I{n == 0 ? center_a : center_b}, + measurement_id.id, functions_of_time); + for (size_t i = 0; i < 3; i++) { + inertial_data_to_write.push_back(inertial_center.get(i)); + } + } + auto& writer_proxy = Parallel::get_parallel_component< observers::ObserverWriter>(cache); Parallel::threaded_action< observers::ThreadedActions::WriteReductionDataRow>( // Node 0 is always the writer - writer_proxy[0], subfile_path_, legend_, - std::make_tuple(std::move(data_to_write))); + writer_proxy[0], grid_subfile_path_, legend_, + std::make_tuple(std::move(grid_data_to_write))); + Parallel::threaded_action< + observers::ThreadedActions::WriteReductionDataRow>( + // Node 0 is always the writer + writer_proxy[0], inertial_subfile_path_, legend_, + std::make_tuple(std::move(inertial_data_to_write))); } } @@ -292,7 +347,10 @@ struct PostReductionSendBNSStarCentersToControlSystem { const static inline std::vector legend_{ "Time", "Center_A_x", "Center_A_y", "Center_A_z", "Center_B_x", "Center_B_y", "Center_B_z"}; - const static inline std::string subfile_path_{"/ControlSystems/BnsCenters"}; + const static inline std::string grid_subfile_path_{ + "/ControlSystems/BnsGridCenters"}; + const static inline std::string inertial_subfile_path_{ + "/ControlSystems/BnsInertialCenters"}; }; } // namespace measurements } // namespace control_system diff --git a/tests/Unit/ControlSystem/Test_Measurements.cpp b/tests/Unit/ControlSystem/Test_Measurements.cpp index f8e7885e6cf8..679c7966adfd 100644 --- a/tests/Unit/ControlSystem/Test_Measurements.cpp +++ b/tests/Unit/ControlSystem/Test_Measurements.cpp @@ -12,6 +12,13 @@ #include "ControlSystem/Measurements/SingleHorizon.hpp" #include "ControlSystem/Tags/SystemTags.hpp" #include "DataStructures/DataVector.hpp" +#include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/Creators/BinaryCompactObject.hpp" +#include "Domain/Creators/RegisterDerivedWithCharm.hpp" +#include "Domain/Creators/Tags/Domain.hpp" +#include "Domain/Creators/Tags/FunctionsOfTime.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" +#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Framework/ActionTesting.hpp" #include "Helpers/ControlSystem/Examples.hpp" @@ -24,6 +31,7 @@ #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp" #include "ParallelAlgorithms/Interpolation/Protocols/InterpolationTargetTag.hpp" #include "Time/Tags/Time.hpp" +#include "Utilities/Gsl.hpp" #include "Utilities/ProtocolHelpers.hpp" #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" #include "Utilities/TMPL.hpp" @@ -139,15 +147,32 @@ void test_serialization() { } // Check that the center of mass is at the expected location for this test. -void check_centers(const double center_a_x, const double center_a_y, - const double center_a_z, const double center_b_x, - const double center_b_y, const double center_b_z) { - CHECK(center_a_x == 1.0); - CHECK(center_a_y == 4.5 / 8.0); - CHECK(center_a_z == 3.0 / 8.0); - CHECK(center_b_x == -1.0); - CHECK(center_b_y == 1.0); - CHECK(center_b_z == 0.0); +void check_centers(const std::array& center_a, + const std::array& center_b, + const bool grid_centers = true) { + const std::array expected_grid_center_a{1.0, 4.5 / 8.0, 3.0 / 8.0}; + const std::array expected_grid_center_b{-1.0, 1.0, 0.0}; + + if (grid_centers) { + CHECK(center_a == expected_grid_center_a); + CHECK(center_b == expected_grid_center_b); + } else { + const std::array expected_inertial_center_a{ + cos(0.1) * expected_grid_center_a[0] + + -sin(0.1) * expected_grid_center_a[1], + sin(0.1) * expected_grid_center_a[0] + + cos(0.1) * expected_grid_center_a[1], + expected_grid_center_a[2]}; + const std::array expected_inertial_center_b{ + cos(0.1) * expected_grid_center_b[0] + + -sin(0.1) * expected_grid_center_b[1], + sin(0.1) * expected_grid_center_b[0] + + cos(0.1) * expected_grid_center_b[1], + expected_grid_center_b[2]}; + + CHECK_ITERABLE_APPROX(center_a, expected_inertial_center_a); + CHECK_ITERABLE_APPROX(center_b, expected_inertial_center_b); + } } struct MockControlSystem @@ -181,8 +206,7 @@ struct MockControlSystem const std::array center_b, Parallel::GlobalCache& /*cache*/, const LinkedMessageId& /*measurement_id*/) { - check_centers(center_a[0], center_a[1], center_a[2], center_b[0], - center_b[1], center_b[2]); + check_centers(center_a, center_b); // Avoid unused variable warning for deriv_order, which is required // as part of the control_system protocol. CHECK(2 == deriv_order); @@ -194,7 +218,10 @@ template struct MockControlSystemComponent { using component_being_mocked = ControlComponent; - using const_global_cache_tags = tmpl::list; + using const_global_cache_tags = + tmpl::list>; + using mutable_global_cache_tags = + tmpl::list; using metavariables = Metavariables; using chare_type = ActionTesting::MockSingletonChare; using array_index = int; @@ -213,7 +240,8 @@ struct MockWriteReductionDataRow { const std::string& subfile_name, const std::vector& file_legend, const std::tuple>& data_row) { - CHECK(subfile_name == "/ControlSystems/BnsCenters.dat"); + CHECK((subfile_name == "/ControlSystems/BnsGridCenters" or + subfile_name == "/ControlSystems/BnsInertialCenters")); CHECK(file_legend == std::vector{ "Time", "Center_A_x", "Center_A_y", "Center_A_z", "Center_B_x", "Center_B_y", "Center_B_z"}); @@ -222,7 +250,8 @@ struct MockWriteReductionDataRow { CHECK(data[0] == 1.0); CHECK(data.size() == file_legend.size()); - check_centers(data[1], data[2], data[3], data[4], data[5], data[6]); + check_centers({data[1], data[2], data[3]}, {data[4], data[5], data[6]}, + subfile_name == "/ControlSystems/BnsGridCenters"); } }; @@ -257,6 +286,8 @@ struct Metavariables { SPECTRE_TEST_CASE("Unit.ControlSystem.FindTwoCenters", "[ControlSystem][Unit]") { test_serialization(); + domain::FunctionsOfTime::register_derived_with_charm(); + domain::creators::register_derived_with_charm(); // Part 1 of the test: calculation of the relevant integrals // within a (mock) element. @@ -287,11 +318,39 @@ SPECTRE_TEST_CASE("Unit.ControlSystem.FindTwoCenters", CHECK(first_moment_a == std::array{8.0, 4.5, 3.0}); CHECK(first_moment_b == std::array{-2.5, 2.5, 0.0}); + // Use BNS domain with rotation for realistic block structure, but simple time + // dependent maps for testing + const domain::creators::bco::TimeDependentMapOptions time_dep_opts{ + 0.0, + std::nullopt, + domain::creators::bco::TimeDependentMapOptions::RotationMapOptions{ + std::array{0.0, 0.0, 0.1}}, + std::nullopt, + std::nullopt, + std::nullopt}; + const domain::creators::BinaryCompactObject binary_compact_object{ + domain::creators::BinaryCompactObject::CartesianCubeAtXCoord{20.0}, + domain::creators::BinaryCompactObject::CartesianCubeAtXCoord{ + -20.0}, + {0.0, 0.0}, + 100.0, + 500.0, + 1.0, + {0_st}, + {4_st}, + true, + domain::CoordinateMaps::Distribution::Projective, + domain::CoordinateMaps::Distribution::Linear, + 90.0, + time_dep_opts}; + using MockRuntimeSystem = ActionTesting::MockRuntimeSystem; using control_system_component = MockControlSystemComponent; using obs_writer = MockObserverWriter; - MockRuntimeSystem runner{{true, ::Verbosity::Silent}}; + MockRuntimeSystem runner{ + {true, ::Verbosity::Silent, binary_compact_object.create_domain()}, + {binary_compact_object.functions_of_time()}}; ActionTesting::emplace_singleton_component( make_not_null(&runner), ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}); @@ -312,4 +371,12 @@ SPECTRE_TEST_CASE("Unit.ControlSystem.FindTwoCenters", template apply(box, cache, unused_element_id, measurement_id, mass_a, mass_b, first_moment_a, first_moment_b); + + // One for grid centers, one for inertial centers + CHECK(ActionTesting::number_of_queued_threaded_actions(runner, + 0) == 2); + for (size_t i = 0; i < 2; i++) { + ActionTesting::invoke_queued_threaded_action( + make_not_null(&runner), 0); + } }