Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <vector>
       7             : 
       8             : #include "DataStructures/DataVector.hpp"
       9             : #include "DataStructures/Tensor/Tensor.hpp"
      10             : #include "Utilities/Gsl.hpp"
      11             : 
      12             : /// Items related to the evolution of particles
      13             : /// Items related to Monte-Carlo radiation transport
      14             : namespace Particles::MonteCarlo {
      15             : 
      16             : /// Struct representing a single Monte Carlo packet of neutrinos
      17           1 : struct Packet {
      18             :   /// Constructor
      19           1 :   Packet(const size_t& species_,
      20             :          const double& number_of_neutrinos_,
      21             :          const size_t& index_of_closest_grid_point_, const double& time_,
      22             :          const double& coord_x_, const double& coord_y_, const double& coord_z_,
      23             :          const double& p_upper_t_, const double& p_x_, const double& p_y_,
      24             :          const double& p_z_)
      25             :       : species(species_),
      26             :         number_of_neutrinos(number_of_neutrinos_),
      27             :         index_of_closest_grid_point(index_of_closest_grid_point_),
      28             :         time(time_),
      29             :         momentum_upper_t(p_upper_t_) {
      30             :     coordinates[0] = coord_x_;
      31             :     coordinates[1] = coord_y_;
      32             :     coordinates[2] = coord_z_;
      33             :     momentum[0] = p_x_;
      34             :     momentum[1] = p_y_;
      35             :     momentum[2] = p_z_;
      36             :   }
      37             : 
      38             :   /// Default constructor needed to make Packet puppable
      39             :   /// For the same reason, we have default initialization
      40             :   /// for a number of member variables to unphysical values.
      41           1 :   Packet()
      42             :       : species{0},
      43             :         number_of_neutrinos{0},
      44             :         index_of_closest_grid_point{0},
      45             :         time{-1.0},
      46             :         momentum_upper_t{0.0} {
      47             :     coordinates[0] = 0.0;
      48             :     coordinates[1] = 0.0;
      49             :     coordinates[2] = 0.0;
      50             :     momentum[0] = 0.0;
      51             :     momentum[1] = 0.0;
      52             :     momentum[2] = 0.0;
      53             :   }
      54             : 
      55             :   /// Species of neutrinos (in the code, just an index used to access the
      56             :   /// right interaction rates; typically \f$0=\nu_e, 1=\nu_a, 2=\nu_x\f$)
      57           1 :   size_t species = std::numeric_limits<size_t>::max();
      58             : 
      59             :   /// Number of neutrinos represented by current packet
      60             :   /// Note that this number is rescaled so that
      61             :   /// `Energy_of_packet = N * Energy_of_neutrinos`
      62             :   /// with the packet energy in G=Msun=c=1 units but
      63             :   /// the neutrino energy in MeV!
      64           1 :   double number_of_neutrinos = std::numeric_limits<double>::signaling_NaN();
      65             : 
      66             :   /// Index of the closest point on the FD grid.
      67           1 :   size_t index_of_closest_grid_point = std::numeric_limits<size_t>::max();
      68             : 
      69             :   /// Current time
      70           1 :   double time = std::numeric_limits<double>::signaling_NaN();
      71             : 
      72             :   /// Stores \f$p^t\f$
      73           1 :   double momentum_upper_t = std::numeric_limits<double>::signaling_NaN();
      74             : 
      75             :   /// Coordinates of the packet, in element logical coordinates
      76           1 :   tnsr::I<double, 3, Frame::ElementLogical> coordinates;
      77             : 
      78             :   /// Spatial components of the 4-momentum \f$p_i\f$, in Inertial coordinates
      79           1 :   tnsr::i<double, 3, Frame::Inertial> momentum;
      80             : 
      81             :   /*!
      82             :    * Recalculte \f$p^t\f$ using the fact that the 4-momentum is a null vector
      83             :    * \f{align}{
      84             :    * p^t = \sqrt{\gamma^{ij} p_i p_j}/\alpha
      85             :    * \f}
      86             :    */
      87           1 :   void renormalize_momentum(
      88             :       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
      89             :       const Scalar<DataVector>& lapse);
      90             : 
      91           0 :   void pup(PUP::er& p);
      92             : 
      93             :   /// Overloaded comparison operator. Useful for test; in an actual simulation
      94             :   /// two distinct packets are not truly "identical" as they represent different
      95             :   /// particles.
      96           1 :   bool operator==(const Packet& rhs) const {
      97             :     return (this->species == rhs.species) and
      98             :            (this->number_of_neutrinos == rhs.number_of_neutrinos) and
      99             :            (this->index_of_closest_grid_point ==
     100             :             rhs.index_of_closest_grid_point) and
     101             :            (this->time == rhs.time) and
     102             :            (this->momentum_upper_t == rhs.momentum_upper_t) and
     103             :            (this->coordinates == rhs.coordinates) and
     104             :            (this->momentum == rhs.momentum);
     105             :   };
     106             : };
     107             : 
     108             : /*!
     109             :  * Calculate energy of neutrinos in a frame comoving with the fluid
     110             :  *
     111             :  * \f{align}{
     112             :  * E = W \alpha p^t - \gamma^{ij} u_i p_j
     113             :  * \f}
     114             :  */
     115           1 : double compute_fluid_frame_energy(
     116             :     const Packet& packet, const Scalar<DataVector>& lorentz_factor,
     117             :     const tnsr::i<DataVector, 3, Frame::Inertial>& lower_spatial_four_velocity,
     118             :     const Scalar<DataVector>& lapse,
     119             :     const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric);
     120             : 
     121             : }  // namespace Particles::MonteCarlo