MC executable running a few steps (communication errors remain)

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

@@ -128,11 +128,9 @@ struct EvolutionMetavars {
     static constexpr bool subcell_enabled_at_external_boundary = true;
   };
 
-  using observe_fields = tmpl::list<
-      // tmpl::push_back<
-      // tmpl::append<typename system::mc_variables_tag::tags_list>,
-      domain::Tags::Coordinates<volume_dim, Frame::Grid>,
-      domain::Tags::Coordinates<volume_dim, Frame::Inertial>>;
+  using observe_fields =
+      tmpl::list<domain::Tags::Coordinates<volume_dim, Frame::Grid>,
+                 domain::Tags::Coordinates<volume_dim, Frame::Inertial>>;
   using non_tensor_compute_tags =
       tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>,
                  ::Events::Tags::ObserverDetInvJacobianCompute<
@@ -147,63 +145,35 @@ 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,
-                dg::Events::field_observations<
-                    volume_dim, observe_fields,
-                    non_tensor_compute_tags>  //,
-                                              // Events::time_events<system>
-                >>>,
+        tmpl::pair<Event, tmpl::flatten<tmpl::list<Events::Completion>>>,
         tmpl::pair<evolution::initial_data::InitialData, initial_data_list>,
         tmpl::pair<
             grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField,
             grmhd::AnalyticData::InitialMagneticFields::
                 initial_magnetic_fields>,
-        tmpl::pair<PhaseChange, PhaseControl::factory_creatable_classes>,
+        tmpl::pair<PhaseChange,
+                   tmpl::list<PhaseControl::CheckpointAndExitAfterWallclock>>,
         tmpl::pair<StepChooser<StepChooserUse::Slab>,
                    StepChoosers::standard_slab_choosers<system, false, false>>,
         tmpl::pair<TimeSequence<double>,
                    TimeSequences::all_time_sequences<double>>,
         tmpl::pair<TimeSequence<std::uint64_t>,
                    TimeSequences::all_time_sequences<std::uint64_t>>,
         tmpl::pair<TimeStepper, TimeSteppers::time_steppers>,
-        tmpl::pair<Trigger, tmpl::append<Triggers::logical_triggers,
-                                         Triggers::time_triggers>>>;
+        tmpl::pair<Trigger, Triggers::time_triggers>>;
   };
 
   using observed_reduction_data_tags =
       observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
           tmpl::at<typename factory_creation::factory_classes, Event>>>>;
 
-  using step_actions = tmpl::flatten<tmpl::list<
-      // Actions::RecordTimeStepperData<system>,
-      // Particles::MonteCarlo::Actions::SendDataForMcCommunication<
-      //   volume_dim,
-      //  No local time stepping
-      //    false, Particles::MonteCarlo::CommunicationStep::PreStep>,
-      // Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication<
-      //    volume_dim, Particles::MonteCarlo::CommunicationStep::PreStep>,
-      Particles::MonteCarlo::Actions::TakeTimeStep<EnergyBins, NeutrinoSpecies>,
-      // Particles::MonteCarlo::Actions::SendDataForMcCommunication<
-      //     volume_dim,
-      //  No local time stepping
-      //    false, Particles::MonteCarlo::CommunicationStep::PostStep>,
-      // Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication<
-      //    volume_dim, Particles::MonteCarlo::CommunicationStep::PostStep>,
-      // Actions::RecordTimeStepperData<system>,
-      // evolution::Actions::RunEventsAndDenseTriggers<tmpl::list<>>,
-      Actions::CleanHistory<system, false>>>;
-
   using dg_registration_list =
       tmpl::list<observers::Actions::RegisterEventsWithObservers>;
 
   using initialization_actions = tmpl::list<
       Initialization::Actions::InitializeItems<
           Initialization::TimeStepping<EvolutionMetavars, TimeStepperBase>,
           evolution::dg::Initialization::Domain<volume_dim>  //,
-          // Initialization::TimeStepperHistory<EvolutionMetavars>
           >,
       Initialization::Actions::AddSimpleTags<
           evolution::dg::BackgroundGrVars<system, EvolutionMetavars, true>>,
@@ -219,27 +189,23 @@ struct EvolutionMetavars {
           hydro::Tags::LowerSpatialFourVelocityCompute,
           Particles::MonteCarlo::InverseJacobianInertialToFluidCompute,
           domain::Tags::JacobianCompute<4, Frame::Inertial, Frame::Fluid>>>,
-      // Initialization::TimeStepperHistory<EvolutionMetavars>
+      evolution::Actions::InitializeRunEventsAndDenseTriggers,
       Parallel::Actions::TerminatePhase>;
 
   using dg_element_array = DgElementArray<
       EvolutionMetavars,
       tmpl::list<
           Parallel::PhaseActions<Parallel::Phase::Initialization,
                                  initialization_actions>,
-          // Parallel::PhaseActions<
-          //     Parallel::Phase::InitializeTimeStepperHistory,
-          //     SelfStart::self_start_procedure<step_actions, system>>,
-
-          //  Parallel::PhaseActions<Parallel::Phase::Register,
-          //                         tmpl::list<dg_registration_list,
-          //                 Parallel::Actions::TerminatePhase>>,
-
           Parallel::PhaseActions<
               Parallel::Phase::Evolve,
               tmpl::list<
-                  Actions::AdvanceTime, Actions::AdvanceTime,
                   Actions::AdvanceTime,
+                  evolution::Actions::RunEventsAndTriggers,
+                  Actions::AdvanceTime,
+                  evolution::Actions::RunEventsAndTriggers,
+                  Actions::AdvanceTime,
+                  evolution::Actions::RunEventsAndTriggers,
                   Particles::MonteCarlo::Actions::SendDataForMcCommunication<
                       volume_dim,
                       // No local time stepping
@@ -257,12 +223,7 @@ struct EvolutionMetavars {
                   Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication<
                       volume_dim,
                       Particles::MonteCarlo::CommunicationStep::PostStep>,
-                  Parallel::Actions::TerminatePhase>>
-          //      tmpl::list<evolution::Actions::RunEventsAndTriggers,
-          //                 Actions::ChangeSlabSize, //step_actions,
-          //                 Actions::AdvanceTime,
-          //                 PhaseControl::Actions::ExecutePhaseChange>>
-          >>;
+                  PhaseControl::Actions::ExecutePhaseChange>>>>;
 
   struct registration
       : tt::ConformsTo<Parallel::protocols::RegistrationMetavariables> {

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

@@ -3,9 +3,11 @@
 
 #pragma once
 
+#include <random>
 #include <tuple>
 
 #include "DataStructures/DataBox/DataBox.hpp"
+#include "Domain/Structure/Element.hpp"
 #include "Domain/Tags.hpp"
 #include "Evolution/Initialization/InitialData.hpp"
 #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
@@ -115,8 +117,10 @@ struct InitializeMCTags {
         tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>>(
         make_not_null(&box), std::move(all_packets));
 
-    // Currently seeds with 0 for testing...
-    typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(0);
+    unsigned long seed =
+        std::random_device{}();  // static_cast<unsigned long>(time(NULL));
+    typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(seed);
+
     Initialization::mutate_assign<
         tmpl::list<Particles::MonteCarlo::Tags::RandomNumberGenerator>>(
         make_not_null(&box), std::move(rng));

src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp

@@ -162,6 +162,22 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
     initial_time = packet.time;
     dt_end_step = final_time - initial_time;
 
+    // Find closest grid point to packet at current time, using
+    // extents for live points only.
+    {
+      std::array<size_t, 3> closest_point_index_3d{0, 0, 0};
+      for (size_t d = 0; d < 3; 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] +
+          extents[0] * (closest_point_index_3d[1] +
+                        extents[1] * closest_point_index_3d[2]);
+    }
+
     // Get quantities that we do NOT update if the packet
     // changes cell.
     // local_idx is the index on the mesh without ghost zones

src/ParallelAlgorithms/Events/ObserveTimeStep.hpp

@@ -165,17 +165,18 @@ class ObserveTimeStep : public Event {
   // We obtain the grid size from the variables, rather than the mesh,
   // so that this observer is not DG-specific.
   using argument_tags =
-      tmpl::list<::Tags::TimeStep, typename System::variables_tag>;
+      tmpl::list<::Tags::TimeStep>;  //, typename System::variables_tag>;
 
   template <typename ArrayIndex, typename ParallelComponent,
             typename Metavariables>
   void operator()(const TimeDelta& time_step,
-                  const typename System::variables_tag::type& variables,
+                  // const typename System::variables_tag::type& variables,
                   Parallel::GlobalCache<Metavariables>& cache,
                   const ArrayIndex& array_index,
                   const ParallelComponent* const /*meta*/,
                   const ObservationValue& observation_value) const {
-    const size_t number_of_grid_points = variables.number_of_grid_points();
+    const size_t number_of_grid_points =
+        1;  // variables.number_of_grid_points();
     const double slab_size = time_step.slab().duration().value();
     const double step_size = abs(time_step.value());
     const double wall_time = sys::wall_time();