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