Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include "DataStructures/DataVector.hpp"
7 : #include "DataStructures/Tensor/Tensor.hpp"
8 : #include "Utilities/Gsl.hpp"
9 :
10 : namespace Frame {
11 : struct Fluid;
12 : } // namespace Frame
13 :
14 : /// Items related to the evolution of particles
15 : /// Items related to Monte-Carlo radiation transport
16 : namespace Particles::MonteCarlo {
17 :
18 : struct Packet;
19 :
20 : /// Precomputed quantities useful for the diffusion approximation
21 : /// in high-scattering opacity regions.
22 : /// We follow \cite Foucart2017mbt
23 : /// Note that r_d in that manuscript should be
24 : /// (distance diffused)/(time elapsed)
25 : /// i.e. the paper is missing a normalization
26 : /// of r_d by delta_t. The upper bound of the integral in Eq (30)
27 : /// should just be sqrt(3*tau/4)
28 : /// All quantities in this struct are fixed; we could also just
29 : /// hard-code them to save us the on-the-fly calculation, or
30 : /// have a single instance of this struct in the global cache.
31 1 : struct DiffusionMonteCarloParameters {
32 0 : DiffusionMonteCarloParameters();
33 :
34 0 : const double MinimumOpacityForDiffusion = 3.0;
35 0 : const double MaximumOpacityForCorrection = 11.0;
36 : /// Definition of the vector B_i, equation (35)
37 1 : const std::array<double, 21> BvsRhoForScattering{
38 : 1000., 18.74, 7.52, 4.75, 3.51, 2.78, 2.32, 2.00, 1.77, 1.60, 1.47,
39 : 1.36, 1.28, 1.21, 1.15, 1.10, 1.07, 1.04, 1.019, 1.0027, 1.0000001};
40 : /// Storage for the function r(P) implicitly defined by Eq (29)
41 : /// ScatteringRofP[i] = r(0.001*i)
42 : /// Calculation performed in constructor
43 1 : std::array<double, 1001> ScatteringRofP;
44 : /// Storage for the opacity dependent correction on the
45 : /// right-hand side of Eq (30)
46 : /// We use 101 points between the min and max opacities
47 : /// defined above.
48 : /// Calculation performed in constructor
49 1 : std::array<double, 101> OpacityDependentCorrectionToRofP;
50 : };
51 :
52 : /// Precompute quantities needed for evolving packets in the diffusion
53 : /// approximation, i.e. the transport velocity in logical coordinates
54 : /// and the coefficients defined by Eq (31-33)
55 1 : void DiffusionPrecomputeForElement(
56 : gsl::not_null<DataVector*> prefactor_diffusion_time_vector,
57 : gsl::not_null<DataVector*> prefactor_diffusion_four_velocity,
58 : gsl::not_null<DataVector*> prefactor_diffusion_time_step,
59 : const Scalar<DataVector>& lorentz_factor,
60 : const tnsr::i<DataVector, 3, Frame::Inertial>& lower_spatial_four_velocity,
61 : const Scalar<DataVector>& lapse,
62 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
63 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric);
64 :
65 : /// Elastic scattering, which is just redrawing the momentum from
66 : /// an isotropic distribution in the fluid frame, at constant
67 : /// fluid frame energy.
68 1 : std::array<double, 4> scatter_packet(
69 : gsl::not_null<Packet*> packet,
70 : gsl::not_null<std::mt19937*> random_number_generator,
71 : const double& fluid_frame_energy,
72 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
73 : inertial_to_fluid_jacobian,
74 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
75 : inertial_to_fluid_inverse_jacobian);
76 :
77 : /// Evolve a packet for dt = time_step, assuming that we can use
78 : /// the diffusion approximation.
79 : ///
80 : /// The coupling terms include ghost zone information, while other
81 : /// tensors only include live points.
82 : /// extended_idx is the index of the packet in DataVectors including
83 : /// ghost zones.
84 1 : void diffuse_packet(
85 : gsl::not_null<Packet*> packet,
86 : gsl::not_null<std::mt19937*> random_number_generator,
87 : gsl::not_null<double*> neutrino_energy,
88 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
89 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> coupling_tilde_s,
90 : gsl::not_null<Scalar<DataVector>*> coupling_rho_ye, size_t extended_idx,
91 : double time_step, const DiffusionMonteCarloParameters& diffusion_params,
92 : double absorption_opacity, double scattering_opacity,
93 : const Scalar<DataVector>& lorentz_factor,
94 : const tnsr::i<DataVector, 3, Frame::Inertial>& lower_spatial_four_velocity,
95 : const Scalar<DataVector>& lapse,
96 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
97 : const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
98 : const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
99 : const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
100 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
101 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
102 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>& mesh_velocity,
103 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
104 : Frame::Inertial>&
105 : inverse_jacobian_logical_to_inertial,
106 : const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
107 : inertial_to_fluid_jacobian,
108 : const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
109 : inertial_to_fluid_inverse_jacobian,
110 : const DataVector& prefactor_diffusion_time_step,
111 : const DataVector& prefactor_diffusion_four_velocity,
112 : const DataVector& prefactor_diffusion_time_vector);
113 :
114 : } // namespace Particles::MonteCarlo
|
