Infrastructure for MC communication

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
   )

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

src/Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp

@@ -0,0 +1,96 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <map>
+#include <optional>
+#include <utility>
+
+#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 <size_t Dim>
+struct McGhostZoneData {
+  void pup(PUP::er& p) {
+    p | ghost_zone_hydro_variables;
+    p | packets_entering_this_element;
+  }
+
+  DataVector ghost_zone_hydro_variables{};
+  std::optional<std::vector<Particles::MonteCarlo::Packet>>
+      packets_entering_this_element{};
+};
+
+/// Inbox tag to be used before MC step
+template <size_t Dim, CommunicationStep CommStep>
+struct McGhostZoneDataInboxTag {
+  using stored_type = McGhostZoneData<Dim>;
+
+  using temporal_id = TimeStepId;
+  using type = std::map<TimeStepId, DirectionalIdMap<Dim, stored_type>>;
+  using value_type = type;
+
+  CommunicationStep comm_step = CommStep;
+
+  template <typename ReceiveDataType>
+  static size_t insert_into_inbox(const gsl::not_null<type*> inbox,
+                                  const temporal_id& time_step_id,
+                                  ReceiveDataType&& data) {
+    auto& current_inbox = (*inbox)[time_step_id];
+    if (not current_inbox.insert(std::forward<ReceiveDataType>(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 <size_t Dim>
+struct McGhostZoneDataTag : db::SimpleTag {
+  using type = DirectionalIdMap<Dim, McGhostZoneData<Dim>>;
+};
+
+}  // namespace Tags
+
+} // namespace Particles::MonteCarlo

src/Evolution/Particles/MonteCarlo/MortarData.hpp

@@ -0,0 +1,47 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <optional>
+
+#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 <size_t Dim>
+struct MortarData {
+  DirectionalIdMap<Dim, std::optional<DataVector>> rest_mass_density{};
+  DirectionalIdMap<Dim, std::optional<DataVector>> electron_fraction{};
+  DirectionalIdMap<Dim, std::optional<DataVector>> temperature{};
+  DirectionalIdMap<Dim, std::optional<DataVector>> 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 <size_t Dim>
+struct MortarDataTag : db::SimpleTag {
+  using type = MortarData<Dim>;
+};
+
+}  // namespace Tags
+
+}  // namespace Particles::MonteCarlo

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<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
     const Scalar<DataVector>& lapse) {

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<size_t>::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<double>::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<size_t>::max();
 
   /// Current time
-  double time;
+  double time = std::numeric_limits<double>::signaling_NaN();
 
   /// Stores \f$p^t\f$
-  double momentum_upper_t;
+  double momentum_upper_t = std::numeric_limits<double>::signaling_NaN();
 
   /// Coordinates of the packet, in element logical coordinates
   tnsr::I<double, 3, Frame::ElementLogical> coordinates;
@@ -70,6 +87,22 @@ struct Packet {
   void renormalize_momentum(
       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
       const Scalar<DataVector>& 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);
+  };
 };
 
 /*!

src/Evolution/Particles/MonteCarlo/Tags.hpp

@@ -0,0 +1,31 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <vector>
+
+#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<Particles::MonteCarlo::Packet>;
+};
+
+/// Simple tag containing an approximation of the light
+/// crossing time for each cell (the shortest time among
+/// all coordinate axis directions).
+template <typename DataType>
+struct CellLightCrossingTime : db::SimpleTag {
+  using type = Scalar<DataType>;
+};
+
+}  // namespace Particles::MonteCarlo::Tags

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) {

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,11 +27,13 @@ target_link_libraries(
   ${LIBRARY}
   PRIVATE
   DataStructures
+  DgSubcell
   GeneralRelativity
   GeneralRelativityHelpers
   H5
   Hydro
   HydroHelpers
   Informer
   MonteCarlo
+  Time
   )