Merge pull request #5851 from ffoucart/DiffusionScatteringMC

Diffusion scattering mc

docs/References.bib

@@ -969,6 +969,21 @@ @article{Foucart2015vpa
   SLACcitation   = "%%CITATION = ARXIV:1502.04146;%%"
 }
 
+@article{Foucart:2017mbt,
+    author = "Foucart, Francois",
+    title = "{Monte Carlo closure for moment-based transport schemes in
+              general relativistic radiation hydrodynamic simulations}",
+    eprint = "1708.08452",
+    archivePrefix = "arXiv",
+    primaryClass = "astro-ph.HE",
+    doi = "10.1093/mnras/sty108",
+    journal = "Mon. Not. Roy. Astron. Soc.",
+    volume = "475",
+    number = "3",
+    pages = "4186--4207",
+    year = "2018"
+}
+
 @article{Foucart:2021mcb,
     author = "Foucart, Francois and Duez, Matthew D. and Hebert, Francois and
               Kidder, Lawrence E. and Kovarik, Phillip and Pfeiffer, Harald P.

src/Evolution/Particles/MonteCarlo/EvolvePackets.hpp

@@ -12,7 +12,7 @@
 
 namespace Frame {
 struct Fluid;
-} // namespace Frame
+}  // namespace Frame
 
 /// Items related to the evolution of particles
 /// Items related to Monte-Carlo radiation transport

src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp

@@ -91,6 +91,7 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
     const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
     const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
     const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
+    const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
     const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
     const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>& mesh_velocity,
     const InverseJacobian<DataVector, 3, Frame::ElementLogical,
@@ -103,6 +104,20 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
   // RNG from uniform distribution in [eps=1.e-100,1[
   // Eps used to avoid log(0).
   std::uniform_real_distribution<double> rng_uniform_eps_to_one(1.e-100, 1.0);
+
+  // Struct used for diffusion pre-computations
+  const DiffusionMonteCarloParameters diffusion_params;
+  auto prefactor_diffusion_time_step = make_with_value<DataVector>(lapse, 0.0);
+  auto prefactor_diffusion_four_velocity =
+      make_with_value<DataVector>(lapse, 0.0);
+  auto prefactor_diffusion_time_vector =
+      make_with_value<DataVector>(lapse, 0.0);
+
+  DiffusionPrecomputeForElement(
+      &prefactor_diffusion_time_vector, &prefactor_diffusion_four_velocity,
+      &prefactor_diffusion_time_step, lorentz_factor,
+      lower_spatial_four_velocity, lapse, shift, spatial_metric);
+
   // Mesh information. Currently assumes uniform grid without map
   const Index<3>& extents = mesh.extents();
   const std::array<double, 3> bottom_coord_mesh{mesh_coordinates.get(0)[0],
@@ -213,7 +228,15 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
         // The scatterig depth of 3.0 was found to be sufficient for diffusion
         // to be accurate (see Foucart 2018, 10.1093/mnras/sty108)
         if (scattering_optical_depth > 3.0) {
-          diffuse_packet(&packet, dt_min);
+          diffuse_packet(
+              &packet, random_number_generator, &fluid_frame_energy, dt_min,
+              diffusion_params, scattering_opacity, lorentz_factor,
+              lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift,
+              d_inv_spatial_metric, spatial_metric, inv_spatial_metric,
+              mesh_velocity, inverse_jacobian_logical_to_inertial,
+              inertial_to_fluid_jacobian, inertial_to_fluid_inverse_jacobian,
+              prefactor_diffusion_time_step, prefactor_diffusion_four_velocity,
+              prefactor_diffusion_time_vector);
         } else {
           // Low optical depth; perform scatterings one by one.
           do {
@@ -265,6 +288,11 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
                          gsl::at(dx_mesh, d) +
                      0.5);
     }
+    // In SpEC, we update a packet index after a full time step;
+    // only the opacities are updated mid-step. Decide whether to do
+    // the same here once we handle ghost zone (the main reason not to
+    // update is to limit the amount of ghost zone information we
+    // store.
     packet.index_of_closest_grid_point =
         closest_point_index_3d[0] +
         extents[0] * (closest_point_index_3d[1] +