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