finished comment incorporation from code review

examples/xml/pincell/settings.xml

@@ -3,9 +3,9 @@
 
   <!-- Define how many particles to run and for how many batches -->
   <run_mode>eigenvalue</run_mode>
-  <batches>20</batches>
+  <batches>100</batches>
   <inactive>10</inactive>
-  <particles>100000</particles>
+  <particles>1000</particles>
 
   <!-- The starting source is a uniform distribution over the entire pin
        cell. Note that since this is effectively a 2D model, the z coordinates

include/openmc/particle.h

@@ -176,9 +176,9 @@ public:
   
   //! Saved ("banked") state of a particle, for nu-fission tallying
   struct NuBank {
-    double E;
-    double wgt;
-    int delayed_group;
+    double E;  //!< particle energy
+    double wgt; //!< particle weight
+    int delayed_group; //!< particle delayed group
   };
 
   //==========================================================================
@@ -343,10 +343,10 @@ public:
   std::vector<NuBank> nu_bank_; // bank of most recently fissioned particles
 
   // Global tally accumulators
-  double tally_absorption_ {0.0};
-  double tally_collision_ {0.0};
-  double tally_tracklength_ {0.0};
-  double tally_leakage_ {0.0};
+  double keff_tally_absorption_ {0.0};
+  double keff_tally_collision_ {0.0};
+  double keff_tally_tracklength_ {0.0};
+  double keff_tally_leakage_ {0.0};
 
   bool trace_ {false};     //!< flag to show debug information
 

include/openmc/simulation.h

@@ -63,8 +63,12 @@ void initialize_batch();
 //! Initialize a fission generation
 void initialize_generation();
 
+//! Full initialization of a particle history
 void initialize_history(Particle* p, int64_t index_source);
 
+//! Helper function for initialize_history() that is called independently elsewhere
+void initialize_history_partial(Particle* p);
+
 //! Finalize a batch
 //!
 //! Handles synchronization and accumulation of tallies, calculation of Shannon

src/eigenvalue.cpp

@@ -138,7 +138,7 @@ void synchronize_bank()
   int64_t index_temp = 0;
   std::vector<Particle::Bank> temp_sites(3*simulation::work_per_rank);
 
-  for (uint64_t i = 0; i < simulation::fission_bank_length; i++ ) {
+  for (int64_t i = 0; i < simulation::fission_bank_length; i++ ) {
     const auto& site = simulation::fission_bank[i];
 
     // If there are less than n_particles particles banked, automatically add

src/particle.cpp

@@ -132,67 +132,67 @@ Particle::from_source(const Bank* src)
 void
 Particle::event_calculate_xs()
 {
-    // Set the random number stream
-    if (type_ == Particle::Type::neutron) {
-      stream_ = STREAM_TRACKING;
-    } else {
-      stream_ = STREAM_PHOTON;
+  // Set the random number stream
+  if (type_ == Particle::Type::neutron) {
+    stream_ = STREAM_TRACKING;
+  } else {
+    stream_ = STREAM_PHOTON;
+  }
+
+  // Store pre-collision particle properties
+  wgt_last_ = wgt_;
+  E_last_ = E_;
+  u_last_ = this->u();
+  r_last_ = this->r();
+
+  // Reset event variables
+  event_ = EVENT_KILL;
+  event_nuclide_ = NUCLIDE_NONE;
+  event_mt_ = REACTION_NONE;
+
+  // If the cell hasn't been determined based on the particle's location,
+  // initiate a search for the current cell. This generally happens at the
+  // beginning of the history and again for any secondary particles
+  if (coord_[n_coord_ - 1].cell == C_NONE) {
+    if (!find_cell(this, false)) {
+      this->mark_as_lost("Could not find the cell containing particle "
+        + std::to_string(id_));
+      return;
     }
 
-    // Store pre-collision particle properties
-    wgt_last_ = wgt_;
-    E_last_ = E_;
-    u_last_ = this->u();
-    r_last_ = this->r();
+    // Set birth cell attribute
+    if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell;
+  }
 
-    // Reset event variables
-    event_ = EVENT_KILL;
-    event_nuclide_ = NUCLIDE_NONE;
-    event_mt_ = REACTION_NONE;
+  // Write particle track.
+  if (write_track_) write_particle_track(*this);
 
-    // If the cell hasn't been determined based on the particle's location,
-    // initiate a search for the current cell. This generally happens at the
-    // beginning of the history and again for any secondary particles
-    if (coord_[n_coord_ - 1].cell == C_NONE) {
-      if (!find_cell(this, false)) {
-        this->mark_as_lost("Could not find the cell containing particle "
-          + std::to_string(id_));
-        return;
-      }
-
-      // Set birth cell attribute
-      if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell;
-    }
-
-    // Write particle track.
-    if (write_track_) write_particle_track(*this);
-
-    if (settings::check_overlaps) check_cell_overlap(this);
-    
-    // Calculate microscopic and macroscopic cross sections
-    if (material_ != MATERIAL_VOID) {
-      if (settings::run_CE) {
-        if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
-          // If the material is the same as the last material and the
-          // temperature hasn't changed, we don't need to lookup cross
-          // sections again.
-          model::materials[material_]->calculate_xs(*this);
-        }
-      } else {
-        // Get the MG data; unlike the CE case above, we have to re-calculate
-        // cross sections for every collision since the cross sections may
-        // be angle-dependent
-        data::mg.macro_xs_[material_].calculate_xs(*this);
-
-        // Update the particle's group while we know we are multi-group
-        g_last_ = g_;
+  if (settings::check_overlaps) check_cell_overlap(this);
+  
+  // Calculate microscopic and macroscopic cross sections
+  if (material_ != MATERIAL_VOID) {
+    if (settings::run_CE) {
+      if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
+        // If the material is the same as the last material and the
+        // temperature hasn't changed, we don't need to lookup cross
+        // sections again.
+        model::materials[material_]->calculate_xs(*this);
       }
     } else {
-      macro_xs_.total      = 0.0;
-      macro_xs_.absorption = 0.0;
-      macro_xs_.fission    = 0.0;
-      macro_xs_.nu_fission = 0.0;
+      // Get the MG data; unlike the CE case above, we have to re-calculate
+      // cross sections for every collision since the cross sections may
+      // be angle-dependent
+      data::mg.macro_xs_[material_].calculate_xs(*this);
+
+      // Update the particle's group while we know we are multi-group
+      g_last_ = g_;
     }
+  } else {
+    macro_xs_.total      = 0.0;
+    macro_xs_.absorption = 0.0;
+    macro_xs_.fission    = 0.0;
+    macro_xs_.nu_fission = 0.0;
+  }
 }
 
 void
@@ -227,7 +227,7 @@ Particle::event_advance()
   // Score track-length estimate of k-eff
   if (settings::run_mode == RUN_MODE_EIGENVALUE &&
       type_ == Particle::Type::neutron) {
-    tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission;
+    keff_tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission;
   }
 
   // Score flux derivative accumulators for differential tallies.
@@ -272,7 +272,7 @@ Particle::event_collide()
   // Score collision estimate of keff
   if (settings::run_mode == RUN_MODE_EIGENVALUE &&
       type_ == Particle::Type::neutron) {
-    tally_collision_ += wgt_ * macro_xs_.nu_fission
+    keff_tally_collision_ += wgt_ * macro_xs_.nu_fission
       / macro_xs_.total;
   }
 
@@ -378,19 +378,19 @@ Particle::event_death()
 
   // Contribute tally reduction variables to global accumulator
   #pragma omp atomic
-  global_tally_absorption += tally_absorption_;
+  global_tally_absorption += keff_tally_absorption_;
   #pragma omp atomic
-  global_tally_collision += tally_collision_;
+  global_tally_collision += keff_tally_collision_;
   #pragma omp atomic
-  global_tally_tracklength += tally_tracklength_;
+  global_tally_tracklength += keff_tally_tracklength_;
   #pragma omp atomic
-  global_tally_leakage += tally_leakage_;
+  global_tally_leakage += keff_tally_leakage_;
 
   // Reset particle tallies once accumulated
-  tally_absorption_  = 0.0;
-  tally_collision_   = 0.0;
-  tally_tracklength_ = 0.0;
-  tally_leakage_     = 0.0;
+  keff_tally_absorption_  = 0.0;
+  keff_tally_collision_   = 0.0;
+  keff_tally_tracklength_ = 0.0;
+  keff_tally_leakage_     = 0.0;
 }
 
 
@@ -424,7 +424,7 @@ Particle::cross_surface()
     }
 
     // Score to global leakage tally
-    tally_leakage_ += wgt_;
+    keff_tally_leakage_ += wgt_;
 
     // Display message
     if (settings::verbosity >= 10 || trace_) {

src/particle_restart.cpp

@@ -111,8 +111,7 @@ void run_particle_restart()
   if (p.write_track_) add_particle_track(p);
 
   // Every particle starts with no accumulated flux derivative.
-  if (!model::active_tallies.empty())
-  {
+  if (!model::active_tallies.empty()) {
     p.flux_derivs_.resize(model::tally_derivs.size(), 0.0);
     std::fill(p.flux_derivs_.begin(), p.flux_derivs_.end(), 0.0);
   }

src/physics.cpp

@@ -176,7 +176,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
 
   // Determine whether to place fission sites into the shared fission bank
   // or the secondary particle bank.
-  bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE)
+  bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE);
 
   for (int i = 0; i < nu; ++i) {
     Particle::Bank* site;
@@ -592,7 +592,7 @@ void absorption(Particle* p, int i_nuclide)
 
     // Score implicit absorption estimate of keff
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      p->tally_absorption_ += p->wgt_absorb_ * p->neutron_xs_[
+      p->keff_tally_absorption_ += p->wgt_absorb_ * p->neutron_xs_[
         i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
     }
   } else {
@@ -601,7 +601,7 @@ void absorption(Particle* p, int i_nuclide)
         prn(p->current_seed()) * p->neutron_xs_[i_nuclide].total) {
       // Score absorption estimate of keff
       if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-        p->tally_absorption_ += p->wgt_ * p->neutron_xs_[
+        p->keff_tally_absorption_ += p->wgt_ * p->neutron_xs_[
           i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
       }
 

src/physics_mg.cpp

@@ -122,7 +122,7 @@ create_fission_sites(Particle* p)
   
   // Determine whether to place fission sites into the shared fission bank
   // or the secondary particle bank.
-  bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE)
+  bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE);
 
   for (int i = 0; i < nu; ++i) {
     Particle::Bank* site;
@@ -222,11 +222,11 @@ absorption(Particle* p)
     p->wgt_last_ = p->wgt_;
 
     // Score implicit absorpion estimate of keff
-    p->tally_absorption_ += p->wgt_absorb_ * p->macro_xs_.nu_fission /
+    p->keff_tally_absorption_ += p->wgt_absorb_ * p->macro_xs_.nu_fission /
         p->macro_xs_.absorption;
   } else {
     if (p->macro_xs_.absorption > prn(p->current_seed()) * p->macro_xs_.total) {
-      p->tally_absorption_ += p->wgt_ * p->macro_xs_.nu_fission /
+      p->keff_tally_absorption_ += p->wgt_ * p->macro_xs_.nu_fission /
            p->macro_xs_.absorption;
       p->alive_ = false;
       p->event_ = EVENT_ABSORB;

src/simulation.cpp

@@ -469,6 +469,11 @@ void initialize_history(Particle* p, int64_t index_source)
   #pragma omp atomic
   simulation::total_weight += p->wgt_;
 
+  initialize_history_partial(p);
+}
+
+void initialize_history_partial(Particle* p)
+{
   // Force calculation of cross-sections by setting last energy to zero
   if (settings::run_CE) {
     for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;

src/tallies/tally_scoring.cpp

@@ -55,7 +55,8 @@ FilterBinIter::FilterBinIter(const Tally& tally, Particle* p)
   this->compute_index_weight();
 }
 
-FilterBinIter::FilterBinIter(const Tally& tally, bool end, std::vector<FilterMatch>* particle_filter_matches)
+FilterBinIter::FilterBinIter(const Tally& tally, bool end,
+    std::vector<FilterMatch>* particle_filter_matches)
   : tally_{tally}, filter_matches_{*particle_filter_matches}
 {
   // Handle the special case for an iterator that points to the end.
@@ -144,7 +145,8 @@ FilterBinIter::compute_index_weight()
 //! Helper function used to increment tallies with a delayed group filter.
 
 void
-score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index, std::vector<FilterMatch>& filter_matches)
+score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index,
+    std::vector<FilterMatch>& filter_matches)
 {
   // Save the original delayed group bin
   auto& tally {*model::tallies[i_tally]};