Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <optional>
8 : #include <random>
9 : #include <vector>
10 :
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 :
13 : /// \cond
14 : class DataVector;
15 :
16 : template <size_t Dim, typename T>
17 : class DirectionalIdMap;
18 :
19 : namespace Frame {
20 : struct Fluid;
21 : } // namespace Frame
22 :
23 : namespace gsl {
24 : template <typename T>
25 : class not_null;
26 : } // namespace gsl
27 :
28 : template<size_t Dim>
29 : class Mesh;
30 :
31 : namespace Particles::MonteCarlo {
32 : template<size_t EnergyBins,size_t NeutrinoSpecies>
33 : class NeutrinoInteractionTable;
34 :
35 : struct Packet;
36 : } // namespace Particles::MonteCarlo
37 :
38 : namespace EquationsOfState {
39 : template <bool IsRelativistic, size_t ThermodynamicDim>
40 : class EquationOfState;
41 : } // namespace EquationsOfState
42 : /// \endcond
43 :
44 : /// Items related to the evolution of particles
45 : /// Items related to Monte-Carlo radiation transport
46 : namespace Particles::MonteCarlo {
47 :
48 : namespace detail {
49 :
50 : void draw_single_packet(gsl::not_null<double*> time,
51 : gsl::not_null<std::array<double, 3>*> coord,
52 : gsl::not_null<std::array<double, 3>*> momentum,
53 : gsl::not_null<std::mt19937*> random_number_generator);
54 : } // namespace detail
55 :
56 :
57 : /// Structure containing Monte-Carlo function templated on EnergyBins
58 : /// and/or NeutrinoSpecies
59 : template <size_t EnergyBins, size_t NeutrinoSpecies>
60 1 : struct TemplatedLocalFunctions {
61 : /*!
62 : * \brief Function to take a single Monte Carlo time step on a
63 : * finite difference element.
64 : */
65 1 : void take_time_step_on_element(
66 : gsl::not_null<std::vector<Packet>*> packets,
67 : gsl::not_null<std::mt19937*> random_number_generator,
68 : gsl::not_null<std::array<DataVector, NeutrinoSpecies>*>
69 : single_packet_energy,
70 :
71 : double start_time, double target_end_time,
72 : const EquationsOfState::EquationOfState<true, 3>& equation_of_state,
73 : const NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>&
74 : interaction_table,
75 : const Scalar<DataVector>& electron_fraction,
76 : const Scalar<DataVector>& rest_mass_density,
77 : const Scalar<DataVector>& temperature,
78 : const Scalar<DataVector>& lorentz_factor,
79 : const tnsr::i<DataVector, 3, Frame::Inertial>&
80 : lower_spatial_four_velocity,
81 : const Scalar<DataVector>& lapse,
82 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
83 : const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
84 : const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
85 : const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
86 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
87 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
88 : const Scalar<DataVector>& determinant_spatial_metric,
89 : const Scalar<DataVector>& cell_light_crossing_time, const Mesh<3>& mesh,
90 : const tnsr::I<DataVector, 3, Frame::ElementLogical>& mesh_coordinates,
91 : size_t num_ghost_zones,
92 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
93 : mesh_velocity,
94 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
95 : Frame::Inertial>&
96 : inverse_jacobian_logical_to_inertial,
97 : const Scalar<DataVector>& det_jacobian_logical_to_inertial,
98 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
99 : inertial_to_fluid_jacobian,
100 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
101 : inertial_to_fluid_inverse_jacobian,
102 : const DirectionalIdMap<3, std::optional<DataVector>>&
103 : electron_fraction_ghost,
104 : const DirectionalIdMap<3, std::optional<DataVector>>&
105 : baryon_density_ghost,
106 : const DirectionalIdMap<3, std::optional<DataVector>>& temperature_ghost,
107 : const DirectionalIdMap<3, std::optional<DataVector>>&
108 : cell_light_crossing_time_ghost);
109 :
110 : /*!
111 : * \brief Function emitting Monte Carlo packets
112 : *
113 : * We emit a total energy emissivity_in_cell * cell_proper_four_volume
114 : * for each grid
115 : * cell, neutrino species, and neutrino energy bin. We aim for packets of
116 : * energy single_packet_energy in the fluid frame. The number of packets
117 : * is, in each bin b and for each species s,
118 : * `N = emission_in_cell[s][b] / single_packet_energy[s]`
119 : * and randomly rounded up or down to get integer number of packets (i.e.
120 : * if we want N=2.6 packets, there is a 60 percent chance of creating a 3rd
121 : * packet). The position of the packets is drawn from a homogeneous
122 : * distribution in the logical frame, and the direction of propagation of the
123 : * neutrinos is drawn from an isotropic distribution in the fluid frame.
124 : * Specifically, the 4-momentum of a packet of energy nu in the fluid frame
125 : * is
126 : * \f{align}
127 : * p^t &= \nu \\
128 : * p^x &= \nu * sin(\theta) * cos(\phi) \\
129 : * p^y &= \nu * sin(\theta) * sin(\phi) \\
130 : * p^z &= \nu * cos(theta)
131 : * \f}
132 : * with \f$cos(\theta)\f$ drawn from a uniform distribution in [-1,1] and
133 : * \f$\phi\f$ from a uniform distribution in \f$[0,2*\pi]\f$. We
134 : * transform to the inertial frame \f$p_t\f$ and \f$p^x,p^y,p^z\f$ using the
135 : * jacobian/inverse jacobian passed as option. The number of neutrinos in
136 : * each packet is defined as
137 : * `n = single_packet_energy[s] / energy_at_bin_center[b]`
138 : * Note that the packet energy is in code units and energy of a bin in MeV.
139 : *
140 : * All tensors are assumed to correspond to live points only, except for
141 : * the coupling terms and emissivity, which include ghost zones.
142 : */
143 1 : void emit_packets(
144 : gsl::not_null<std::vector<Packet>*> packets,
145 : gsl::not_null<std::mt19937*> random_number_generator,
146 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
147 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> coupling_tilde_s,
148 : gsl::not_null<Scalar<DataVector>*> coupling_rho_ye,
149 : const double& time_start_step, const double& time_step,
150 : const Mesh<3>& mesh, size_t num_ghost_zones,
151 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
152 : emissivity_in_cell,
153 : const std::array<DataVector, NeutrinoSpecies>& single_packet_energy,
154 : const std::array<double, EnergyBins>& energy_at_bin_center,
155 : const Scalar<DataVector>& lorentz_factor,
156 : const tnsr::i<DataVector, 3, Frame::Inertial>&
157 : lower_spatial_four_velocity,
158 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
159 : inertial_to_fluid_jacobian,
160 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
161 : inertial_to_fluid_inverse_jacobian,
162 : const Scalar<DataVector>& cell_proper_four_volume);
163 :
164 : /*!
165 : * \brief Evolve Monte-Carlo packets within an element.
166 : *
167 : * Evolve packets until (approximately) the provided final time
168 : * following the methods of \cite Foucart2021mcb.
169 : * The vector of Packets should contain all MC packets that we wish
170 : * to advance in time. Note that this function only handles
171 : * propagation / absorption / scattering of packets, but not
172 : * emission. The final time of the packet may differ from the
173 : * desired final time by up to 5 percent, due to the fact that when
174 : * using the diffusion approximation (for large scattering
175 : * opacities) we take fixed time steps in the fluid frame,
176 : * leading to unpredictable time steps in the inertial frame.
177 : *
178 : * The absorption and scattering opacity tables include ghost
179 : * zones, and so do the coupling terms. Other variables are
180 : * only using live points.
181 : */
182 1 : void evolve_packets(
183 : gsl::not_null<std::vector<Packet>*> packets,
184 : gsl::not_null<std::mt19937*> random_number_generator,
185 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
186 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> coupling_tilde_s,
187 : gsl::not_null<Scalar<DataVector>*> coupling_rho_ye, double final_time,
188 : const Mesh<3>& mesh,
189 : const tnsr::I<DataVector, 3, Frame::ElementLogical>& mesh_coordinates,
190 : size_t num_ghost_zones,
191 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
192 : absorption_opacity_table,
193 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
194 : scattering_opacity_table,
195 : const std::array<double, EnergyBins>& energy_at_bin_center,
196 : const Scalar<DataVector>& lorentz_factor,
197 : const tnsr::i<DataVector, 3, Frame::Inertial>&
198 : lower_spatial_four_velocity,
199 : const Scalar<DataVector>& lapse,
200 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
201 : const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
202 : const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
203 : const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
204 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
205 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
206 : const Scalar<DataVector>& cell_light_crossing_time,
207 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
208 : mesh_velocity,
209 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
210 : Frame::Inertial>&
211 : inverse_jacobian_logical_to_inertial,
212 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
213 : inertial_to_fluid_jacobian,
214 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
215 : inertial_to_fluid_inverse_jacobian);
216 :
217 : /// Interpolate the opacities from values tabulated as a function of
218 : /// neutrino energy in the fluid frame, to the value at the given
219 : /// fluid frame energy.
220 : /// We interpolate the logarithm of the opacity, linearly in fluid
221 : /// frame energy.
222 1 : void interpolate_opacities_at_fluid_energy(
223 : gsl::not_null<double*> absorption_opacity_packet,
224 : gsl::not_null<double*> scattering_opacity_packet,
225 : double fluid_frame_energy, size_t species, size_t index,
226 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
227 : absorption_opacity_table,
228 : const std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>&
229 : scattering_opacity_table,
230 : const std::array<double, EnergyBins>& energy_at_bin_center);
231 :
232 : /// Function responsible to correct emissivity, absorption, and scattering
233 : /// rates in regions where there is a stiff coupling between the neutrinos
234 : /// and the fluid. This is done by transferring a fraction
235 : /// fraction_ka_to_ks of the absorption opacity to scattering opacity,
236 : /// while multiplying the emissivity by ( 1 - fraction_ka_to_ks ) to
237 : /// keep the equilibrium energy density constant.
238 1 : void implicit_monte_carlo_interaction_rates(
239 : gsl::not_null<
240 : std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
241 : emissivity_in_cell,
242 : gsl::not_null<
243 : std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
244 : absorption_opacity,
245 : gsl::not_null<
246 : std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
247 : scattering_opacity,
248 : gsl::not_null<
249 : std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
250 : fraction_ka_to_ks,
251 : const Scalar<DataVector>& cell_light_crossing_time,
252 : const Scalar<DataVector>& electron_fraction,
253 : const Scalar<DataVector>& rest_mass_density,
254 : const Scalar<DataVector>& temperature, double minimum_temperature,
255 : const NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>&
256 : interaction_table,
257 : const EquationsOfState::EquationOfState<true, 3>& equation_of_state);
258 :
259 : // Floor for opacity values (absorption or scattering).
260 : // This is used in two places. First, when interpolating opacities
261 : // in energy space, we interpolate log(max(kappa,opacity_floor)),
262 : // with kappa the tabulated value. Second, if kappa = opacity_floor
263 : // for some interaction, we assume that the interaction never happens.
264 0 : const double opacity_floor = 1.e-100;
265 : };
266 :
267 : } // namespace Particles::MonteCarlo
|