MonteCarlo emission and TemplatedFunction structure

src/Evolution/Particles/MonteCarlo/CMakeLists.txt

@@ -8,18 +8,21 @@ add_spectre_library(${LIBRARY})
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
-  EvolveMonteCarloPackets.cpp
+  EvolvePackets.cpp
   InverseJacobianInertialToFluidCompute.cpp
-  MonteCarloPacket.cpp
+  TemplatedLocalFunctions.cpp
+  Packet.cpp
   )
 
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
-  EvolveMonteCarloPackets.hpp
+  EmitPackets.tpp
+  EvolvePackets.hpp
   InverseJacobianInertialToFluidCompute.hpp
-  MonteCarloPacket.hpp
+  TemplatedLocalFunctions.hpp
+  Packet.hpp
   )
 
 target_link_libraries(

src/Evolution/Particles/MonteCarlo/EmitPackets.tpp

@@ -0,0 +1,143 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include "Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp"
+
+#include <cmath>
+
+#include "Evolution/Particles/MonteCarlo/Packet.hpp"
+
+namespace Particles::MonteCarlo {
+
+template <size_t EnergyBins, size_t NeutrinoSpecies>
+void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::emit_packets(
+    const gsl::not_null<std::vector<Packet>*> packets,
+    const gsl::not_null<std::mt19937*> random_number_generator,
+    const double& time_start_step, const double& time_step,
+    const Mesh<3>& mesh,
+    const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
+        emission_in_cell,
+    const std::array<DataVector, NeutrinoSpecies>& single_packet_energy,
+    const std::array<double, EnergyBins>& energy_at_bin_center,
+    const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
+        inertial_to_fluid_jacobian,
+    const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
+        inertial_to_fluid_inverse_jacobian) {
+  std::uniform_real_distribution<double> rng_uniform_zero_to_one(0.0, 1.0);
+  // We begin by determining how many packets to create for each cell, neutrino
+  // species, and energy group.
+  const size_t grid_size = mesh.number_of_grid_points();
+  std::array<std::array<std::vector<size_t>, EnergyBins>, NeutrinoSpecies>
+      number_of_packets_to_create_per_cell;
+  size_t number_of_packets_to_create_total = 0;
+  for (size_t s = 0; s < NeutrinoSpecies; s++) {
+    for (size_t g = 0; g < EnergyBins; g++) {
+      gsl::at(gsl::at(number_of_packets_to_create_per_cell, s), g)
+          .resize(grid_size);
+      for (size_t i = 0; i < grid_size; i++) {
+        const double packets_to_create_double =
+            gsl::at(gsl::at(emission_in_cell, s), g)[i] /
+            gsl::at(single_packet_energy, s)[i];
+        size_t packets_to_create_int = floor(packets_to_create_double);
+        if (rng_uniform_zero_to_one(*random_number_generator) <
+            packets_to_create_double -
+                static_cast<double>(packets_to_create_int)) {
+          packets_to_create_int++;
+        }
+        gsl::at(gsl::at(number_of_packets_to_create_per_cell, s), g)[i] =
+            packets_to_create_int;
+        number_of_packets_to_create_total += packets_to_create_int;
+      }
+    }
+  }
+
+  // With the number of packets known, resize the vector holding the packets
+  Packet default_packet(0, 0.0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+  const size_t initial_size = packets->size();
+  packets->resize(initial_size + number_of_packets_to_create_total,
+                  default_packet);
+
+  // Allocate memory for some temporary data.
+  double time_normalized = 0.0;
+  std::array<double, 3> coord_normalized{{0.0, 0.0, 0.0}};
+  std::array<double, 3> three_momentum_normalized{{0.0, 0.0, 0.0}};
+  std::array<double, 4> four_momentum_fluid_frame{{0.0, 0.0, 0.0, 0.0}};
+  const Index<3> extents = mesh.extents();
+  const std::array<double, 3> logical_dx{2.0 / static_cast<double>(extents[0]),
+                                         2.0 / static_cast<double>(extents[1]),
+                                         2.0 / static_cast<double>(extents[2])};
+  std::array<double, 3> coord_bottom_cell{{0.0, 0.0, 0.0}};
+
+  // Now create the packets
+  size_t next_packet_index = initial_size;
+  for (size_t x_idx = 0; x_idx < extents[0]; x_idx++) {
+    gsl::at(coord_bottom_cell, 0) =
+        -1.0 + 2.0 * static_cast<double>(x_idx) * gsl::at(logical_dx, 0);
+    for (size_t y_idx = 0; y_idx < extents[1]; y_idx++) {
+      gsl::at(coord_bottom_cell, 1) =
+          -1.0 + 2.0 * static_cast<double>(y_idx) * gsl::at(logical_dx, 1);
+      for (size_t z_idx = 0; z_idx < extents[2]; z_idx++) {
+        const size_t idx = mesh.storage_index(Index<3>{{x_idx, y_idx, z_idx}});
+        gsl::at(coord_bottom_cell, 2) =
+            -1.0 + 2.0 * static_cast<double>(z_idx) * gsl::at(logical_dx, 2);
+        for (size_t s = 0; s < NeutrinoSpecies; s++) {
+          for (size_t g = 0; g < EnergyBins; g++) {
+            for (size_t p = 0;
+                 p < gsl::at(gsl::at(number_of_packets_to_create_per_cell, s),
+                             g)[idx];
+                 p++) {
+              detail::draw_single_packet(&time_normalized, &coord_normalized,
+                                         &three_momentum_normalized,
+                                         random_number_generator);
+              Packet& current_packet = (*packets)[next_packet_index];
+              current_packet.species = s;
+              current_packet.time =
+                  time_start_step + time_normalized * time_step;
+              current_packet.number_of_neutrinos =
+                  gsl::at(gsl::at(single_packet_energy, s), idx) /
+                  gsl::at(energy_at_bin_center, g);
+              current_packet.index_of_closest_grid_point = idx;
+              // 4_momentum_fluid_frame contains p^mu in an orthornormal
+              // coordinate system moving with the fluid.
+              gsl::at(four_momentum_fluid_frame, 0) =
+                  gsl::at(energy_at_bin_center, g);
+              for (size_t d = 0; d < 3; d++) {
+                current_packet.coordinates.get(d) =
+                    gsl::at(coord_bottom_cell, d) +
+                    gsl::at(coord_normalized, d) * gsl::at(logical_dx, d);
+                gsl::at(four_momentum_fluid_frame, d + 1) =
+                    gsl::at(three_momentum_normalized, d) *
+                    gsl::at(four_momentum_fluid_frame, 0);
+              }
+
+              // The packet stores p^t and p_i in the inertial frame. We
+              // transform the momentum accordingly.
+              current_packet.momentum_upper_t = 0.0;
+              for (size_t d = 0; d < 4; d++) {
+                current_packet.momentum_upper_t +=
+                    inertial_to_fluid_inverse_jacobian.get(0, d)[idx] *
+                    gsl::at(four_momentum_fluid_frame, d);
+              }
+              for (size_t d = 0; d < 3; d++) {
+                // Multiply by -1.0 because p_t = -p^t in an orthonormal frame
+                current_packet.momentum.get(d) =
+                    (-1.0) * gsl::at(four_momentum_fluid_frame, 0) *
+                    inertial_to_fluid_jacobian.get(d + 1, 0)[idx];
+                for (size_t dd = 0; dd < 3; dd++) {
+                  current_packet.momentum.get(d) +=
+                      gsl::at(four_momentum_fluid_frame, dd + 1) *
+                      inertial_to_fluid_jacobian.get(d + 1, dd + 1)[idx];
+                }
+              }
+              next_packet_index++;
+            }
+          }
+        }
+      }
+    }
+  }
+}
+
+}  // namespace Particles::MonteCarlo

src/Evolution/Particles/MonteCarlo/EvolveMonteCarloPackets.cpp → src/Evolution/Particles/MonteCarlo/EvolvePackets.cpp

@@ -1,9 +1,9 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
-#include "EvolveMonteCarloPackets.hpp"
+#include "EvolvePackets.hpp"
 
-#include "Evolution/Particles/MonteCarlo/MonteCarloPacket.hpp"
+#include "Evolution/Particles/MonteCarlo/Packet.hpp"
 
 namespace Particles::MonteCarlo {
 

src/Evolution/Particles/MonteCarlo/EvolveMonteCarloPackets.hpp → src/Evolution/Particles/MonteCarlo/EvolvePackets.hpp

@@ -42,6 +42,8 @@ void evolve_single_packet_on_geodesic(
 
 }  // namespace detail
 
+/// Evolve all packets in the provided std::vector (approximately) to the
+/// end of the current time step.
 void evolve_packets(
     gsl::not_null<std::vector<Packet>*> packets, const double& time_step,
     const Mesh<3>& mesh,

src/Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.hpp

@@ -17,9 +17,10 @@ struct Inertial;
 
 namespace Particles::MonteCarlo {
 
-// Inverse Jacobian of the map from inertial coordinate to an orthonormal frame
-// comoving with the fluid. That frame uses the 4-velocity as its time axis, and
-// constructs the other members of the tetrads using Gram-Schmidt's algorithm.
+/// Inverse Jacobian of the map from inertial coordinate to an orthonormal frame
+/// comoving with the fluid. That frame uses the 4-velocity as its time axis,
+/// and constructs the other members of the tetrads using Gram-Schmidt's
+/// algorithm.
 struct InverseJacobianInertialToFluidCompute
     : domain::Tags::InverseJacobian<4, Frame::Inertial, Frame::Fluid>,
       db::ComputeTag {

src/Evolution/Particles/MonteCarlo/MonteCarloPacket.cpp → src/Evolution/Particles/MonteCarlo/Packet.cpp

@@ -1,7 +1,7 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
-#include "MonteCarloPacket.hpp"
+#include "Evolution/Particles/MonteCarlo/Packet.hpp"
 
 namespace Particles::MonteCarlo {
 

src/Evolution/Particles/MonteCarlo/MonteCarloPacket.hpp → src/Evolution/Particles/MonteCarlo/Packet.hpp

@@ -13,14 +13,17 @@
 /// Items related to Monte-Carlo radiation transport
 namespace Particles::MonteCarlo {
 
+/// Struct representing a single Monte Carlo packet of neutrinos
 struct Packet {
-  // Constructor
-  Packet(const double& number_of_neutrinos_,
+  /// Constructor
+  Packet(const size_t& species_,
+         const double& number_of_neutrinos_,
          const size_t& index_of_closest_grid_point_, const double& time_,
          const double& coord_x_, const double& coord_y_, const double& coord_z_,
          const double& p_upper_t_, const double& p_x_, const double& p_y_,
          const double& p_z_)
-      : number_of_neutrinos(number_of_neutrinos_),
+      : species(species_),
+        number_of_neutrinos(number_of_neutrinos_),
         index_of_closest_grid_point(index_of_closest_grid_point_),
         time(time_),
         momentum_upper_t(p_upper_t_) {
@@ -32,29 +35,33 @@ struct Packet {
     momentum[2] = p_z_;
   }
 
-  // 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!
+  /// 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;
+
+  /// 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;
 
-  // Index of the closest point on the FD grid.
+  /// Index of the closest point on the FD grid.
   size_t index_of_closest_grid_point;
 
-  // Current time
+  /// Current time
   double time;
 
-  // p^t
+  /// Stores \f$p^t\f$
   double momentum_upper_t;
 
-  // Coordinates of the packet, currently in Inertial coordinates
+  /// Coordinates of the packet, in element logical coordinates
   tnsr::I<double, 3, Frame::ElementLogical> coordinates;
 
-  // Spatial components of the 4-momentum, also in Inertial coordinates
+  /// Spatial components of the 4-momentum \f$p_i\f$, in Inertial coordinates
   tnsr::i<double, 3, Frame::Inertial> momentum;
 
-  // Recalculte p^t using the fact that the 4-momentum is a null vector
+  /// Recalculte \f$p^t\f$ using the fact that the 4-momentum is a null vector
   void renormalize_momentum(
       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
       const Scalar<DataVector>& lapse);

src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.cpp

@@ -0,0 +1,37 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp"
+
+#include <cmath>
+
+#include "Evolution/Particles/MonteCarlo/EmitPackets.tpp"
+#include "Evolution/Particles/MonteCarlo/Packet.hpp"
+
+namespace Particles::MonteCarlo::detail {
+// Draw a single packet with homogeneouse spatial distribution
+// in [0,1]x[0,1]x[0,1], energy of 1, isotropic momentum
+// distribution, and time in [0,1].
+void draw_single_packet(
+    const gsl::not_null<double*> time,
+    const gsl::not_null<std::array<double, 3>*> coord,
+    const gsl::not_null<std::array<double, 3>*> momentum,
+    const gsl::not_null<std::mt19937*> random_number_generator) {
+  std::uniform_real_distribution<double> rng_uniform_zero_to_one(0.0, 1.0);
+
+  *time = rng_uniform_zero_to_one(*random_number_generator);
+  for (size_t i = 0; i < 3; i++) {
+    gsl::at(*coord, i) = rng_uniform_zero_to_one(*random_number_generator);
+  }
+  const double cos_theta =
+      -1.0 + 2.0 * rng_uniform_zero_to_one(*random_number_generator);
+  const double phi =
+      2.0 * M_PI * rng_uniform_zero_to_one(*random_number_generator);
+  const double sin_theta = sqrt(1.0 - cos_theta * cos_theta);
+  gsl::at(*momentum, 0) = sin_theta * cos(phi);
+  gsl::at(*momentum, 1) = sin_theta * sin(phi);
+  gsl::at(*momentum, 2) = cos_theta;
+}
+}  // namespace Particles::MonteCarlo::detail
+
+template struct Particles::MonteCarlo::TemplatedLocalFunctions<2, 2>;