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/DataBox/DataBox.hpp"
       7             : #include "DataStructures/DataVector.hpp"
       8             : #include "DataStructures/Index.hpp"
       9             : #include "DataStructures/Tensor/Tensor.hpp"
      10             : #include "Evolution/DgSubcell/ActiveGrid.hpp"
      11             : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
      12             : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
      13             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      14             : #include "Evolution/Particles/MonteCarlo/CellVolume.hpp"
      15             : #include "Evolution/Particles/MonteCarlo/Tags.hpp"
      16             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
      17             : #include "Parallel/AlgorithmExecution.hpp"
      18             : #include "Parallel/GlobalCache.hpp"
      19             : 
      20           0 : namespace Particles::MonteCarlo {
      21             : 
      22             : /// Mutator adding the Monte-Carlo contribution
      23             : /// to the evolution of the fluid.
      24           1 : struct FluidCouplingMutator {
      25           0 :   static const size_t Dim = 3;
      26             : 
      27             :   // We modify the fluid evolved variables, and reset the coupling
      28             :   // terms to zero.
      29           0 :   using return_tags =
      30             :       tmpl::list<grmhd::ValenciaDivClean::Tags::TildeTau,
      31             :                  grmhd::ValenciaDivClean::Tags::TildeYe,
      32             :                  grmhd::ValenciaDivClean::Tags::TildeS<Frame::Inertial>,
      33             :                  Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>,
      34             :                  Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>,
      35             :                  Particles::MonteCarlo::Tags::CouplingTildeS<DataVector, Dim>>;
      36           0 :   using argument_tags = tmpl::list<
      37             :       evolution::dg::subcell::Tags::Mesh<Dim>,
      38             :       evolution::dg::subcell::Tags::ActiveGrid,
      39             :       evolution::dg::subcell::fd::Tags::DetInverseJacobianLogicalToInertial>;
      40             : 
      41           0 :   static void apply(
      42             :       const gsl::not_null<Scalar<DataVector>*> tilde_tau,
      43             :       const gsl::not_null<Scalar<DataVector>*> tilde_ye,
      44             :       const gsl::not_null<tnsr::i<DataVector, Dim>*> tilde_s,
      45             :       const gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
      46             :       const gsl::not_null<Scalar<DataVector>*> coupling_tilde_rho_ye,
      47             :       const gsl::not_null<tnsr::i<DataVector, Dim>*> coupling_tilde_s,
      48             :       const Mesh<Dim>& mesh,
      49             :       const evolution::dg::subcell::ActiveGrid& active_grid,
      50             :       const Scalar<DataVector>& det_inverse_jacobian_logical_to_inertial) {
      51             :     // Currently, MC skips all non-communication actions when not using
      52             :     // Subcell.
      53             :     if (active_grid != evolution::dg::subcell::ActiveGrid::Subcell) {
      54             :       return;
      55             :     }
      56             : 
      57             :     // Coupling terms are on the mesh with ghost zones, while the evolved
      58             :     // variables are not. We also need to normalize the coupling terms with the
      59             :     // cell volume to get the evolution of energy/momentum density
      60             :     const Index<3>& extents = mesh.extents();
      61             :     const size_t num_ghost_zones = 1;
      62             :     const Index<3> extents_with_ghost{extents[0] + 2 * num_ghost_zones,
      63             :                                       extents[1] + 2 * num_ghost_zones,
      64             :                                       extents[2] + 2 * num_ghost_zones};
      65             : 
      66             :     Scalar<DataVector> det_jacobian_logical_to_inertial(*tilde_tau);
      67             :     get(det_jacobian_logical_to_inertial) =
      68             :         1.0 / get(det_inverse_jacobian_logical_to_inertial);
      69             :     Scalar<DataVector> cell_inertial_three_volume =
      70             :         make_with_value<Scalar<DataVector>>(*tilde_tau, 0.0);
      71             :     cell_inertial_coordinate_three_volume_finite_difference(
      72             :         &cell_inertial_three_volume, mesh, det_jacobian_logical_to_inertial);
      73             : 
      74             :     for (size_t i = 0; i < extents[0]; i++) {
      75             :       for (size_t j = 0; j < extents[1]; j++) {
      76             :         for (size_t k = 0; k < extents[2]; k++) {
      77             :           // The coupling terms are computed on a grid with ghost zone points,
      78             :           // the fluid variables are on the grid without GZ points
      79             :           // Index without GZ
      80             :           const size_t local_idx = collapsed_index(Index<3>{i, j, k}, extents);
      81             :           // Index with GZ
      82             :           const size_t extended_idx =
      83             :               collapsed_index(Index<3>{i + num_ghost_zones, j + num_ghost_zones,
      84             :                                        k + num_ghost_zones},
      85             :                               extents_with_ghost);
      86             :           // The MC coupling calculates the change in energy, momentum, and
      87             :           // lepton number. We need to divide by the cell 3-volume to get
      88             :           // the change in the evolved variables.
      89             :           const double& volume = get(cell_inertial_three_volume)[local_idx];
      90             :           get(*tilde_tau)[local_idx] +=
      91             :               get(*coupling_tilde_tau)[extended_idx] / volume;
      92             :           get(*tilde_ye)[local_idx] +=
      93             :               get(*coupling_tilde_rho_ye)[extended_idx] / volume;
      94             :           for (size_t d = 0; d < 3; d++) {
      95             :             tilde_s->get(d)[local_idx] +=
      96             :                 coupling_tilde_s->get(d)[extended_idx] / volume;
      97             :           }
      98             :         }
      99             :       }
     100             :     }
     101             :     // Reset coupling terms to 0 after use
     102             :     get(*coupling_tilde_tau) = 0.0;
     103             :     get(*coupling_tilde_rho_ye) = 0.0;
     104             :     for (size_t d = 0; d < 3; d++) {
     105             :       coupling_tilde_s->get(d) = 0.0;
     106             :     }
     107             :   }
     108             : };
     109             : 
     110             : }  // namespace Particles::MonteCarlo