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