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