Optimize Mapping of Random Ray Source Regions to Tallies (#3465)

Author: jtramm

include/openmc/random_ray/flat_source_domain.h

@@ -34,7 +34,7 @@ class FlatSourceDomain {
 
   int64_t add_source_to_scalar_flux();
   virtual void batch_reset();
-  void convert_source_regions_to_tallies();
+  void convert_source_regions_to_tallies(int64_t start_sr_id);
   void reset_tally_volumes();
   void random_ray_tally();
   virtual void accumulate_iteration_flux();

src/random_ray/flat_source_domain.cpp

@@ -421,7 +421,12 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
 // be passed back to the caller to alert them that this function doesn't
 // need to be called for the remainder of the simulation.
 
-void FlatSourceDomain::convert_source_regions_to_tallies()
+// It takes as an argument the starting index in the source region array,
+// and it will operate from that index until the end of the array. This
+// is useful as it can be called for both explicit user source regions or
+// when a source region mesh is overlaid.
+
+void FlatSourceDomain::convert_source_regions_to_tallies(int64_t start_sr_id)
 {
   openmc::simulation::time_tallies.start();
 
@@ -430,7 +435,7 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
 
 // Attempt to generate mapping for all source regions
 #pragma omp parallel for
-  for (int64_t sr = 0; sr < n_source_regions(); sr++) {
+  for (int64_t sr = start_sr_id; sr < n_source_regions(); sr++) {
 
     // If this source region has not been hit by a ray yet, then
     // we aren't going to be able to map it, so skip it.
@@ -468,7 +473,7 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
       // Loop over all active tallies. This logic is essentially identical
       // to what happens when scanning for applicable tallies during
       // MC transport.
-      for (auto i_tally : model::active_tallies) {
+      for (int i_tally = 0; i_tally < model::tallies.size(); i_tally++) {
         Tally& tally {*model::tallies[i_tally]};
 
         // Initialize an iterator over valid filter bin combinations.
@@ -1525,16 +1530,18 @@ void FlatSourceDomain::finalize_discovered_source_regions()
     // order due to shared memory threading.
     std::sort(keys.begin(), keys.end());
 
+    // Remember the index of the first new source region
+    int64_t start_sr_id = source_regions_.n_source_regions();
+
     // Append the source regions in the sorted key order.
     for (const auto& key : keys) {
       const SourceRegion& sr = discovered_source_regions_[key];
       source_region_map_[key] = source_regions_.n_source_regions();
       source_regions_.push_back(sr);
     }
 
-    // If any new source regions were discovered, we need to update the
-    // tally mapping between source regions and tally bins.
-    mapped_all_tallies_ = false;
+    // Map all new source regions to tallies
+    convert_source_regions_to_tallies(start_sr_id);
   }
 
   discovered_source_regions_.clear();

src/random_ray/random_ray_simulation.cpp

@@ -508,8 +508,9 @@ void RandomRaySimulation::simulate()
         domain_->accumulate_iteration_flux();
 
         // Generate mapping between source regions and tallies
-        if (!domain_->mapped_all_tallies_) {
-          domain_->convert_source_regions_to_tallies();
+        if (!domain_->mapped_all_tallies_ &&
+            !RandomRay::mesh_subdivision_enabled_) {
+          domain_->convert_source_regions_to_tallies(0);
         }
 
         // Use above mapping to contribute FSR flux data to appropriate

src/random_ray/source_region.cpp

@@ -259,15 +259,11 @@ void SourceRegionContainer::adjoint_reset()
     MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
   std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(),
     MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
-  for (auto& task_set : volume_task_) {
-    task_set.clear();
-  }
   std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0);
   std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0);
   std::fill(scalar_flux_final_.begin(), scalar_flux_final_.end(), 0.0);
   std::fill(source_.begin(), source_.end(), 0.0f);
   std::fill(external_source_.begin(), external_source_.end(), 0.0f);
-
   std::fill(source_gradients_.begin(), source_gradients_.end(),
     MomentArray {0.0, 0.0, 0.0});
   std::fill(flux_moments_old_.begin(), flux_moments_old_.end(),
@@ -276,10 +272,6 @@ void SourceRegionContainer::adjoint_reset()
     MomentArray {0.0, 0.0, 0.0});
   std::fill(flux_moments_t_.begin(), flux_moments_t_.end(),
     MomentArray {0.0, 0.0, 0.0});
-
-  for (auto& task_set : tally_task_) {
-    task_set.clear();
-  }
 }
 
 } // namespace openmc