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 Foucart:2017mbt
      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           1 : void diffuse_packet(
      80             :     gsl::not_null<Packet*> packet,
      81             :     gsl::not_null<std::mt19937*> random_number_generator,
      82             :     gsl::not_null<double*> neutrino_energy,
      83             :     gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
      84             :     gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
      85             :       coupling_tilde_s,
      86             :     gsl::not_null<Scalar<DataVector>*> coupling_rho_ye,
      87             :     double time_step, const DiffusionMonteCarloParameters& diffusion_params,
      88             :     double absorption_opacity, double scattering_opacity,
      89             :     const Scalar<DataVector>& lorentz_factor,
      90             :     const tnsr::i<DataVector, 3, Frame::Inertial>& lower_spatial_four_velocity,
      91             :     const Scalar<DataVector>& lapse,
      92             :     const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
      93             :     const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
      94             :     const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
      95             :     const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
      96             :     const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
      97             :     const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
      98             :     const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>& mesh_velocity,
      99             :     const InverseJacobian<DataVector, 3, Frame::ElementLogical,
     100             :                           Frame::Inertial>&
     101             :         inverse_jacobian_logical_to_inertial,
     102             :     const Jacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
     103             :         inertial_to_fluid_jacobian,
     104             :     const InverseJacobian<DataVector, 4, Frame::Inertial, Frame::Fluid>&
     105             :         inertial_to_fluid_inverse_jacobian,
     106             :     const DataVector& prefactor_diffusion_time_step,
     107             :     const DataVector& prefactor_diffusion_four_velocity,
     108             :     const DataVector& prefactor_diffusion_time_vector);
     109             : 
     110             : }  // namespace Particles::MonteCarlo