Mc optically thick sphere

src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt

@@ -33,6 +33,7 @@ target_link_libraries(
   LinearOperators
   MathFunctions
   MonteCarlo
+  MonteCarloSolutions
   Observer
   Options
   Parallel

src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp

@@ -20,6 +20,12 @@ void register_neutrino_tables() {
                  Particles::MonteCarlo::NeutrinoInteractionTable<16, 3>>{});
 }
 
+void register_mc_options() {
+  register_classes_with_charm(
+      tmpl::list<Particles::MonteCarlo::MonteCarloOptions<2>,
+                 Particles::MonteCarlo::MonteCarloOptions<3>>{});
+}
+
 extern "C" void CkRegisterMainModule() {
   Parallel::charmxx::register_main_module<EvolutionMetavars<4, 3>>();
   Parallel::charmxx::register_init_node_and_proc(
@@ -28,6 +34,6 @@ extern "C" void CkRegisterMainModule() {
        &domain::FunctionsOfTime::register_derived_with_charm,
        &EquationsOfState::register_derived_with_charm,
        &register_factory_classes_with_charm<EvolutionMetavars<4, 3>>,
-       &register_neutrino_tables},
+       &register_neutrino_tables, &register_mc_options},
       {});
 }

src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp

@@ -13,6 +13,8 @@
 #include "Evolution/ComputeTags.hpp"
 #include "Evolution/DgSubcell/Actions/Initialize.hpp"
 #include "Evolution/DgSubcell/BackgroundGrVars.hpp"
+#include "Evolution/DgSubcell/Tags/ObserverCoordinates.hpp"
+#include "Evolution/DgSubcell/Tags/ObserverMesh.hpp"
 #include "Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp"
 #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp"
 #include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp"
@@ -31,8 +33,11 @@
 #include "Evolution/Particles/MonteCarlo/Actions/Labels.hpp"
 #include "Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp"
 #include "Evolution/Particles/MonteCarlo/Actions/TriggerMonteCarloEvolution.hpp"
+#include "Evolution/Particles/MonteCarlo/CellCrossingTime.hpp"
 #include "Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp"
 #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
+#include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp"
+#include "Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp"
 #include "Evolution/Particles/MonteCarlo/System.hpp"
 #include "Evolution/Particles/MonteCarlo/Tags.hpp"
 #include "Evolution/Systems/GrMhd/ValenciaDivClean/AllSolutions.hpp"
@@ -59,6 +64,7 @@
 #include "ParallelAlgorithms/Actions/MutateApply.hpp"
 #include "ParallelAlgorithms/Actions/TerminatePhase.hpp"
 #include "ParallelAlgorithms/Events/Factory.hpp"
+#include "ParallelAlgorithms/Events/ObserveFields.hpp"
 #include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTrigger.hpp"
 #include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTriggers/Factory.hpp"
 #include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp"
@@ -70,7 +76,7 @@
 #include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
 #include "PointwiseFunctions/AnalyticData/Tags.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
-#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/ConstantM1.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Factory.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
 #include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp"
 #include "PointwiseFunctions/Hydro/Tags.hpp"
@@ -120,7 +126,7 @@ struct EvolutionMetavars {
   static constexpr bool use_dg_subcell = true;
 
   using initial_data_list =
-      grmhd::ValenciaDivClean::InitialData::initial_data_list;
+      RadiationTransport::MonteCarlo::Solutions::all_solutions;
   using equation_of_state_tag = hydro::Tags::GrmhdEquationOfState;
 
   struct SubcellOptions {
@@ -131,12 +137,20 @@ struct EvolutionMetavars {
   };
 
   using observe_fields =
-      tmpl::list<domain::Tags::Coordinates<volume_dim, Frame::Grid>,
-                 domain::Tags::Coordinates<volume_dim, Frame::Inertial>>;
+      tmpl::list<evolution::dg::subcell::Tags::ObserverCoordinatesCompute<
+                     volume_dim, Frame::ElementLogical>,
+                 evolution::dg::subcell::Tags::ObserverCoordinatesCompute<
+                     volume_dim, Frame::Grid>,
+                 evolution::dg::subcell::Tags::ObserverCoordinatesCompute<
+                     volume_dim, Frame::Inertial>,
+                 Particles::MonteCarlo::Tags::InertialFrameEnergyDensity>;
   using non_tensor_compute_tags =
-      tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>,
-                 ::Events::Tags::ObserverDetInvJacobianCompute<
-                     Frame::ElementLogical, Frame::Inertial>>;
+      tmpl::list<evolution::dg::subcell::Tags::ObserverMeshCompute<volume_dim>,
+                 evolution::dg::subcell::Tags::
+                  ObserverJacobianAndDetInvJacobianCompute<
+                      volume_dim, Frame::ElementLogical, Frame::Inertial>,
+                 evolution::dg::subcell::Tags::ObserverInverseJacobianCompute<
+                     volume_dim, Frame::ElementLogical, Frame::Inertial>>;
 
   using analytic_variables_tags = typename system::variables_tag::tags_list;
   using analytic_compute = evolution::Tags::AnalyticSolutionsCompute<
@@ -147,7 +161,13 @@ struct EvolutionMetavars {
     using factory_classes = tmpl::map<
         tmpl::pair<DenseTrigger, DenseTriggers::standard_dense_triggers>,
         tmpl::pair<DomainCreator<volume_dim>, domain_creators<volume_dim>>,
-        tmpl::pair<Event, tmpl::flatten<tmpl::list<Events::Completion>>>,
+        tmpl::pair<Event,
+                   tmpl::flatten<tmpl::list<
+                       Events::Completion,
+                       ::Events::ObserveNorms<observe_fields,
+                                              non_tensor_compute_tags>,
+                       dg::Events::ObserveFields<volume_dim, observe_fields,
+                                                 non_tensor_compute_tags>>>>,
         tmpl::pair<evolution::initial_data::InitialData, initial_data_list>,
         tmpl::pair<
             grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField,
@@ -188,6 +208,8 @@ struct EvolutionMetavars {
                                                 NeutrinoSpecies>,
       Initialization::Actions::AddComputeTags<tmpl::list<
           hydro::Tags::LowerSpatialFourVelocityCompute,
+          Particles::MonteCarlo::CellLightCrossingTimeCompute,
+          Particles::MonteCarlo::InertialFrameEnergyDensityCompute,
           Particles::MonteCarlo::InverseJacobianInertialToFluidCompute,
           domain::Tags::JacobianCompute<4, Frame::Inertial, Frame::Fluid>>>,
       evolution::Actions::InitializeRunEventsAndDenseTriggers,
@@ -198,6 +220,9 @@ struct EvolutionMetavars {
       tmpl::list<
           Parallel::PhaseActions<Parallel::Phase::Initialization,
                                  initialization_actions>,
+          Parallel::PhaseActions<Parallel::Phase::Register,
+                                 tmpl::list<dg_registration_list,
+                                            Parallel::Actions::TerminatePhase>>,
           Parallel::PhaseActions<
               Parallel::Phase::Evolve,
               tmpl::list<
@@ -241,6 +266,7 @@ struct EvolutionMetavars {
                  dg_element_array>;
 
   using const_global_cache_tags = tmpl::list<
+      Particles::MonteCarlo::Tags::MonteCarloOptions<NeutrinoSpecies>,
       equation_of_state_tag, evolution::initial_data::Tags::InitialData,
       Particles::MonteCarlo::Tags::InteractionRatesTable<EnergyBins,
                                                          NeutrinoSpecies>>;

src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

@@ -15,6 +15,7 @@
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
 #include "Evolution/Initialization/InitialData.hpp"
 #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
+#include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp"
 #include "Evolution/Particles/MonteCarlo/MortarData.hpp"
 #include "Evolution/Particles/MonteCarlo/Packet.hpp"
 #include "Evolution/Particles/MonteCarlo/Tags.hpp"
@@ -50,14 +51,16 @@ namespace Initialization::Actions {
 /// - evolution::dg::subcell::Tags::Mesh<dim>
 /// - evolution::dg::subcell::Tags::Coordinates<dim, Frame::Inertial>
 /// - evolution::initial_data::Tags::InitialData
+/// - Particles::MonteCarlo::Tags::MonteCarloOptions<EnergyBins,
+/// NeutrinoSpecies>
+/// - domain::Tags::Element<dim>
 ///
 /// DataBox changes:
 /// - Adds:
 ///   * Particles::MonteCarlo::Tags::PacketsOnElement
 ///   * Particles::MonteCarlo::Tags::RandomNumberGenerator
 ///   * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
 ///                                  NeutrinoSpecies>
-///   * Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>
 ///   * Background hydro variables
 ///   * Particles::MonteCarlo::Tags::MortarDataTag<dim>
 ///   * Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>
@@ -75,7 +78,6 @@ struct InitializeMCTags {
                  Particles::MonteCarlo::Tags::RandomNumberGenerator,
                  Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
                      NeutrinoSpecies>,
-                 Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
                  hydro_variables_tag,
                  Particles::MonteCarlo::Tags::MortarDataTag<dim>,
                  Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>,
@@ -96,9 +98,9 @@ struct InitializeMCTags {
         evolution::dg::subcell::ActiveGrid::Subcell) {
       ERROR("MC requires all elements to use Subcell");
     }
-    const size_t num_grid_points =
-        db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box)
-            .number_of_grid_points();
+    const Mesh<dim>& mesh =
+        db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box);
+    const size_t num_grid_points = mesh.number_of_grid_points();
     using derived_classes =
         tmpl::at<typename Metavariables::factory_creation::factory_classes,
                  evolution::initial_data::InitialData>;
@@ -120,6 +122,11 @@ struct InitializeMCTags {
               make_not_null(&box), std::move(hydro_variables));
         });
 
+    // Read global options for Monte-Carlo evolution
+    const auto mc_options = db::get<
+        Particles::MonteCarlo::Tags::MonteCarloOptions<NeutrinoSpecies>>(box);
+    const auto& initial_packet_energy = mc_options.get_initial_packet_energy();
+
     typename Particles::MonteCarlo::Tags::PacketsOnElement::type all_packets;
     Initialization::mutate_assign<
         tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>>(
@@ -132,26 +139,38 @@ struct InitializeMCTags {
         tmpl::list<Particles::MonteCarlo::Tags::RandomNumberGenerator>>(
         make_not_null(&box), std::move(rng));
 
-    // Currently hard-code energy at emission; should be set by option
+    // Initial energy of packets, read from MC options
     typename Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
         NeutrinoSpecies>::type packet_energy_at_emission =
         make_with_value<std::array<DataVector, NeutrinoSpecies>>(
-            DataVector{num_grid_points}, 1.e-12);
+            DataVector{num_grid_points}, 0.0);
+    for (size_t s = 0; s < NeutrinoSpecies; s++) {
+      packet_energy_at_emission[s] = initial_packet_energy[s];
+    }
     Initialization::mutate_assign<
         tmpl::list<Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
             NeutrinoSpecies>>>(make_not_null(&box),
                                std::move(packet_energy_at_emission));
 
-    typename Particles::MonteCarlo::Tags::CellLightCrossingTime<
-        DataVector>::type cell_light_crossing_time(num_grid_points, 1.0);
-    Initialization::mutate_assign<tmpl::list<
-        Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>>>(
-        make_not_null(&box), std::move(cell_light_crossing_time));
-
-    // Initialize empty mortar data; do we need more at initialization stage?
+    // Initialize mortar data. Currently assumes a single neighbor on each
+    // face (i.e. no h-refinement)
     using MortarData =
         typename Particles::MonteCarlo::Tags::MortarDataTag<dim>::type;
     MortarData mortar_data;
+    const Element<dim>& element = db::get<::domain::Tags::Element<dim>>(box);
+    for (const auto& [direction, neighbors] : element.neighbors()) {
+      const size_t sliced_mesh_size =
+          mesh.slice_away(direction.dimension()).number_of_grid_points();
+      DataVector zero_dv_ghost_zones(sliced_mesh_size, 0.0);
+      for (const auto& neighbor : neighbors) {
+        const DirectionalId<dim> mortar_id{direction, neighbor};
+        mortar_data.rest_mass_density.emplace(mortar_id, zero_dv_ghost_zones);
+        mortar_data.electron_fraction.emplace(mortar_id, zero_dv_ghost_zones);
+        mortar_data.temperature.emplace(mortar_id, zero_dv_ghost_zones);
+        mortar_data.cell_light_crossing_time.emplace(mortar_id,
+                                                     zero_dv_ghost_zones);
+      }
+    }
     Initialization::mutate_assign<
         tmpl::list<Particles::MonteCarlo::Tags::MortarDataTag<dim>>>(
         make_not_null(&box), std::move(mortar_data));

src/Evolution/Particles/MonteCarlo/CMakeLists.txt

@@ -13,7 +13,9 @@ spectre_target_sources(
   CouplingTermsForPropagation.cpp
   EvolvePackets.cpp
   InverseJacobianInertialToFluidCompute.cpp
+  MonteCarloOptions.cpp
   NeutrinoInteractionTable.cpp
+  NeutrinoMomentsFromMC.cpp
   Packet.cpp
   Scattering.cpp
   TemplatedLocalFunctions.cpp
@@ -34,8 +36,10 @@ spectre_target_headers(
   GhostZoneCommunicationTags.hpp
   ImplicitMonteCarloCorrections.tpp
   InverseJacobianInertialToFluidCompute.hpp
+  MonteCarloOptions.hpp
   MortarData.hpp
   NeutrinoInteractionTable.hpp
+  NeutrinoMomentsFromMC.hpp
   Packet.hpp
   Scattering.hpp
   System.hpp

src/Evolution/Particles/MonteCarlo/CellCrossingTime.cpp

@@ -25,6 +25,7 @@ void cell_light_crossing_time(
       inertial_coordinates.get(1)[step[1]] - inertial_coordinates.get(1)[0],
       inertial_coordinates.get(2)[step[2]] - inertial_coordinates.get(2)[0]};
 
+  get(*cell_light_crossing_time) = DataVector(n_pts,0.0);
   // Estimate light-crossing time in the cell.
   for (size_t i = 0; i < n_pts; i++) {
     double& min_crossing_time = get(*cell_light_crossing_time)[i];

src/Evolution/Particles/MonteCarlo/CellCrossingTime.hpp

@@ -4,6 +4,12 @@
 #pragma once
 
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/DgSubcell/Tags/Coordinates.hpp"
+#include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/Particles/MonteCarlo/Tags.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "Utilities/Gsl.hpp"
 
 /// \cond
 namespace gsl {
@@ -27,4 +33,28 @@ void cell_light_crossing_time(
     const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
     const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric);
 
+struct CellLightCrossingTimeCompute : Tags::CellLightCrossingTime<DataVector>,
+                                           db::ComputeTag {
+  using base = Tags::CellLightCrossingTime<DataVector>;
+  using return_type = typename base::type;
+  using argument_tags = tmpl::list<
+    evolution::dg::subcell::Tags::Mesh<3>,
+    evolution::dg::subcell::Tags::Coordinates<3, Frame::Inertial>,
+    gr::Tags::Lapse<DataVector>,
+    gr::Tags::Shift<DataVector, 3, Frame::Inertial>,
+    gr::Tags::InverseSpatialMetric<DataVector, 3, Frame::Inertial>>;
+
+  static void function(
+    gsl::not_null<return_type*> cell_light_crossing_time_,
+    const Mesh<3>& mesh,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coordinates,
+    const Scalar<DataVector>& lapse,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
+    const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric){
+    cell_light_crossing_time(cell_light_crossing_time_, mesh,
+        inertial_coordinates, lapse, shift, inv_spatial_metric);
+  }
+};
+
+
 }  // namespace Particles::MonteCarlo

src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp

@@ -23,6 +23,7 @@
 #include "Evolution/DgSubcell/GhostData.hpp"
 #include "Evolution/DgSubcell/SliceData.hpp"
 #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
+#include "Evolution/DgSubcell/Tags/Coordinates.hpp"
 #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
 #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp"
@@ -362,16 +363,31 @@ struct ReceiveDataForMcCommunication {
           },
           make_not_null(&box));
     } else {
+      const Mesh<Dim>& subcell_mesh =
+          db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box);
+      auto& mesh_coordinates =
+          db::get<evolution::dg::subcell::Tags::Coordinates<
+              Dim, Frame::ElementLogical>>(box);
       // TO DO: Deal with data coupling neutrinos back to fluid evolution
       db::mutate<Particles::MonteCarlo::Tags::PacketsOnElement>(
-          [&element, &received_data](
+          [&element, &received_data, &subcell_mesh, &mesh_coordinates](
               const gsl::not_null<std::vector<Particles::MonteCarlo::Packet>*>
                   packet_list) {
+            const Index<Dim>& extents = subcell_mesh.extents();
+            std::array<double, Dim> bottom_coord_mesh;
+            std::array<size_t, Dim> step;
+            std::array<double, Dim> dx_mesh;
+            for (size_t d = 0; d < Dim; d++) {
+              bottom_coord_mesh[d] = mesh_coordinates.get(d)[0];
+              step[d] = (d == 0) ? 1 : step[d - 1] * extents[d - 1];
+              dx_mesh[d] =
+                  mesh_coordinates.get(d)[step[d]] - bottom_coord_mesh[d];
+            }
             for (const auto& [direction, neighbors_in_direction] :
                  element.neighbors()) {
               for (const auto& neighbor : neighbors_in_direction) {
                 DirectionalId<Dim> directional_element_id{direction, neighbor};
-                const std::optional<std::vector<Particles::MonteCarlo::Packet>>&
+                std::optional<std::vector<Particles::MonteCarlo::Packet>>&
                     received_data_packets =
                         received_data[directional_element_id]
                             .packets_entering_this_element;
@@ -381,6 +397,26 @@ struct ReceiveDataForMcCommunication {
                 } else {
                   const size_t n_packets = received_data_packets.value().size();
                   for (size_t p = 0; p < n_packets; p++) {
+                    // Find closest grid point to received packet at current
+                    // time, using extents for live points only.
+                    Packet& packet = received_data_packets.value()[p];
+                    std::array<size_t, Dim> closest_point_index_3d;
+                    for (size_t d = 0; d < Dim; d++) {
+                      gsl::at(closest_point_index_3d, d) =
+                          std::floor((packet.coordinates[d] -
+                                      gsl::at(bottom_coord_mesh, d)) /
+                                         gsl::at(dx_mesh, d) +
+                                     0.5);
+                    }
+                    packet.index_of_closest_grid_point =
+                        closest_point_index_3d[0];
+                    size_t shift = extents[0];
+                    for (size_t d = 1; d < Dim; d++) {
+                      packet.index_of_closest_grid_point +=
+                          shift * closest_point_index_3d[d];
+                      shift *= extents[d];
+                    }
+                    // Now add packets to list in current element
                     packet_list->push_back(received_data_packets.value()[p]);
                   }
                 }
@@ -389,7 +425,6 @@ struct ReceiveDataForMcCommunication {
           },
           make_not_null(&box));
     }
-
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };