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 "NumericalAlgorithms/Spectral/Mesh.hpp"
11 : #include "Utilities/Gsl.hpp"
12 :
13 : /// \cond
14 : namespace Frame {
15 : struct Fluid;
16 : } // namespace Frame
17 : /// \endcond
18 :
19 : /// Items related to the evolution of particles
20 : /// Items related to Monte-Carlo radiation transport
21 : namespace Particles::MonteCarlo {
22 :
23 : struct Packet;
24 :
25 : namespace detail {
26 :
27 : void draw_single_packet(
28 : gsl::not_null<double*> time,
29 : gsl::not_null<std::array<double, 3>*> coord,
30 : gsl::not_null<std::array<double, 3>*> momentum,
31 : gsl::not_null<std::mt19937*> random_number_generator);
32 : } // namespace detail
33 :
34 :
35 : /// Structure containing Monte-Carlo function templated on EnergyBins
36 : /// and/or NeutrinoSpecies
37 : template <size_t EnergyBins, size_t NeutrinoSpecies>
38 1 : struct TemplatedLocalFunctions {
39 : /*!
40 : * \brief Function emitting Monte Carlo packets
41 : *
42 : * We emit a total energy emission_in_cell for each grid
43 : * cell, neutrino species, and neutrino energy bin. We aim for packets of
44 : * energy single_packet_energy in the fluid frame. The number of packets
45 : * is, in each bin b and for each species s,
46 : * `N = emission_in_cell[s][b] / single_packet_energy[s]`
47 : * and randomly rounded up or down to get integer number of packets (i.e.
48 : * if we want N=2.6 packets, there is a 60 percent chance of creating a 3rd
49 : * packet). The position of the packets is drawn from a homogeneous
50 : * distribution in the logical frame, and the direction of propagation of the
51 : * neutrinos is drawn from an isotropic distribution in the fluid frame.
52 : * Specifically, the 4-momentum of a packet of energy nu in the fluid frame
53 : * is
54 : * \f{align}
55 : * p^t &= \nu \\
56 : * p^x &= \nu * sin(\theta) * cos(\phi) \\
57 : * p^y &= \nu * sin(\theta) * sin(\phi) \\
58 : * p^z &= \nu * cos(theta)
59 : * \f}
60 : * with \f$cos(\theta)\f$ drawn from a uniform distribution in [-1,1] and
61 : * \f$\phi\f$ from a uniform distribution in \f$[0,2*\pi]\f$. We
62 : * transform to the inertial frame \f$p_t\f$ and \f$p^x,p^y,p^z\f$ using the
63 : * jacobian/inverse jacobian passed as option. The number of neutrinos in
64 : * each packet is defined as
65 : * `n = single_packet_energy[s] / energy_at_bin_center[b]`
66 : * Note that the packet energy is in code units and energy of a bin in MeV.
67 : */
68 1 : void emit_packets(
69 : gsl::not_null<std::vector<Packet>*> packets,
70 : gsl::not_null<std::mt19937*> random_number_generator,
71 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
72 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> coupling_tilde_s,
73 : gsl::not_null<Scalar<DataVector>*> coupling_rho_ye,
74 : const double& time_start_step, const double& time_step,
75 : const Mesh<3>& mesh,
76 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
77 : emission_in_cell,
78 : const std::array<DataVector, NeutrinoSpecies>& single_packet_energy,
79 : const std::array<double, EnergyBins>& energy_at_bin_center,
80 : const Scalar<DataVector>& lorentz_factor,
81 : const tnsr::i<DataVector, 3, Frame::Inertial>&
82 : lower_spatial_four_velocity,
83 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
84 : inertial_to_fluid_jacobian,
85 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
86 : inertial_to_fluid_inverse_jacobian);
87 :
88 : /*!
89 : * \brief Evolve Monte-Carlo packets within an element.
90 : *
91 : * Evolve packets until (approximately) the provided final time
92 : * following the methods of \cite Foucart:2021mcb.
93 : * The vector of Packets should contain all MC packets that we wish
94 : * to advance in time. Note that this function only handles
95 : * propagation / absorption / scattering of packets, but not
96 : * emission. The final time of the packet may differ from the
97 : * desired final time by up to 5 percent, due to the fact that when
98 : * using the diffusion approximation (for large scattering
99 : * opacities) we take fixed time steps in the fluid frame,
100 : * leading to unpredictable time steps in the inertial frame.
101 : */
102 1 : void evolve_packets(
103 : gsl::not_null<std::vector<Packet>*> packets,
104 : gsl::not_null<std::mt19937*> random_number_generator,
105 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
106 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
107 : coupling_tilde_s,
108 : gsl::not_null<Scalar<DataVector>*> coupling_rho_ye,
109 : double final_time, const Mesh<3>& mesh,
110 : const tnsr::I<DataVector, 3, Frame::ElementLogical>& mesh_coordinates,
111 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
112 : absorption_opacity_table,
113 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
114 : scattering_opacity_table,
115 : const std::array<double, EnergyBins>& energy_at_bin_center,
116 : const Scalar<DataVector>& lorentz_factor,
117 : const tnsr::i<DataVector, 3, Frame::Inertial>&
118 : lower_spatial_four_velocity,
119 : const Scalar<DataVector>& lapse,
120 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
121 : const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
122 : const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
123 : const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
124 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
125 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
126 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
127 : mesh_velocity,
128 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
129 : Frame::Inertial>&
130 : inverse_jacobian_logical_to_inertial,
131 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
132 : inertial_to_fluid_jacobian,
133 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
134 : inertial_to_fluid_inverse_jacobian);
135 :
136 : /// Interpolate the opacities from values tabulated as a function of
137 : /// neutrino energy in the fluid frame, to the value at the given
138 : /// fluid frame energy.
139 : /// We interpolate the logarithm of the opacity, linearly in fluid
140 : /// frame energy.
141 1 : void interpolate_opacities_at_fluid_energy(
142 : gsl::not_null<double*> absorption_opacity_packet,
143 : gsl::not_null<double*> scattering_opacity_packet,
144 : double fluid_frame_energy, size_t species, size_t index,
145 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
146 : absorption_opacity_table,
147 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
148 : scattering_opacity_table,
149 : const std::array<double, EnergyBins>& energy_at_bin_center);
150 :
151 : // Floor for opacity values (absorption or scattering).
152 : // This is used in two places. First, when interpolating opacities
153 : // in energy space, we interpolate log(max(kappa,opacity_floor)),
154 : // with kappa the tabulated value. Second, if kappa = opacity_floor
155 : // for some interaction, we assume that the interaction never happens.
156 0 : const double opacity_floor = 1.e-100;
157 : };
158 :
159 : } // namespace Particles::MonteCarlo
|
