From 8e243dff8dfe33ff72c571dcc39a6dacdb8d2983 Mon Sep 17 00:00:00 2001 From: Francois Foucart Date: Mon, 1 Jul 2024 15:49:38 -0400 Subject: [PATCH] Infrastructure for MC communication --- .../Particles/MonteCarlo/CMakeLists.txt | 5 + .../MonteCarlo/GhostZoneCommunication.hpp | 391 +++++++++++ .../MonteCarlo/GhostZoneCommunicationStep.hpp | 30 + .../MonteCarlo/GhostZoneCommunicationTags.hpp | 96 +++ .../Particles/MonteCarlo/MortarData.hpp | 47 ++ src/Evolution/Particles/MonteCarlo/Packet.cpp | 10 + src/Evolution/Particles/MonteCarlo/Packet.hpp | 43 +- src/Evolution/Particles/MonteCarlo/Tags.hpp | 31 + .../Particles/MonteCarlo/TakeTimeStep.tpp | 2 +- .../Particles/MonteCarlo/CMakeLists.txt | 3 + .../MonteCarlo/Test_CommunicationTags.cpp | 638 ++++++++++++++++++ 11 files changed, 1290 insertions(+), 6 deletions(-) create mode 100644 src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/MortarData.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/Tags.hpp create mode 100644 tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp diff --git a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt index 14bcb54c018b..0ac636bf619c 100644 --- a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt +++ b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt @@ -29,11 +29,16 @@ spectre_target_headers( EmitPackets.tpp EvolvePackets.hpp EvolvePacketsInElement.tpp + GhostZoneCommunication.hpp + GhostZoneCommunicationStep.hpp + GhostZoneCommunicationTags.hpp ImplicitMonteCarloCorrections.tpp InverseJacobianInertialToFluidCompute.hpp + MortarData.hpp NeutrinoInteractionTable.hpp Packet.hpp Scattering.hpp + Tags.hpp TakeTimeStep.tpp TemplatedLocalFunctions.hpp ) diff --git a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp new file mode 100644 index 000000000000..2f9422a5848a --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp @@ -0,0 +1,391 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Index.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalId.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Structure/OrientationMapHelpers.hpp" +#include "Domain/Tags.hpp" +#include "Domain/Tags/NeighborMesh.hpp" +#include "Evolution/DgSubcell/ActiveGrid.hpp" +#include "Evolution/DgSubcell/GhostData.hpp" +#include "Evolution/DgSubcell/SliceData.hpp" +#include "Evolution/DgSubcell/Tags/ActiveGrid.hpp" +#include "Evolution/DgSubcell/Tags/Interpolators.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp" +#include "Evolution/Particles/MonteCarlo/MortarData.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "Parallel/AlgorithmExecution.hpp" +#include "Parallel/GlobalCache.hpp" +#include "PointwiseFunctions/Hydro/TagsDeclarations.hpp" +#include "Time/TimeStepId.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Gsl.hpp" + +/// \cond +namespace Tags { +struct TimeStepId; +} // namespace Tags +/// \endcond + +namespace Particles::MonteCarlo::Actions { + +/// Mutator to get required volume data for communication +/// before a MC step; i.e. data sent from live points +/// to ghost points in neighbors. +struct GhostDataMutatorPreStep { + using return_tags = tmpl::list<>; + using argument_tags = tmpl::list< + hydro::Tags::RestMassDensity, + hydro::Tags::ElectronFraction, + hydro::Tags::Temperature, + Particles::MonteCarlo::Tags::CellLightCrossingTime>; + static const size_t number_of_vars = 4; + + static DataVector apply(const Scalar& rest_mass_density, + const Scalar& electron_fraction, + const Scalar& temperature, + const Scalar& cell_light_crossing_time) { + const size_t dv_size = get(rest_mass_density).size(); + DataVector buffer{dv_size * number_of_vars}; + std::copy( + get(rest_mass_density).data(), + std::next(get(rest_mass_density).data(), static_cast(dv_size)), + buffer.data()); + std::copy( + get(electron_fraction).data(), + std::next(get(electron_fraction).data(), static_cast(dv_size)), + std::next(buffer.data(), static_cast(dv_size))); + std::copy(get(temperature).data(), + std::next(get(temperature).data(), static_cast(dv_size)), + std::next(buffer.data(), static_cast(dv_size * 2))); + std::copy(get(cell_light_crossing_time).data(), + std::next(get(cell_light_crossing_time).data(), + static_cast(dv_size)), + std::next(buffer.data(), static_cast(dv_size * 3))); + return buffer; + } +}; + +/// Mutator that returns packets currently in ghost zones in a +/// DirectionMap> +/// and remove them from the list of packets of the current +/// element +template +struct GhostDataMcPackets { + using return_tags = tmpl::list; + using argument_tags = tmpl::list<::domain::Tags::Element>; + + static DirectionMap> apply( + const gsl::not_null*> packets, + const Element& element) { + DirectionMap> output{}; + for (const auto& [direction, neighbors_in_direction] : + element.neighbors()) { + output[direction] = std::vector{}; + } + // Loop over packets. If a packet is out of the current element, + // remove it from the element and copy it to the appropriate + // neighbor list if it exists. + // Note: At this point, we only handle simple boundaries + // (one neighbor per direction, with the logical coordinates + // being transformed by adding/removing 2). We assume that + // packets that go out of the element but do not match any + // neighbor have reached the domain boundary and can be removed. + size_t n_packets = packets->size(); + for (size_t p = 0; p < n_packets; p++) { + Packet& packet = (*packets)[p]; + std::optional> max_distance_direction = std::nullopt; + double max_distance = 0.0; + // TO DO: Deal with dimensions higher than the grid dimension + // e.g. for axisymmetric evolution + for (size_t d = 0; d < Dim; d++) { + if (packet.coordinates.get(d) < -1.0) { + const double distance = -1.0 - packet.coordinates.get(d); + if (max_distance < distance) { + max_distance = distance; + max_distance_direction = Direction(d, Side::Lower); + } + } + if (packet.coordinates.get(d) > 1.0) { + const double distance = packet.coordinates.get(d) - 1.0; + if (max_distance < distance) { + max_distance = distance; + max_distance_direction = Direction(d, Side::Upper); + } + } + } + // If a packet should be moved along an edge / corner, we move + // it into the direct neighbor it is effectively closest to + // (in topological coordinates). This may need improvements. + if (max_distance_direction != std::nullopt) { + const size_t d = max_distance_direction.value().dimension(); + const Side side = max_distance_direction.value().side(); + packet.coordinates.get(d) += (side == Side::Lower) ? (2.0) : (-2.0); + if (output.contains(max_distance_direction.value())) { + output[max_distance_direction.value()].push_back(packet); + } + std::swap((*packets)[p], (*packets)[n_packets - 1]); + packets->pop_back(); + p--; + n_packets--; + } + } + return output; + } +}; + +/// Action responsible for the Send operation of ghost +/// zone communication for Monte-Carlo transport. +/// If CommStep == PreStep, this sends the fluid and +/// metric variables needed for evolution of MC packets. +/// If CommStep == PostStep, this sends packets that +/// have moved from one element to another. +/// The current code does not yet deal with data required +/// to calculate the backreaction on the fluid. +template +struct SendDataForMcCommunication { + using inbox_tags = + tmpl::list>; + + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, tuples::TaggedTuple& /*inboxes*/, + Parallel::GlobalCache& cache, + const ArrayIndex& /*array_index*/, const ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + static_assert(not LocalTimeStepping, + "Monte Carlo is not tested with local time stepping."); + + ASSERT(db::get(box) == + evolution::dg::subcell::ActiveGrid::Subcell, + "The SendDataForReconstructionPreStep action in MC " + "can only be called when Subcell is the active scheme."); + db::mutate>( + [](const auto ghost_data_ptr) { + // Clear the previous neighbor data and add current local data + ghost_data_ptr->clear(); + }, + make_not_null(&box)); + + auto& receiver_proxy = + Parallel::get_parallel_component(cache); + const size_t ghost_zone_size = 1; + + const Element& element = db::get<::domain::Tags::Element>(box); + const Mesh& subcell_mesh = + db::get>(box); + const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box); + + // Fill volume data that should be fed to neighbor ghost zones + DataVector volume_data_to_slice = + CommStep == Particles::MonteCarlo::CommunicationStep::PreStep + ? db::mutate_apply(GhostDataMutatorPreStep{}, make_not_null(&box)) + : DataVector{}; + const DirectionMap all_sliced_data = + CommStep == Particles::MonteCarlo::CommunicationStep::PreStep + ? evolution::dg::subcell::slice_data( + volume_data_to_slice, subcell_mesh.extents(), ghost_zone_size, + element.internal_boundaries(), 0, + db::get>(box)) + : DirectionMap{}; + + const DirectionMap> + all_packets_ghost_zone = + CommStep == Particles::MonteCarlo::CommunicationStep::PostStep + ? db::mutate_apply(GhostDataMcPackets{}, + make_not_null(&box)) + : DirectionMap>{}; + + // TO DO: + // Need to deal with coupling data, which is brought back from + // neighbor's ghost zones to the processor with the live cell! + + for (const auto& [direction, neighbors_in_direction] : + element.neighbors()) { + const auto& orientation = neighbors_in_direction.orientation(); + const auto direction_from_neighbor = orientation(direction.opposite()); + for (const ElementId& neighbor : neighbors_in_direction) { + std::optional> + packets_to_send = std::nullopt; + DataVector subcell_data_to_send{}; + + if (CommStep == Particles::MonteCarlo::CommunicationStep::PreStep) { + const auto& sliced_data_in_direction = all_sliced_data.at(direction); + subcell_data_to_send = DataVector{sliced_data_in_direction.size()}; + std::copy(sliced_data_in_direction.begin(), + sliced_data_in_direction.end(), + subcell_data_to_send.begin()); + } + if (CommStep == Particles::MonteCarlo::CommunicationStep::PostStep) { + packets_to_send = all_packets_ghost_zone.at(direction); + if (packets_to_send.value().empty()) { + packets_to_send = std::nullopt; + } + } + + McGhostZoneData data{subcell_data_to_send, packets_to_send}; + + Parallel::receive_data< + Particles::MonteCarlo::McGhostZoneDataInboxTag>( + receiver_proxy[neighbor], time_step_id, + std::pair{DirectionalId{direction_from_neighbor, element.id()}, + std::move(data)}); + } + } + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; + +/// Action responsible for the Receive operation of ghost +/// zone communication for Monte-Carlo transport. +/// If CommStep == PreStep, this gets the fluid and +/// metric variables needed for evolution of MC packets. +/// If CommStep == PostStep, this gets packets that +/// have moved into the current element. +/// The current code does not yet deal with data required +/// to calculate the backreaction on the fluid. +template +struct ReceiveDataForMcCommunication { + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, tuples::TaggedTuple& inboxes, + const Parallel::GlobalCache& /*cache*/, + const ArrayIndex& /*array_index*/, const ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + const Element& element = db::get<::domain::Tags::Element>(box); + const auto number_of_expected_messages = element.neighbors().size(); + if (UNLIKELY(number_of_expected_messages == 0)) { + // We have no neighbors, so just continue without doing any work + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } + + const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box); + std::map>>& + inbox = tuples::get>(inboxes); + const auto& received = inbox.find(current_time_step_id); + // Check we have at least some data from correct time, and then check that + // we have received all data + if (received == inbox.end() or + received->second.size() != number_of_expected_messages) { + return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + } + + // Now that we have received all the data, copy it over as needed. + DirectionalIdMap> + received_data = std::move(inbox[current_time_step_id]); + inbox.erase(current_time_step_id); + + // Copy the received PreStep data in the MortarData container + if (CommStep == Particles::MonteCarlo::CommunicationStep::PreStep) { + db::mutate>( + [&element, + &received_data](const gsl::not_null*> mortar_data) { + for (const auto& [direction, neighbors_in_direction] : + element.neighbors()) { + for (const auto& neighbor : neighbors_in_direction) { + DirectionalId directional_element_id{direction, neighbor}; + const DataVector& received_data_direction = + received_data[directional_element_id] + .ghost_zone_hydro_variables; + REQUIRE(received_data[directional_element_id] + .packets_entering_this_element == std::nullopt); + if (mortar_data->rest_mass_density[directional_element_id] == + std::nullopt) { + continue; + } + const size_t mortar_data_size = + mortar_data->rest_mass_density[directional_element_id] + .value() + .size(); + std::copy(received_data_direction.data(), + std::next(received_data_direction.data(), + static_cast(mortar_data_size)), + mortar_data->rest_mass_density[directional_element_id] + .value() + .data()); + std::copy(std::next(received_data_direction.data(), + static_cast(mortar_data_size)), + std::next(received_data_direction.data(), + static_cast(mortar_data_size * 2)), + mortar_data->electron_fraction[directional_element_id] + .value() + .data()); + std::copy(std::next(received_data_direction.data(), + static_cast(mortar_data_size * 2)), + std::next(received_data_direction.data(), + static_cast(mortar_data_size * 3)), + mortar_data->temperature[directional_element_id] + .value() + .data()); + std::copy(std::next(received_data_direction.data(), + static_cast(mortar_data_size * 3)), + std::next(received_data_direction.data(), + static_cast(mortar_data_size * 4)), + mortar_data + ->cell_light_crossing_time[directional_element_id] + .value() + .data()); + } + } + }, + make_not_null(&box)); + } else { + // TO DO: Deal with data coupling neutrinos back to fluid evolution + db::mutate( + [&element, &received_data]( + const gsl::not_null*> + packet_list) { + for (const auto& [direction, neighbors_in_direction] : + element.neighbors()) { + for (const auto& neighbor : neighbors_in_direction) { + DirectionalId directional_element_id{direction, neighbor}; + const std::optional>& + received_data_packets = + received_data[directional_element_id] + .packets_entering_this_element; + // Temporary: currently no data for coupling to the fluid + REQUIRE(received_data[directional_element_id] + .ghost_zone_hydro_variables.size() == 0); + if (received_data_packets == std::nullopt) { + continue; + } else { + const size_t n_packets = received_data_packets.value().size(); + for (size_t p = 0; p < n_packets; p++) { + packet_list->push_back(received_data_packets.value()[p]); + } + } + } + } + }, + make_not_null(&box)); + } + + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; + +} // namespace Particles::MonteCarlo::Actions diff --git a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp new file mode 100644 index 000000000000..ee794df9cc2e --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp @@ -0,0 +1,30 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +/// Items related to the evolution of particles +/// Items related to Monte-Carlo radiation transport +namespace Particles::MonteCarlo { + +/// This class is used to template communication actions +/// in the Monte-Carlo code, to know which of the two +/// communication steps we need to take care of +/// (just before or just after the MC step). +enum class CommunicationStep { + /// PreStep should occur just before the MC step. + /// It should send to the ghost zones of each + /// element the fluid and metric data. + PreStep, + /// PostStep should occur just after a MC step. + /// It sends packets that have left their current + /// element to the element they now belong to. + /// NOT CODE YET: This is also when information + /// about energy/momentum/lepton number exchance + /// in the ghost zones should be communicated back + /// to the corresponding live point for coupling to + /// the fluid evolution. + PostStep +}; + +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp new file mode 100644 index 000000000000..576a3861e1db --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp @@ -0,0 +1,96 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include + +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp" +#include "Evolution/Particles/MonteCarlo/Packet.hpp" +#include "Time/Slab.hpp" +#include "Time/Tags/TimeStepId.hpp" +#include "Time/Time.hpp" +#include "Time/TimeStepId.hpp" +#include "Utilities/Gsl.hpp" + +/// \cond +class DataVector; +/// \endcond + +/// Items related to the evolution of particles +/// Items related to Monte-Carlo radiation transport +namespace Particles::MonteCarlo { + +/*! + * \brief The container for ghost zone data to be communicated + * in the Monte-Carlo algorithm + * + * The stored data consists of + * (1) The fluid variables in those ghost cells, in the form + * of a DataVector. The data should contain the + * scalars (baryon density, temperature, electron fraction, + * and cell light-crossing time) in pre-step communication + * and (coupling_tilde_tau, coupling_tilde_s, coupling_rho_ye) + * in post-step communication + * (2) The packets moved from a neighbor to this element. This + * can be std::null_ptr in pre-step communication + */ +template +struct McGhostZoneData { + void pup(PUP::er& p) { + p | ghost_zone_hydro_variables; + p | packets_entering_this_element; + } + + DataVector ghost_zone_hydro_variables{}; + std::optional> + packets_entering_this_element{}; +}; + +/// Inbox tag to be used before MC step +template +struct McGhostZoneDataInboxTag { + using stored_type = McGhostZoneData; + + using temporal_id = TimeStepId; + using type = std::map>; + using value_type = type; + + CommunicationStep comm_step = CommStep; + + template + static size_t insert_into_inbox(const gsl::not_null inbox, + const temporal_id& time_step_id, + ReceiveDataType&& data) { + auto& current_inbox = (*inbox)[time_step_id]; + if (not current_inbox.insert(std::forward(data)) + .second) { + ERROR("Failed to insert data to receive at instance '" + << time_step_id + << " in McGhostZonePreStepDataInboxTag.\n"); + } + return current_inbox.size(); + } + + void pup(PUP::er& /*p*/) {} +}; + +// Tags to put in DataBox for MC communication. +namespace Tags { + +/// Simple tag for the structure containing ghost zone data +/// for Monte Carlo radiation transport. +template +struct McGhostZoneDataTag : db::SimpleTag { + using type = DirectionalIdMap>; +}; + +} // namespace Tags + +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/MortarData.hpp b/src/Evolution/Particles/MonteCarlo/MortarData.hpp new file mode 100644 index 000000000000..e13066d61057 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/MortarData.hpp @@ -0,0 +1,47 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/DataBox/Tag.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" + +/// \cond +class DataVector; +/// \endcond + +namespace Particles::MonteCarlo { + +/// Structure used to gather ghost zone data for Monte-Carlo evolution. +/// We need the rest mass density, electron fraction, temperature, and +/// an estimate of the light-crossing time one cell deep within each +/// neighboring element. +template +struct MortarData { + DirectionalIdMap> rest_mass_density{}; + DirectionalIdMap> electron_fraction{}; + DirectionalIdMap> temperature{}; + DirectionalIdMap> cell_light_crossing_time{}; + + void pup(PUP::er& p) { + p | rest_mass_density; + p | electron_fraction; + p | temperature; + p | cell_light_crossing_time; + } +}; + +namespace Tags { + +/// Simple tag containing the fluid and metric data in the ghost zones +/// for Monte-Carlo packets evolution. +template +struct MortarDataTag : db::SimpleTag { + using type = MortarData; +}; + +} // namespace Tags + +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/Packet.cpp b/src/Evolution/Particles/MonteCarlo/Packet.cpp index ce664c1dae0f..b2b3de5ae70d 100644 --- a/src/Evolution/Particles/MonteCarlo/Packet.cpp +++ b/src/Evolution/Particles/MonteCarlo/Packet.cpp @@ -5,6 +5,16 @@ namespace Particles::MonteCarlo { +void Packet::pup(PUP::er& p) { + p | species; + p | number_of_neutrinos; + p | index_of_closest_grid_point; + p | time; + p | momentum_upper_t; + p | coordinates; + p | momentum; +} + void Packet::renormalize_momentum( const tnsr::II& inv_spatial_metric, const Scalar& lapse) { diff --git a/src/Evolution/Particles/MonteCarlo/Packet.hpp b/src/Evolution/Particles/MonteCarlo/Packet.hpp index 450c7a327229..c4b911a5baee 100644 --- a/src/Evolution/Particles/MonteCarlo/Packet.hpp +++ b/src/Evolution/Particles/MonteCarlo/Packet.hpp @@ -35,25 +35,42 @@ struct Packet { momentum[2] = p_z_; } + /// Default constructor needed to make Packet puppable + /// For the same reason, we have default initialization + /// for a number of member variables to unphysical values. + Packet() + : species{0}, + number_of_neutrinos{0}, + index_of_closest_grid_point{0}, + time{-1.0}, + momentum_upper_t{0.0} { + coordinates[0] = 0.0; + coordinates[1] = 0.0; + coordinates[2] = 0.0; + momentum[0] = 0.0; + momentum[1] = 0.0; + momentum[2] = 0.0; + } + /// Species of neutrinos (in the code, just an index used to access the /// right interaction rates; typically \f$0=\nu_e, 1=\nu_a, 2=\nu_x\f$) - size_t species; + size_t species = std::numeric_limits::max(); /// Number of neutrinos represented by current packet /// Note that this number is rescaled so that /// `Energy_of_packet = N * Energy_of_neutrinos` /// with the packet energy in G=Msun=c=1 units but /// the neutrino energy in MeV! - double number_of_neutrinos; + double number_of_neutrinos = std::numeric_limits::signaling_NaN(); /// Index of the closest point on the FD grid. - size_t index_of_closest_grid_point; + size_t index_of_closest_grid_point = std::numeric_limits::max(); /// Current time - double time; + double time = std::numeric_limits::signaling_NaN(); /// Stores \f$p^t\f$ - double momentum_upper_t; + double momentum_upper_t = std::numeric_limits::signaling_NaN(); /// Coordinates of the packet, in element logical coordinates tnsr::I coordinates; @@ -70,6 +87,22 @@ struct Packet { void renormalize_momentum( const tnsr::II& inv_spatial_metric, const Scalar& lapse); + + void pup(PUP::er& p); + + /// Overloaded comparison operator. Useful for test; in an actual simulation + /// two distinct packets are not truly "identical" as they represent different + /// particles. + bool operator==(const Packet& rhs) const { + return (this->species == rhs.species) and + (this->number_of_neutrinos == rhs.number_of_neutrinos) and + (this->index_of_closest_grid_point == + rhs.index_of_closest_grid_point) and + (this->time == rhs.time) and + (this->momentum_upper_t == rhs.momentum_upper_t) and + (this->coordinates == rhs.coordinates) and + (this->momentum == rhs.momentum); + }; }; /*! diff --git a/src/Evolution/Particles/MonteCarlo/Tags.hpp b/src/Evolution/Particles/MonteCarlo/Tags.hpp new file mode 100644 index 000000000000..0b09a0a56f8a --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/Tags.hpp @@ -0,0 +1,31 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Evolution/Particles/MonteCarlo/Packet.hpp" + +/// Items related to the evolution of particles +/// Items related to Monte-Carlo radiation transport +/// Tags for MC +namespace Particles::MonteCarlo::Tags { + +/// Simple tag containing the vector of Monte-Carlo +/// packets belonging to an element. +struct PacketsOnElement : db::SimpleTag { + using type = std::vector; +}; + +/// Simple tag containing an approximation of the light +/// crossing time for each cell (the shortest time among +/// all coordinate axis directions). +template +struct CellLightCrossingTime : db::SimpleTag { + using type = Scalar; +}; + +} // namespace Particles::MonteCarlo::Tags diff --git a/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp b/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp index 41150439a133..00838bf3f48b 100644 --- a/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp +++ b/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp @@ -37,7 +37,7 @@ void combine_ghost_data( } } } - // Loop over each direction. We assume at most one neighber in each + // Loop over each direction. We assume at most one neighbor in each // direction. for (auto& [direction_id, ghost_data_dir] : ghost_data) { if (ghost_data_dir) { diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt b/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt index 9617f6673567..e4b1ed21fdc6 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt +++ b/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt @@ -6,6 +6,7 @@ set(LIBRARY "Test_MonteCarlo") set(LIBRARY_SOURCES Test_CellCrossingTime.cpp Test_CellVolume.cpp + Test_CommunicationTags.cpp Test_EmitPackets.cpp Test_EvolvePackets.cpp Test_ImplicitMonteCarloCorrections.cpp @@ -26,6 +27,7 @@ target_link_libraries( ${LIBRARY} PRIVATE DataStructures + DgSubcell GeneralRelativity GeneralRelativityHelpers H5 @@ -33,4 +35,5 @@ target_link_libraries( HydroHelpers Informer MonteCarlo + Time ) diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp new file mode 100644 index 000000000000..a91b6f4c8f5e --- /dev/null +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp @@ -0,0 +1,638 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Domain/Tags/NeighborMesh.hpp" +#include "Evolution/DgSubcell/ActiveGrid.hpp" +#include "Evolution/DgSubcell/Mesh.hpp" +#include "Evolution/DgSubcell/SliceData.hpp" +#include "Evolution/DgSubcell/SubcellOptions.hpp" +#include "Evolution/DgSubcell/Tags/ActiveGrid.hpp" +#include "Evolution/DgSubcell/Tags/Interpolators.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp" +#include "Evolution/Particles/MonteCarlo/MortarData.hpp" +#include "Evolution/Particles/MonteCarlo/Packet.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "Framework/ActionTesting.hpp" +#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Parallel/Phase.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Time/Slab.hpp" +#include "Time/Tags/TimeStepId.hpp" +#include "Time/Time.hpp" +#include "Time/TimeStepId.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/Numeric.hpp" +#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace { + +struct Var1 : db::SimpleTag { + using type = Scalar; +}; + +template +struct System { + static constexpr size_t volume_dim = Dim; + using variables_tag = ::Tags::Variables>; + using flux_variables = tmpl::list; +}; + +template +struct component { + using metavariables = Metavariables; + using chare_type = ActionTesting::MockArrayChare; + using array_index = ElementId; + + using initial_tags = tmpl::list< + ::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>, + domain::Tags::Mesh, evolution::dg::subcell::Tags::Mesh, + evolution::dg::subcell::Tags::ActiveGrid, domain::Tags::Element, + Particles::MonteCarlo::Tags::McGhostZoneDataTag, + Particles::MonteCarlo::Tags::PacketsOnElement, + ::Tags::Variables>, + Particles::MonteCarlo::Tags::MortarDataTag, + hydro::Tags::RestMassDensity, + hydro::Tags::ElectronFraction, + hydro::Tags::Temperature, + Particles::MonteCarlo::Tags::CellLightCrossingTime, + domain::Tags::NeighborMesh, + evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd>; + + using phase_dependent_action_list = tmpl::list, + Particles::MonteCarlo::Actions::SendDataForMcCommunication< + Dim, + // No local time stepping + false, Particles::MonteCarlo::CommunicationStep::PreStep>, + Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication< + Dim, Particles::MonteCarlo::CommunicationStep::PreStep>, + Particles::MonteCarlo::Actions::SendDataForMcCommunication< + Dim, + // No local time stepping + false, Particles::MonteCarlo::CommunicationStep::PostStep>, + Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication< + Dim, Particles::MonteCarlo::CommunicationStep::PostStep>>>>; +}; + +template +struct Metavariables { + static constexpr size_t volume_dim = Dim; + using component_list = tmpl::list>; + using system = System; + using const_global_cache_tags = tmpl::list<>; +}; + +template +void test_send_receive_actions() { + CAPTURE(Dim); + + using Interps = DirectionalIdMap>>; + using metavars = Metavariables; + using comp = component; + using MockRuntimeSystem = ActionTesting::MockRuntimeSystem; + MockRuntimeSystem runner{{}}; + + const TimeStepId time_step_id{true, 1, Time{Slab{1.0, 2.0}, {0, 10}}}; + const TimeStepId next_time_step_id{true, 1, Time{Slab{1.0, 2.0}, {1, 10}}}; + const Mesh dg_mesh{5, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh subcell_mesh = evolution::dg::subcell::fd::mesh(dg_mesh); + const evolution::dg::subcell::ActiveGrid active_grid = + evolution::dg::subcell::ActiveGrid::Subcell; + + // ^ eta + // +-+-+> xi + // |X| | + // +-+-+ + // | | | + // +-+-+ + // + // The "self_id" for the element that we are considering is marked by an X in + // the diagram. We consider a configuration with one neighbor in the +xi + // direction (east_id), and (in 2d and 3d) one in the -eta (south_id) + // direction. + // + // In 1d there aren't any projections to test, and in 3d we only have 1 + // element in the z-direction. + DirectionMap> neighbors{}; + ElementId self_id{}; + ElementId east_id{}; + // NOLINTNEXTLINE(misc-const-correctness) + ElementId south_id{}; // not used in 1d + // NOLINTNEXTLINE(misc-const-correctness) + OrientationMap orientation{}; + // Note: in 2d and 3d it is the neighbor in the lower eta direction that has a + // non-trivial orientation. + + if constexpr (Dim == 1) { + self_id = ElementId{0, {{{1, 0}}}}; + east_id = ElementId{0, {{{1, 1}}}}; + neighbors[Direction::upper_xi()] = Neighbors{{east_id}, + OrientationMap::create_aligned()}; + } else if constexpr (Dim == 2) { + orientation = OrientationMap{ + std::array{Direction::lower_xi(), Direction::lower_eta()}}; + self_id = ElementId{0, {{{1, 0}, {0, 0}}}}; + east_id = ElementId{0, {{{1, 1}, {0, 0}}}}; + south_id = ElementId{1, {{{0, 0}, {0, 0}}}}; + neighbors[Direction::upper_xi()] = Neighbors{{east_id}, + OrientationMap::create_aligned()}; + neighbors[Direction::lower_eta()] = + Neighbors{{south_id}, orientation}; + } else { + static_assert(Dim == 3, "Only implemented tests in 1, 2, and 3d"); + orientation = OrientationMap{std::array{Direction::lower_xi(), + Direction::lower_eta(), + Direction::upper_zeta()}}; + self_id = ElementId{0, {{{1, 0}, {0, 0}, {0, 0}}}}; + east_id = ElementId{0, {{{1, 1}, {0, 0}, {0, 0}}}}; + south_id = ElementId{1, {{{0, 0}, {0, 0}, {0, 0}}}}; + neighbors[Direction::upper_xi()] = Neighbors{{east_id}, + OrientationMap::create_aligned()}; + neighbors[Direction::lower_eta()] = + Neighbors{{south_id}, orientation}; + } + const Element element{self_id, neighbors}; + + using NeighborDataMap = + DirectionalIdMap>; + NeighborDataMap neighbor_data{}; + const DirectionalId east_neighbor_id{Direction::upper_xi(), + east_id}; + // insert data from one of the neighbors to make sure the send actions clears + // it. + neighbor_data[east_neighbor_id] = {}; + + using evolved_vars_tags = tmpl::list; + Variables evolved_vars{ + subcell_mesh.number_of_grid_points()}; + // Set Var1 to the logical coords, just need some data + get(get(evolved_vars)) = get<0>(logical_coordinates(subcell_mesh)); + + Scalar rest_mass_density( + get<0>(logical_coordinates(subcell_mesh))); + Scalar electron_fraction( + get<0>(logical_coordinates(subcell_mesh)) / 2.0); + Scalar temperature(get<0>(logical_coordinates(subcell_mesh)) * + 3.5); + Scalar cell_light_crossing_time( + get<0>(logical_coordinates(subcell_mesh)) * 2.0); + std::vector packets_on_element{}; + + for (const auto& [direction, neighbor_ids] : neighbors) { + (void)direction; + for (const auto& neighbor_id : neighbor_ids) { + // Initialize neighbors with garbage data. We won't ever run any actions + // on them, we just need to insert them to make sure things are sent to + // the right places. We'll check their inboxes directly. + ActionTesting::emplace_array_component_and_initialize( + &runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}, + neighbor_id, + {time_step_id, next_time_step_id, Mesh{}, Mesh{}, + active_grid, Element{}, NeighborDataMap{}, packets_on_element, + Variables{}, + typename Particles::MonteCarlo::Tags::MortarDataTag::type{}, + Scalar{}, Scalar{}, Scalar{}, + Scalar{}, + typename domain::Tags::NeighborMesh::type{}, Interps{}}); + } + } + + // Initialize MortarData as needed + Particles::MonteCarlo::MortarData mortar_data{}; + { + const size_t number_of_points_in_ghost_zone = + Dim > 1 ? + subcell_mesh.slice_away(0).number_of_grid_points() : + 1; + mortar_data.rest_mass_density[east_neighbor_id] = + DataVector(number_of_points_in_ghost_zone, 0.1); + mortar_data.electron_fraction[east_neighbor_id] = + DataVector(number_of_points_in_ghost_zone, 0.1); + mortar_data.temperature[east_neighbor_id] = + DataVector(number_of_points_in_ghost_zone, 0.1); + mortar_data.cell_light_crossing_time[east_neighbor_id] = + DataVector(number_of_points_in_ghost_zone, 0.1); + } + if constexpr (Dim > 1) { + const DirectionalId south_neighbor_id{Direction::lower_eta(), + south_id}; + const Mesh ghost_mesh = subcell_mesh.slice_away(1); + mortar_data.rest_mass_density[south_neighbor_id] = + DataVector(ghost_mesh.number_of_grid_points(), 0.1); + mortar_data.electron_fraction[south_neighbor_id] = + DataVector(ghost_mesh.number_of_grid_points(), 0.1); + mortar_data.temperature[south_neighbor_id] = + DataVector(ghost_mesh.number_of_grid_points(), 0.1); + mortar_data.cell_light_crossing_time[south_neighbor_id] = + DataVector(ghost_mesh.number_of_grid_points(), 0.1); + } + + const size_t species = 1; + const double number_of_neutrinos = 2.0; + const size_t index_of_closest_grid_point = 0; + const double t0 = 1.2; + const double x0 = 0.3; + const double y0 = 0.5; + const double z0 = -0.7; + const double p_upper_t0 = 1.1; + const double p_x0 = 0.9; + const double p_y0 = 0.7; + const double p_z0 = 0.1; + // Packet to be sent to east element + Particles::MonteCarlo::Packet packet_east( + species, number_of_neutrinos, index_of_closest_grid_point, t0, 1.3, y0, + z0, p_upper_t0, p_x0, p_y0, p_z0); + // Packet to be kept by current element + Particles::MonteCarlo::Packet packet_keep( + species, number_of_neutrinos, index_of_closest_grid_point, t0, x0, y0, z0, + p_upper_t0, p_x0, p_y0, p_z0); + // Packet to be deleted for being out of bounds + const Particles::MonteCarlo::Packet packet_delete( + species, number_of_neutrinos, index_of_closest_grid_point, t0, -1.3, y0, + z0, p_upper_t0, p_x0, p_y0, p_z0); + // Packet to be sent to south element + // Note that the packet extends outside of the current domain in both the y + // and z direction, but should be moved in the direction where it is the + // farthest from the boundary (y, here). + Particles::MonteCarlo::Packet packet_south( + species, number_of_neutrinos, index_of_closest_grid_point, t0, x0, -1.2, + -1.1, p_upper_t0, p_x0, p_y0, p_z0); + packets_on_element.push_back(packet_east); + packets_on_element.push_back(packet_keep); + packets_on_element.push_back(packet_delete); + if constexpr (Dim > 1) { + packets_on_element.push_back(packet_south); + } + + Interps fd_to_neighbor_fd_interpolants{}; + ActionTesting::emplace_array_component_and_initialize( + &runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}, self_id, + {time_step_id, next_time_step_id, dg_mesh, subcell_mesh, active_grid, + element, neighbor_data, packets_on_element, evolved_vars, mortar_data, + rest_mass_density, electron_fraction, temperature, + cell_light_crossing_time, + typename domain::Tags::NeighborMesh::type{}, + fd_to_neighbor_fd_interpolants}); + + using ghost_data_tag = Particles::MonteCarlo::Tags::McGhostZoneDataTag; + using ActionTesting::get_databox_tag; + CHECK(get_databox_tag(runner, self_id).size() == 1); + CHECK(get_databox_tag(runner, self_id) + .count(east_neighbor_id) == 1); + + // Run the SendDataForReconstruction action on self_id (PreStep) + ActionTesting::next_action(make_not_null(&runner), self_id); + CHECK(get_databox_tag(runner, self_id).empty()); + + // Check data sent to neighbors + const size_t ghost_zone_size = 1; + const auto& directions_to_slice = element.internal_boundaries(); + const DirectionMap all_sliced_data = + [&rest_mass_density, &electron_fraction, &temperature, + &cell_light_crossing_time, &subcell_mesh, ghost_zone_size, + &directions_to_slice, &fd_to_neighbor_fd_interpolants]() { + (void)ghost_zone_size; + const size_t dv_size = subcell_mesh.number_of_grid_points(); + const size_t number_of_vars = 4; + DataVector buffer{dv_size * number_of_vars}; + std::copy( + get(rest_mass_density).data(), + std::next(get(rest_mass_density).data(), static_cast(dv_size)), + buffer.data()); + std::copy( + get(electron_fraction).data(), + std::next(get(electron_fraction).data(), static_cast(dv_size)), + std::next(buffer.data(), static_cast(dv_size))); + std::copy(get(temperature).data(), + std::next(get(temperature).data(), static_cast(dv_size)), + std::next(buffer.data(), static_cast(dv_size * 2))); + std::copy(get(cell_light_crossing_time).data(), + std::next(get(cell_light_crossing_time).data(), + static_cast(dv_size)), + std::next(buffer.data(), static_cast(dv_size * 3))); + + return evolution::dg::subcell::slice_data( + buffer, subcell_mesh.extents(), ghost_zone_size, + directions_to_slice, 0, fd_to_neighbor_fd_interpolants); + }(); + + { + const auto& expected_east_data = + all_sliced_data.at(east_neighbor_id.direction()); + const auto& east_data = ActionTesting::get_inbox_tag< + comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PreStep>>( + runner, east_id); + CHECK(east_data.at(time_step_id) + .at(DirectionalId{Direction::lower_xi(), self_id}) + .ghost_zone_hydro_variables == expected_east_data); + CHECK(east_data.at(time_step_id) + .at(DirectionalId{Direction::lower_xi(), self_id}) + .packets_entering_this_element == std::nullopt); + } + + if constexpr (Dim > 1) { + const auto direction = Direction::lower_eta(); + const auto& expected_south_data = all_sliced_data.at(direction); + const auto& south_data = ActionTesting::get_inbox_tag< + comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PreStep>>( + runner, south_id); + CHECK( + south_data.at(time_step_id) + .at(DirectionalId{orientation(direction.opposite()), self_id}) + .ghost_zone_hydro_variables == expected_south_data); + CHECK(south_data.at(time_step_id) + .at(DirectionalId{direction, self_id}) + .packets_entering_this_element == std::nullopt); + } + + // Set the inbox data on self_id and then check that it gets processed + // correctly. + auto& self_inbox_pre = ActionTesting::get_inbox_tag< + comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PreStep>>( + make_not_null(&runner), self_id); + REQUIRE_FALSE(ActionTesting::next_action_if_ready( + make_not_null(&runner), self_id)); + + // Send data from east neighbor + DataVector east_ghost_cells{}; + { + const size_t number_of_vars = 4; + east_ghost_cells = + DataVector{subcell_mesh.slice_away(0).number_of_grid_points() * + ghost_zone_size * number_of_vars}; + alg::iota(east_ghost_cells, 2.0); + const size_t items_in_inbox = + Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PreStep>:: + insert_into_inbox( + make_not_null(&self_inbox_pre), time_step_id, + std::pair{ + DirectionalId{Direction::upper_xi(), east_id}, + Particles::MonteCarlo::McGhostZoneData{ + east_ghost_cells, std::nullopt}}); + CHECK(items_in_inbox == 1); + } + // NOLINTNEXTLINE(misc-const-correctness) + [[maybe_unused]] DataVector south_ghost_cells{}; + + if constexpr (Dim > 1) { + REQUIRE_FALSE(ActionTesting::next_action_if_ready( + make_not_null(&runner), self_id)); + + const size_t number_of_vars = 4; + south_ghost_cells = + DataVector{subcell_mesh.slice_away(1).number_of_grid_points() * + ghost_zone_size * number_of_vars}; + alg::iota(south_ghost_cells, 10000.0); + *std::prev(south_ghost_cells.end()) = -10.0; + const size_t items_in_inbox = + Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PreStep>:: + insert_into_inbox( + make_not_null(&self_inbox_pre), time_step_id, + std::pair{ + DirectionalId{Direction::lower_eta(), south_id}, + Particles::MonteCarlo::McGhostZoneData{ + south_ghost_cells, std::nullopt}}); + CHECK(items_in_inbox == 2); + } + + // Run the ReceiveDataForReconstruction action on self_id (PreStep) + ActionTesting::next_action(make_not_null(&runner), self_id); + + // Check the received data was stored correctly + using mortar_data_tag = Particles::MonteCarlo::Tags::MortarDataTag; + const auto& mortar_data_from_box = + get_databox_tag(runner, self_id); + + { + REQUIRE(mortar_data_from_box.rest_mass_density.find(east_neighbor_id) != + mortar_data_from_box.rest_mass_density.end()); + REQUIRE(mortar_data_from_box.electron_fraction.find(east_neighbor_id) != + mortar_data_from_box.electron_fraction.end()); + REQUIRE(mortar_data_from_box.temperature.find(east_neighbor_id) != + mortar_data_from_box.temperature.end()); + REQUIRE( + mortar_data_from_box.cell_light_crossing_time.find(east_neighbor_id) != + mortar_data_from_box.cell_light_crossing_time.end()); + + const size_t number_of_east_points = east_ghost_cells.size() / 4; + const DataVector rest_mass_density_view{east_ghost_cells.data(), + number_of_east_points}; + const DataVector electron_fraction_view{ + east_ghost_cells.data() + number_of_east_points, number_of_east_points}; + const DataVector temperature_view{ + east_ghost_cells.data() + number_of_east_points * 2, + number_of_east_points}; + const DataVector cell_light_crossing_time_view{ + east_ghost_cells.data() + number_of_east_points * 3, + east_ghost_cells.size() - number_of_east_points * 3}; + + CHECK( + mortar_data_from_box.rest_mass_density.find(east_neighbor_id)->second == + rest_mass_density_view); + CHECK( + mortar_data_from_box.electron_fraction.find(east_neighbor_id)->second == + electron_fraction_view); + CHECK(mortar_data_from_box.temperature.find(east_neighbor_id)->second == + temperature_view); + CHECK(mortar_data_from_box.cell_light_crossing_time.find(east_neighbor_id) + ->second == cell_light_crossing_time_view); + } + if constexpr (Dim > 1) { + const DirectionalId south_neighbor_id{Direction::lower_eta(), + south_id}; + REQUIRE(mortar_data_from_box.rest_mass_density.find(south_neighbor_id) != + mortar_data_from_box.rest_mass_density.end()); + REQUIRE(mortar_data_from_box.electron_fraction.find(south_neighbor_id) != + mortar_data_from_box.electron_fraction.end()); + REQUIRE(mortar_data_from_box.temperature.find(south_neighbor_id) != + mortar_data_from_box.temperature.end()); + REQUIRE( + mortar_data_from_box.cell_light_crossing_time.find(south_neighbor_id) != + mortar_data_from_box.cell_light_crossing_time.end()); + + const size_t number_of_south_points = south_ghost_cells.size() / 4; + const DataVector rest_mass_density_view{south_ghost_cells.data(), + number_of_south_points}; + const DataVector electron_fraction_view{ + south_ghost_cells.data() + number_of_south_points, + number_of_south_points}; + const DataVector temperature_view{ + south_ghost_cells.data() + number_of_south_points * 2, + number_of_south_points}; + const DataVector cell_light_crossing_time_view{ + south_ghost_cells.data() + number_of_south_points * 3, + south_ghost_cells.size() - number_of_south_points * 3}; + + CHECK(mortar_data_from_box.rest_mass_density.find(south_neighbor_id) + ->second == rest_mass_density_view); + CHECK(mortar_data_from_box.electron_fraction.find(south_neighbor_id) + ->second == electron_fraction_view); + CHECK(mortar_data_from_box.temperature.find(south_neighbor_id)->second == + temperature_view); + CHECK(mortar_data_from_box.cell_light_crossing_time.find(south_neighbor_id) + ->second == cell_light_crossing_time_view); + } + { + const auto& packets_from_box = + get_databox_tag( + runner, self_id); + CHECK(packets_from_box.size() == packets_on_element.size()); + } + + // Run the SendDataForReconstruction action on self_id (PostStep) + ActionTesting::next_action(make_not_null(&runner), self_id); + CHECK(get_databox_tag(runner, self_id).empty()); + + { + // Verify that the packets out of the element have been removed. + const auto& packets_from_box = + get_databox_tag( + runner, self_id); + CHECK(packets_from_box.size() == 1); + CHECK(packets_from_box[0] == packet_keep); + } + + // Correct coordinates of packets sent east/south to get in the + // topological coordinate of their new element. + packet_east.coordinates[0] -= 2.0; + packet_south.coordinates[1] += 2.0; + { + const auto& east_data = ActionTesting::get_inbox_tag< + comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PostStep>>( + runner, east_id); + CHECK(east_data.at(time_step_id) + .at(DirectionalId{Direction::lower_xi(), self_id}) + .packets_entering_this_element != std::nullopt); + CHECK(east_data.at(time_step_id) + .at(DirectionalId{Direction::lower_xi(), self_id}) + .packets_entering_this_element.value() + .size() == 1); + const Particles::MonteCarlo::Packet received_packet = + east_data.at(time_step_id) + .at(DirectionalId{Direction::lower_xi(), self_id}) + .packets_entering_this_element.value()[0]; + CHECK(packet_east == received_packet); + } + if constexpr (Dim > 1) { + const auto direction = Direction::lower_eta(); + const auto& south_data = ActionTesting::get_inbox_tag< + comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PostStep>>( + runner, south_id); + CHECK( + south_data.at(time_step_id) + .at(DirectionalId{orientation(direction.opposite()), self_id}) + .packets_entering_this_element != std::nullopt); + CHECK( + south_data.at(time_step_id) + .at(DirectionalId{orientation(direction.opposite()), self_id}) + .packets_entering_this_element.value() + .size() == 1); + const Particles::MonteCarlo::Packet received_packet = + south_data.at(time_step_id) + .at(DirectionalId{orientation(direction.opposite()), self_id}) + .packets_entering_this_element.value()[0]; + CHECK(packet_south == received_packet); + } + + // Set the inbox data on self_id and then check that it gets processed + // correctly. + auto& self_inbox_post = ActionTesting::get_inbox_tag< + comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PostStep>>( + make_not_null(&runner), self_id); + // Check that we are not ready to receive yet (Inboxes unfilled) + REQUIRE_FALSE(ActionTesting::next_action_if_ready( + make_not_null(&runner), self_id)); + + // Set up fake data coming from east neighbor + { + const std::optional> + packets_from_east = + std::vector{packet_east}; + const DataVector east_ghost_cells_post{}; + const size_t items_in_inbox = + Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PostStep>:: + insert_into_inbox( + make_not_null(&self_inbox_post), time_step_id, + std::pair{ + DirectionalId{Direction::upper_xi(), east_id}, + Particles::MonteCarlo::McGhostZoneData{ + east_ghost_cells_post, packets_from_east}}); + CHECK(items_in_inbox == 1); + } + // Set up fake data coming from south neighbor + if constexpr (Dim > 1) { + REQUIRE_FALSE(ActionTesting::next_action_if_ready( + make_not_null(&runner), self_id)); + + const std::optional> + packets_from_south = + std::vector{packet_south}; + const DataVector south_ghost_cells_post{}; + const size_t items_in_inbox = + Particles::MonteCarlo::McGhostZoneDataInboxTag< + Dim, Particles::MonteCarlo::CommunicationStep::PostStep>:: + insert_into_inbox( + make_not_null(&self_inbox_post), time_step_id, + std::pair{ + DirectionalId{Direction::lower_eta(), south_id}, + Particles::MonteCarlo::McGhostZoneData{ + south_ghost_cells_post, packets_from_south}}); + CHECK(items_in_inbox == 2); + } + // Run the ReceiveDataForReconstruction action on self_id (PostStep) + ActionTesting::next_action(make_not_null(&runner), self_id); + + // We should now have 3 packets (2 for Dim=1) + { + const auto& packets_from_box = + get_databox_tag( + runner, self_id); + CHECK(packets_from_box.size() == (Dim > 1 ? 3 : 2)); + CHECK(packets_from_box[0] == packet_keep); + CHECK(packets_from_box[1] == packet_east); + if constexpr (Dim > 1) { + CHECK(packets_from_box[2] == packet_south); + } + } +} + +} // namespace + +SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarlo.CommunicationTags", + "[Evolution][Unit]") { + test_send_receive_actions<1>(); + test_send_receive_actions<2>(); + test_send_receive_actions<3>(); +}