Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <optional>
       9             : #include <random>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/DataBox/DataBox.hpp"
      13             : #include "Domain/Tags.hpp"
      14             : #include "Domain/TagsTimeDependent.hpp"
      15             : #include "Evolution/DgSubcell/ActiveGrid.hpp"
      16             : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
      17             : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
      18             : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
      19             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      20             : #include "Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.hpp"
      21             : #include "Evolution/Particles/MonteCarlo/MortarData.hpp"
      22             : #include "Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp"
      23             : #include "Evolution/Particles/MonteCarlo/Tags.hpp"
      24             : #include "Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp"
      25             : #include "Parallel/AlgorithmExecution.hpp"
      26             : #include "Parallel/GlobalCache.hpp"
      27             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      28             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      29             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      30             : #include "Time/Tags/TimeStepId.hpp"
      31             : #include "Time/Time.hpp"
      32             : #include "Time/TimeStepId.hpp"
      33             : 
      34           0 : namespace Particles::MonteCarlo {
      35             : 
      36             : /// Mutator advancing neutrinos by a single step
      37             : template <size_t EnergyBins, size_t NeutrinoSpecies>
      38           1 : struct TimeStepMutator {
      39           0 :   static const size_t Dim = 3;
      40             : 
      41           0 :   using return_tags =
      42             :       tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement,
      43             :                  Particles::MonteCarlo::Tags::RandomNumberGenerator,
      44             :                  Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<3>>;
      45             :   // To do : check carefully DG vs Subcell quantities... everything should
      46             :   // be on the Subcell grid!
      47           0 :   using argument_tags = tmpl::list<
      48             :       ::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>,
      49             :       hydro::Tags::GrmhdEquationOfState,
      50             :       Particles::MonteCarlo::Tags::InteractionRatesTable<EnergyBins,
      51             :                                                          NeutrinoSpecies>,
      52             :       hydro::Tags::ElectronFraction<DataVector>,
      53             :       hydro::Tags::RestMassDensity<DataVector>,
      54             :       hydro::Tags::Temperature<DataVector>,
      55             :       hydro::Tags::LorentzFactor<DataVector>,
      56             :       hydro::Tags::LowerSpatialFourVelocity<DataVector, Dim, Frame::Inertial>,
      57             :       gr::Tags::Lapse<DataVector>,
      58             :       gr::Tags::Shift<DataVector, Dim, Frame::Inertial>,
      59             :       ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<Dim>,
      60             :                     Frame::Inertial>,
      61             :       ::Tags::deriv<gr::Tags::Shift<DataVector, Dim>, tmpl::size_t<Dim>,
      62             :                     Frame::Inertial>,
      63             :       ::Tags::deriv<gr::Tags::InverseSpatialMetric<DataVector, Dim>,
      64             :                     tmpl::size_t<Dim>, Frame::Inertial>,
      65             :       gr::Tags::SpatialMetric<DataVector, Dim, Frame::Inertial>,
      66             :       gr::Tags::InverseSpatialMetric<DataVector, Dim, Frame::Inertial>,
      67             :       gr::Tags::DetSpatialMetric<DataVector>,
      68             :       Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
      69             :       evolution::dg::subcell::Tags::Mesh<Dim>,
      70             :       evolution::dg::subcell::Tags::Coordinates<Dim, Frame::ElementLogical>,
      71             :       domain::Tags::MeshVelocity<Dim>,
      72             :       evolution::dg::subcell::fd::Tags::InverseJacobianLogicalToInertial<Dim>,
      73             :       evolution::dg::subcell::fd::Tags::DetInverseJacobianLogicalToInertial,
      74             :       domain::Tags::InverseJacobian<Dim + 1, Frame::Inertial, Frame::Fluid>,
      75             :       domain::Tags::Jacobian<Dim + 1, Frame::Inertial, Frame::Fluid>,
      76             :       Particles::MonteCarlo::Tags::MortarDataTag<Dim>>;
      77             : 
      78           0 :   static void apply(
      79             :       const gsl::not_null<std::vector<Packet>*> packets,
      80             :       const gsl::not_null<std::mt19937*> random_number_generator,
      81             :       const gsl::not_null<std::array<DataVector, NeutrinoSpecies>*>
      82             :           single_packet_energy,
      83             :       const TimeStepId& current_step_id, const TimeStepId& next_step_id,
      84             : 
      85             :       const EquationsOfState::EquationOfState<true, 3>& equation_of_state,
      86             :       const NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>&
      87             :           interaction_table,
      88             :       const Scalar<DataVector>& electron_fraction,
      89             :       const Scalar<DataVector>& rest_mass_density,
      90             :       const Scalar<DataVector>& temperature,
      91             :       const Scalar<DataVector>& lorentz_factor,
      92             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
      93             :           lower_spatial_four_velocity,
      94             :       const Scalar<DataVector>& lapse,
      95             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
      96             : 
      97             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& d_lapse,
      98             :       const tnsr::iJ<DataVector, Dim, Frame::Inertial>& d_shift,
      99             :       const tnsr::iJJ<DataVector, Dim, Frame::Inertial>& d_inv_spatial_metric,
     100             :       const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
     101             :       const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
     102             :       const Scalar<DataVector>& determinant_spatial_metric,
     103             :       const Scalar<DataVector>& cell_light_crossing_time, const Mesh<Dim>& mesh,
     104             :       const tnsr::I<DataVector, Dim, Frame::ElementLogical>& mesh_coordinates,
     105             :       const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
     106             :           mesh_velocity,
     107             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     108             :                             Frame::Inertial>&
     109             :           inverse_jacobian_logical_to_inertial,
     110             :       const Scalar<DataVector>& det_inverse_jacobian_logical_to_inertial,
     111             :       const InverseJacobian<DataVector, Dim + 1, Frame::Inertial, Frame::Fluid>&
     112             :           inertial_to_fluid_inverse_jacobian,
     113             :       const Jacobian<DataVector, Dim + 1, Frame::Inertial, Frame::Fluid>&
     114             :           inertial_to_fluid_jacobian,
     115             :       const MortarData<Dim>& mortar_data) {
     116             :     // Number of ghost zones for MC is assumed to be 1 for now.
     117             :     const size_t num_ghost_zones = 1;
     118             :     // Get information stored in various databox containers in
     119             :     // the format expected by take_time_step_on_element
     120             :     const double start_time = current_step_id.step_time().value();
     121             :     const double end_time = next_step_id.step_time().value();
     122             :     Scalar<DataVector> det_jacobian_logical_to_inertial(lapse);
     123             :     get(det_jacobian_logical_to_inertial) =
     124             :         1.0 / get(det_inverse_jacobian_logical_to_inertial);
     125             :     const DirectionalIdMap<Dim, std::optional<DataVector>>&
     126             :         electron_fraction_ghost = mortar_data.electron_fraction;
     127             :     const DirectionalIdMap<Dim, std::optional<DataVector>>&
     128             :         baryon_density_ghost = mortar_data.rest_mass_density;
     129             :     const DirectionalIdMap<Dim, std::optional<DataVector>>& temperature_ghost =
     130             :         mortar_data.temperature;
     131             :     const DirectionalIdMap<Dim, std::optional<DataVector>>&
     132             :         cell_light_crossing_time_ghost = mortar_data.cell_light_crossing_time;
     133             : 
     134             :     TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies> templated_functions;
     135             :     templated_functions.take_time_step_on_element(
     136             :         packets, random_number_generator, single_packet_energy, start_time,
     137             :         end_time, equation_of_state, interaction_table, electron_fraction,
     138             :         rest_mass_density, temperature, lorentz_factor,
     139             :         lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift,
     140             :         d_inv_spatial_metric, spatial_metric, inv_spatial_metric,
     141             :         determinant_spatial_metric, cell_light_crossing_time, mesh,
     142             :         mesh_coordinates, num_ghost_zones, mesh_velocity,
     143             :         inverse_jacobian_logical_to_inertial, det_jacobian_logical_to_inertial,
     144             :         inertial_to_fluid_jacobian, inertial_to_fluid_inverse_jacobian,
     145             :         electron_fraction_ghost, baryon_density_ghost, temperature_ghost,
     146             :         cell_light_crossing_time_ghost);
     147             :   }
     148             : };
     149             : 
     150           0 : namespace Actions {
     151             : 
     152             : /// Action taking a single time step of the Monte-Carlo evolution
     153             : /// algorithm, assuming that the fluid and metric data in the ghost
     154             : /// zones have been communicated and that packets are on the elements
     155             : /// that owns them.
     156             : template <size_t EnergyBins, size_t NeutrinoSpecies>
     157           1 : struct TakeTimeStep {
     158             :   template <typename DbTags, typename... InboxTags, typename ArrayIndex,
     159             :             typename ActionList, typename ParallelComponent,
     160             :             typename Metavariables>
     161           0 :   static Parallel::iterable_action_return_t apply(
     162             :       db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
     163             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
     164             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
     165             :       const ParallelComponent* const /*meta*/) {
     166             :     ASSERT(db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) ==
     167             :                evolution::dg::subcell::ActiveGrid::Subcell,
     168             :            "MC assumes that we are using the Subcell grid!");
     169             : 
     170             :     db::mutate_apply(TimeStepMutator<EnergyBins, NeutrinoSpecies>{},
     171             :                      make_not_null(&box));
     172             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     173             :   }
     174             : };
     175             : 
     176             : }  // namespace Actions
     177             : }  // namespace Particles::MonteCarlo