Refactor Event Based Tally Related Functions

include/openmc/particle.h

@@ -39,6 +39,30 @@ class Particle : public ParticleData {
 
   double speed() const;
 
+  //! get the distance for particle movement
+  //
+  //! sample a distance for next collision based on the xs
+  //! and find the distance to the boundary along the particle's
+  //! direction and return the above two's minimum
+  double get_destination_distance();
+
+  // Advance particle in space and time
+  // Short-term solution until the surface source is revised and we can use
+  // this->move_distance(distance)
+  void is_hit_time_boundary(double distance, bool& hit_time_boundary);
+
+  void set_particle_weight_to_zero_if_it_hit_time_boundary(bool hit_time_boundary);
+
+  void score_non_mesh_track_length_tallies(double distance);
+
+  void event_advance_deprecated();
+
+  //! score tallies corresponding particle movement
+  //
+  //! score tracklength tally, keff tally and differential tallies
+  //! \param distance the distance the particle is moved
+  void score_the_tallies(double distance);
+
   //! moves the particle by the distance length to its next location
   //! \param length the distance the particle is moved
   void move_distance(double length);

src/event.cpp

@@ -103,15 +103,32 @@ void process_calculate_xs_events(SharedArray<EventQueueItem>& queue)
   simulation::time_event_calculate_xs.stop();
 }
 
+void openmc_event_advance_wrapper()
+{
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
+    int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
+    Particle& p = simulation::particles[buffer_idx];
+    p.event_advance();
+  }
+}
+
+void send_particles_to_other_queues();
+
 void process_advance_particle_events()
 {
   simulation::time_event_advance_particle.start();
+  openmc_event_advance_wrapper();
+  send_particles_to_other_queues();
 
+  simulation::time_event_advance_particle.stop();
+}
+void send_particles_to_other_queues()
+{
 #pragma omp parallel for schedule(runtime)
   for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
     int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
     Particle& p = simulation::particles[buffer_idx];
-    p.event_advance();
     if (!p.alive())
       continue;
     if (p.collision_distance() > p.boundary().distance) {
@@ -122,8 +139,6 @@ void process_advance_particle_events()
   }
 
   simulation::advance_particle_queue.resize(0);
-
-  simulation::time_event_advance_particle.stop();
 }
 
 void process_surface_crossing_events()

src/particle.cpp

@@ -209,7 +209,89 @@ void Particle::event_calculate_xs()
   }
 }
 
+
+double Particle::get_destination_distance()
+{
+  // Find the distance to the nearest boundary
+  this->boundary() = distance_to_boundary(*this);
+
+  // Sample a distance to collision
+  if (this->type() == ParticleType::electron || this->type() == ParticleType::positron) {
+    this->collision_distance() = 0.0;
+  } else if (this->macro_xs().total == 0.0) {
+    this->collision_distance() = INFINITY;
+  } else {
+    this->collision_distance() = -log(prn(this->current_seed())) / this->macro_xs().total;
+  }
+
+  // Select smaller of the two distances
+  double distance = std::min(this->boundary().distance, this->collision_distance());
+  return distance;
+}
+
+void Particle::is_hit_time_boundary(double distance, bool& hit_time_boundary)
+{
+  for (int j = 0; j < this->n_coord(); ++j) {
+    this->coord(j).r += distance * this->coord(j).u;
+  }
+  this->time() += distance / this->speed();
+
+  // Kill particle if its time exceeds the cutoff
+  hit_time_boundary = false;
+  double time_cutoff = settings::time_cutoff[static_cast<int>(this->type())];
+  if (this->time() > time_cutoff) {
+    double dt = this->time() - time_cutoff;
+    this->time() = time_cutoff;
+
+    double push_back_distance = this->speed() * dt;
+    this->move_distance(-push_back_distance);
+    hit_time_boundary = true;
+  }
+}
+
+void Particle::set_particle_weight_to_zero_if_it_hit_time_boundary(bool hit_time_boundary)
+{
+  // Set particle weight to zero if it hit the time boundary
+  if (hit_time_boundary) {
+    this->wgt() = 0.0;
+  }
+}
+
+void Particle::score_the_tallies(double distance)
+{
+  if (!model::active_tracklength_tallies.empty()) {
+    score_tracklength_tally(*this, distance);
+  }
+  // Score track-length estimate of k-eff
+  score_non_mesh_track_length_tallies(distance);
+}
+
+void Particle::score_non_mesh_track_length_tallies(double distance)
+{
+  if (settings::run_mode == RunMode::EIGENVALUE &&
+      this->type() == ParticleType::neutron) {
+    this->keff_tally_tracklength() += this->wgt() * distance * this->macro_xs().nu_fission;
+  }
+
+  // Score flux derivative accumulators for differential tallies.
+  if (!model::active_tallies.empty()) {
+    score_track_derivative(*this, distance);
+  }
+}
+
+
 void Particle::event_advance()
+{
+  double distance = get_destination_distance();
+  bool hit_time_boundary = false;
+  is_hit_time_boundary(distance, hit_time_boundary);
+
+  // Score track-length tallies
+  score_the_tallies(distance);
+  set_particle_weight_to_zero_if_it_hit_time_boundary(hit_time_boundary);
+}
+
+void Particle::event_advance_deprecated()
 {
   // Find the distance to the nearest boundary
   boundary() = distance_to_boundary(*this);

src/tallies/tally_scoring.cpp

@@ -2338,14 +2338,77 @@ void score_analog_tally_mg(Particle& p)
     match.bins_present_ = false;
 }
 
-void score_tracklength_tally(Particle& p, double distance)
+
+double get_atom_density(
+  Particle& p, int i_log_union, const Tally& tally, int i_nuclide)
 {
-  // Determine the tracklength estimate of the flux
-  double flux = p.wgt() * distance;
+  double atom_density = 0.;
+  if (i_nuclide >= 0) {
+    if (p.material() != MATERIAL_VOID) {
+      const auto& mat = model::materials[p.material()];
+      auto j = mat->mat_nuclide_index_[i_nuclide];
+      if (j == C_NONE) {
+        // Determine log union grid index
+        if (i_log_union == C_NONE) {
+          int neutron = static_cast<int>(ParticleType::neutron);
+          i_log_union = log(p.E() / data::energy_min[neutron]) /
+                        simulation::log_spacing;
+        }
 
-  // Set 'none' value for log union grid index
-  int i_log_union = C_NONE;
+        // Update micro xs cache
+        if (!tally.multiply_density()) {
+          p.update_neutron_xs(i_nuclide, i_log_union);
+          atom_density = 1.0;
+        }
+      } else {
+        atom_density =
+          tally.multiply_density() ? mat->atom_density_(j) : 1.0;
+      }
+    }
+  }
+  return atom_density;
+}
+
+void loop_over_filters(Particle& p, double flux, int i_log_union,
+  int i_tally, const Tally& tally)
+{
+  // Initialize an iterator over valid filter bin combinations.  If there are
+  // no valid combinations, use a continue statement to ensure we skip the
+  // assume_separate break below.
+  auto filter_iter = FilterBinIter(tally, p);
+  auto end = FilterBinIter(tally, true, &p.filter_matches());
+  if (filter_iter == end)
+    return;
 
+  // Loop over filter bins.
+  for (; filter_iter != end; ++filter_iter) {
+    auto filter_index = filter_iter.index_;
+    auto filter_weight = filter_iter.weight_;
+
+    // Loop over nuclide bins.
+    for (auto i = 0; i < tally.nuclides_.size(); ++i) {
+      auto i_nuclide = tally.nuclides_[i];
+
+      double atom_density =
+        get_atom_density(p, i_log_union, tally, i_nuclide);
+
+      // TODO: consider replacing this "if" with pointers or templates
+      // handles all estimators except tracklength estimator
+      // tracklength estimator tallies are done in the FilterBinIterator call above
+      if (settings::run_CE) {
+        score_general_ce_nonanalog(p, i_tally, i * tally.scores_.size(),
+          filter_index, filter_weight, i_nuclide, atom_density, flux);
+      } else {
+        score_general_mg(p, i_tally, i * tally.scores_.size(), filter_index,
+          filter_weight, i_nuclide, atom_density, flux);
+      }
+    }
+  }
+}
+
+// compute all the active tracklength tallies
+void loop_over_all_tallies(Particle& p, double flux, int i_log_union)
+{
   for (auto i_tally : model::active_tracklength_tallies) {
     const Tally& tally {*model::tallies[i_tally]};
 
@@ -2357,50 +2420,7 @@ void score_tracklength_tally(Particle& p, double distance)
     if (filter_iter == end)
       continue;
 
-    // Loop over filter bins.
-    for (; filter_iter != end; ++filter_iter) {
-      auto filter_index = filter_iter.index_;
-      auto filter_weight = filter_iter.weight_;
-
-      // Loop over nuclide bins.
-      for (auto i = 0; i < tally.nuclides_.size(); ++i) {
-        auto i_nuclide = tally.nuclides_[i];
-
-        double atom_density = 0.;
-        if (i_nuclide >= 0) {
-          if (p.material() != MATERIAL_VOID) {
-            const auto& mat = model::materials[p.material()];
-            auto j = mat->mat_nuclide_index_[i_nuclide];
-            if (j == C_NONE) {
-              // Determine log union grid index
-              if (i_log_union == C_NONE) {
-                int neutron = static_cast<int>(ParticleType::neutron);
-                i_log_union = std::log(p.E() / data::energy_min[neutron]) /
-                              simulation::log_spacing;
-              }
-
-              // Update micro xs cache
-              if (!tally.multiply_density()) {
-                p.update_neutron_xs(i_nuclide, i_log_union);
-                atom_density = 1.0;
-              }
-            } else {
-              atom_density =
-                tally.multiply_density() ? mat->atom_density_(j) : 1.0;
-            }
-          }
-        }
-
-        // TODO: consider replacing this "if" with pointers or templates
-        if (settings::run_CE) {
-          score_general_ce_nonanalog(p, i_tally, i * tally.scores_.size(),
-            filter_index, filter_weight, i_nuclide, atom_density, flux);
-        } else {
-          score_general_mg(p, i_tally, i * tally.scores_.size(), filter_index,
-            filter_weight, i_nuclide, atom_density, flux);
-        }
-      }
-    }
+    loop_over_filters(p, flux, i_log_union, i_tally, tally);
 
     // If the user has specified that we can assume all tallies are spatially
     // separate, this implies that once a tally has been scored to, we needn't
@@ -2409,6 +2429,16 @@ void score_tracklength_tally(Particle& p, double distance)
     if (settings::assume_separate)
       break;
   }
+}
+
+void score_tracklength_tally(Particle& p, double distance)
+{
+  // Determine the tracklength estimate of the flux
+  double flux = p.wgt() * distance;
+
+  // Set 'none' value for log union grid index
+  int i_log_union = C_NONE;
+  loop_over_all_tallies(p, flux, i_log_union);
 
   // Reset all the filter matches for the next tally event.
   for (auto& match : p.filter_matches())