Executable for MC: initialize background

src/DataStructures/Variables.hpp

@@ -594,6 +594,8 @@ class Variables<tmpl::list<>> {
   template <typename T>
   Variables(const T* /*pointer*/, const size_t /*size*/) {}
   static constexpr size_t size() { return 0; }
+  void assign_subset(const Variables<tmpl::list<>>& /*unused*/) {}
+  void assign_subset(const tuples::TaggedTuple<>& /*unused*/) {}
 };
 
 // gcc8 screams when the empty Variables has pup as a member function, so we

src/Evolution/DgSubcell/Actions/Initialize.hpp

@@ -46,6 +46,7 @@
 #include "Utilities/CallWithDynamicType.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/Overloader.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
 
@@ -187,17 +188,10 @@ struct SetSubcellGrid {
         },
         make_not_null(&box));
 
-    db::mutate_apply<
-        tmpl::list<Tags::ActiveGrid, Tags::DidRollback,
-                   typename System::variables_tag, subcell::Tags::TciDecision,
-                   subcell::Tags::TciCallsSinceRollback,
-                   subcell::Tags::StepsSinceTciCall>,
-        tmpl::list<>>(
-        [&cell_is_not_on_external_boundary, &dg_mesh,
-         subcell_allowed_in_element, &subcell_mesh](
+    const auto init_non_vars =
+        [&cell_is_not_on_external_boundary, subcell_allowed_in_element](
             const gsl::not_null<ActiveGrid*> active_grid_ptr,
             const gsl::not_null<bool*> did_rollback_ptr,
-            const auto active_vars_ptr,
             const gsl::not_null<int*> tci_decision_ptr,
             const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
             const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
@@ -210,17 +204,59 @@ struct SetSubcellGrid {
                subcell_enabled_at_external_boundary) and
               subcell_allowed_in_element) {
             *active_grid_ptr = ActiveGrid::Subcell;
-            active_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
           } else {
             *active_grid_ptr = ActiveGrid::Dg;
-            active_vars_ptr->initialize(dg_mesh.number_of_grid_points());
           }
 
           *tci_decision_ptr = 0;
           *tci_calls_since_rollback_ptr = 0;
           *steps_since_tci_call_ptr = 0;
+        };
+
+    db::mutate_apply<
+        tmpl::flatten<tmpl::list<
+            Tags::ActiveGrid, Tags::DidRollback,
+            tmpl::conditional_t<
+                std::is_same_v<typename System::variables_tag::tags_list,
+                               tmpl::list<>>,
+                tmpl::list<>, typename System::variables_tag>,
+            subcell::Tags::TciDecision, subcell::Tags::TciCallsSinceRollback,
+            subcell::Tags::StepsSinceTciCall>>,
+        tmpl::list<>>(
+        Overloader{
+            [&init_non_vars](
+                const gsl::not_null<ActiveGrid*> active_grid_ptr,
+                const gsl::not_null<bool*> did_rollback_ptr,
+                const gsl::not_null<int*> tci_decision_ptr,
+                const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
+                const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
+              init_non_vars(active_grid_ptr, did_rollback_ptr, tci_decision_ptr,
+                            tci_calls_since_rollback_ptr,
+                            steps_since_tci_call_ptr);
+            },
+            [&init_non_vars, &dg_mesh, &subcell_mesh,
+             &cell_is_not_on_external_boundary, subcell_allowed_in_element](
+                const gsl::not_null<ActiveGrid*> active_grid_ptr,
+                const gsl::not_null<bool*> did_rollback_ptr,
+                const auto active_vars_ptr,
+                const gsl::not_null<int*> tci_decision_ptr,
+                const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
+                const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
+              init_non_vars(active_grid_ptr, did_rollback_ptr, tci_decision_ptr,
+                            tci_calls_since_rollback_ptr,
+                            steps_since_tci_call_ptr);
+              if ((cell_is_not_on_external_boundary or
+                   subcell_enabled_at_external_boundary) and
+                  subcell_allowed_in_element) {
+                active_vars_ptr->initialize(
+                    subcell_mesh.number_of_grid_points());
+              } else {
+                active_vars_ptr->initialize(dg_mesh.number_of_grid_points());
+              }
+            },
         },
         make_not_null(&box));
+
     if constexpr (System::has_primitive_and_conservative_vars) {
       db::mutate<typename System::primitive_variables_tag>(
           [&dg_mesh, &subcell_mesh](const auto prim_vars_ptr,

src/Evolution/DgSubcell/BackgroundGrVars.hpp

@@ -105,50 +105,62 @@ struct BackgroundGrVars : tt::ConformsTo<db::protocols::Mutator> {
       const bool did_rollback, const T& solution_or_data) {
     const size_t num_subcell_pts = subcell_mesh.number_of_grid_points();
 
-    if (gsl::at(*subcell_face_gr_vars, 0).number_of_grid_points() != 0) {
-      // Evolution phase
+    bool evolution_phase = false;
+    if constexpr (not std::is_same_v<
+                      typename SubcellFaceGrVars::value_type::tags_list,
+                      tmpl::list<>>) {
+      evolution_phase =
+          (gsl::at(*subcell_face_gr_vars, 0).number_of_grid_points() != 0);
+    }
 
-      // Check if the mesh is actually moving i.e. block coordinate map is
-      // time-dependent. If not, we can skip the evaluation of GR variables
-      // since they may stay with their values assigned at the initialization
-      // phase.
-      const auto& element_id = element.id();
-      const size_t block_id = element_id.block_id();
-      const Block<volume_dim>& block = domain.blocks()[block_id];
+    if (evolution_phase) {
+      if constexpr (not std::is_same_v<
+                        typename SubcellFaceGrVars::value_type::tags_list,
+                        tmpl::list<>>) {
+        // Evolution phase
 
-      if (block.is_time_dependent()) {
-        if (did_rollback or not ComputeOnlyOnRollback) {
-          if (did_rollback) {
-            // Right after rollback, subcell GR vars are stored in the
-            // `inactive` one.
-            ASSERT(inactive_gr_vars->number_of_grid_points() == num_subcell_pts,
-                   "The size of subcell GR variables ("
-                       << inactive_gr_vars->number_of_grid_points()
-                       << ") is not equal to the number of FD grid points ("
-                       << subcell_mesh.number_of_grid_points() << ").");
+        // Check if the mesh is actually moving i.e. block coordinate map is
+        // time-dependent. If not, we can skip the evaluation of GR variables
+        // since they may stay with their values assigned at the initialization
+        // phase.
+        const auto& element_id = element.id();
+        const size_t block_id = element_id.block_id();
+        const Block<volume_dim>& block = domain.blocks()[block_id];
 
-            cell_centered_impl(inactive_gr_vars, time, subcell_inertial_coords,
-                               solution_or_data);
+        if (block.is_time_dependent()) {
+          if (did_rollback or not ComputeOnlyOnRollback) {
+            if (did_rollback) {
+              // Right after rollback, subcell GR vars are stored in the
+              // `inactive` one.
+              ASSERT(
+                  inactive_gr_vars->number_of_grid_points() == num_subcell_pts,
+                  "The size of subcell GR variables ("
+                      << inactive_gr_vars->number_of_grid_points()
+                      << ") is not equal to the number of FD grid points ("
+                      << subcell_mesh.number_of_grid_points() << ").");
 
-          } else {
-            // In this case the element didn't rollback but started from FD.
-            // Therefore subcell GR vars are in the `active` one.
-            ASSERT(active_gr_vars->number_of_grid_points() == num_subcell_pts,
-                   "The size of subcell GR variables ("
-                       << active_gr_vars->number_of_grid_points()
-                       << ") is not equal to the number of FD grid points ("
-                       << subcell_mesh.number_of_grid_points() << ").");
+              cell_centered_impl(inactive_gr_vars, time,
+                                 subcell_inertial_coords, solution_or_data);
 
-            cell_centered_impl(active_gr_vars, time, subcell_inertial_coords,
-                               solution_or_data);
-          }
+            } else {
+              // In this case the element didn't rollback but started from FD.
+              // Therefore subcell GR vars are in the `active` one.
+              ASSERT(active_gr_vars->number_of_grid_points() == num_subcell_pts,
+                     "The size of subcell GR variables ("
+                         << active_gr_vars->number_of_grid_points()
+                         << ") is not equal to the number of FD grid points ("
+                         << subcell_mesh.number_of_grid_points() << ").");
 
-          face_centered_impl(subcell_face_gr_vars, time, functions_of_time,
-                             logical_to_grid_map, grid_to_inertial_map,
-                             subcell_mesh, solution_or_data);
+              cell_centered_impl(active_gr_vars, time, subcell_inertial_coords,
+                                 solution_or_data);
+            }
+
+            face_centered_impl(subcell_face_gr_vars, time, functions_of_time,
+                               logical_to_grid_map, grid_to_inertial_map,
+                               subcell_mesh, solution_or_data);
+          }
         }
       }
-
     } else {
       // Initialization phase
       (*inactive_gr_vars).initialize(num_subcell_pts);
@@ -158,19 +170,22 @@ struct BackgroundGrVars : tt::ConformsTo<db::protocols::Mutator> {
              "The subcell mesh must have isotropic basis, quadrature. and "
              "extents but got "
                  << subcell_mesh);
-      const size_t num_face_centered_mesh_grid_pts =
-          (subcell_mesh.extents(0) + 1) * subcell_mesh.extents(1) *
-          subcell_mesh.extents(2);
-      for (size_t d = 0; d < volume_dim; ++d) {
-        gsl::at(*subcell_face_gr_vars, d)
-            .initialize(num_face_centered_mesh_grid_pts);
+      if constexpr (not std::is_same_v<
+                        typename SubcellFaceGrVars::value_type::tags_list,
+                        tmpl::list<>>) {
+        const size_t num_face_centered_mesh_grid_pts =
+            (subcell_mesh.extents(0) + 1) * subcell_mesh.extents(1) *
+            subcell_mesh.extents(2);
+        for (size_t d = 0; d < volume_dim; ++d) {
+          gsl::at(*subcell_face_gr_vars, d)
+              .initialize(num_face_centered_mesh_grid_pts);
+        }
+        face_centered_impl(subcell_face_gr_vars, time, functions_of_time,
+                           logical_to_grid_map, grid_to_inertial_map,
+                           subcell_mesh, solution_or_data);
       }
-
       cell_centered_impl(inactive_gr_vars, time, subcell_inertial_coords,
                          solution_or_data);
-      face_centered_impl(subcell_face_gr_vars, time, functions_of_time,
-                         logical_to_grid_map, grid_to_inertial_map,
-                         subcell_mesh, solution_or_data);
     }
   }
 

src/Evolution/Executables/RadiationTransport/CMakeLists.txt

@@ -2,3 +2,4 @@
 # See LICENSE.txt for details.
 
 add_subdirectory(M1Grey)
+add_subdirectory(MonteCarlo)

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

@@ -0,0 +1,45 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(EXECUTABLE EvolveMonteCarlo)
+
+add_spectre_executable(
+  ${EXECUTABLE}
+  EXCLUDE_FROM_ALL
+  EvolveMonteCarlo.cpp
+  )
+
+target_link_libraries(
+  ${EXECUTABLE}
+  PRIVATE
+  Actions
+  Charmxx::main
+  CoordinateMaps
+  DataStructures
+  DgSubcell
+  DiscontinuousGalerkin
+  DomainCreators
+  Events
+  EventsAndDenseTriggers
+  EventsAndTriggers
+  Evolution
+  GeneralRelativity
+  GeneralRelativitySolutions
+  GrMhdAnalyticData
+  GrMhdSolutions
+  H5
+  Hydro
+  Informer
+  LinearOperators
+  MathFunctions
+  MonteCarlo
+  Observer
+  Options
+  Parallel
+  PhaseControl
+  Serialization
+  Time
+  Utilities
+  ValenciaDivClean
+  )
+

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

@@ -0,0 +1,33 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp"
+
+#include <vector>
+
+#include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
+#include "Parallel/CharmMain.tpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp"
+#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
+
+void register_neutrino_tables() {
+  register_classes_with_charm(
+      tmpl::list<Particles::MonteCarlo::NeutrinoInteractionTable<2, 2>,
+                 Particles::MonteCarlo::NeutrinoInteractionTable<2, 3>,
+                 Particles::MonteCarlo::NeutrinoInteractionTable<4, 3>,
+                 Particles::MonteCarlo::NeutrinoInteractionTable<16, 3>>{});
+}
+
+extern "C" void CkRegisterMainModule() {
+  Parallel::charmxx::register_main_module<EvolutionMetavars<4, 3>>();
+  Parallel::charmxx::register_init_node_and_proc(
+      {&domain::creators::register_derived_with_charm,
+       &domain::creators::time_dependence::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
+       &EquationsOfState::register_derived_with_charm,
+       &register_factory_classes_with_charm<EvolutionMetavars<4, 3>>,
+       &register_neutrino_tables},
+      {});
+}