Fix analog fission tallies for sites rejected by the time cutoff

docs/source/io_formats/settings.rst

@@ -206,7 +206,9 @@ time.
     *Default*: 0.0
 
   :time_neutron
-    The time above which neutrons will be killed.
+    The time above which neutrons will be killed. Delayed fission neutrons that
+    would be emitted after this time are not banked, but analog fission tallies
+    count them as produced, as the track-length and collision estimators do.
 
     *Default*: Infinity
 

include/openmc/particle_data.h

@@ -669,16 +669,19 @@ class ParticleData : public GeometryState {
   } // energy of last reaction secondaries
   const double& bank_second_E() const { return bank_second_E_; }
 
-  int& n_bank() { return n_bank_; }        // number of banked fission sites
-  double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites
+  // Fission neutrons recorded in nu_bank() for analog fission tallies: those
+  // stored in the fission or secondary bank and delayed neutrons that are not
+  // banked because they would be emitted after the time cutoff
+  int& n_bank() { return n_bank_; }        // number of fission neutrons
+  double& wgt_bank() { return wgt_bank_; } // weight of fission neutrons
   int* n_delayed_bank()
   {
     return n_delayed_bank_;
-  } // number of delayed fission sites
+  } // number of delayed fission neutrons
   int& n_delayed_bank(int i)
   {
     return n_delayed_bank_[i];
-  } // number of delayed fission sites
+  } // number of delayed fission neutrons
 
   // Index of cell particle is born in
   int& cell_born() { return cell_born_; }

openmc/settings.py

@@ -84,7 +84,10 @@ class Settings:
         assigned to particles that are not killed after Russian roulette. Value
         of energy should be a float indicating energy in eV below which particle
         type will be killed. Value of time should be a float in seconds.
-        Particles will be killed exactly at the specified time. Value for
+        Particles will be killed exactly at the specified time. Delayed fission
+        neutrons that would be emitted after the 'time_neutron' cutoff are not
+        banked, but analog fission tallies count them as produced, as the
+        track-length and collision estimators do. Value for
         'survival_normalization' is a bool indicating whether or not the weight
         cutoff parameters will be applied relative to the particle's starting
         weight or to its current weight.

src/physics.cpp

@@ -204,11 +204,21 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
   // or the secondary particle bank.
   bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
 
-  // Counter for the number of fission sites successfully stored to the shared
-  // fission bank or the secondary particle bank
-  int n_sites_stored;
+  // Record a fission neutron for analog fission tallies
+  auto record_fission_neutron = [&](const SourceSite& site) {
+    // Increment the number of neutrons born delayed
+    if (site.delayed_group > 0) {
+      nu_d[site.delayed_group - 1]++;
+    }
+
+    // Write fission particles to nuBank
+    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
+    nu_bank_entry.wgt = site.wgt;
+    nu_bank_entry.E = site.E;
+    nu_bank_entry.delayed_group = site.delayed_group;
+  };
 
-  for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
+  for (int i = 0; i < nu; ++i) {
     // Initialize fission site object with particle data
     SourceSite site;
     site.r = p.r();
@@ -220,10 +230,13 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, &site, p);
 
-    // Reject site if it exceeds time cutoff
+    // Reject site if it exceeds time cutoff. The neutron would be killed at
+    // birth, so it is not banked, but it was produced by this fission and is
+    // recorded for analog fission tallies like any other fission neutron.
     if (site.delayed_group > 0) {
       double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
       if (site.time > t_cutoff) {
+        record_fission_neutron(site);
         continue;
       }
     }
@@ -261,30 +274,22 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
       p.n_secondaries()++;
     }
 
-    // Increment the number of neutrons born delayed
-    if (site.delayed_group > 0) {
-      nu_d[site.delayed_group - 1]++;
-    }
-
-    // Write fission particles to nuBank
-    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
-    nu_bank_entry.wgt = site.wgt;
-    nu_bank_entry.E = site.E;
-    nu_bank_entry.delayed_group = site.delayed_group;
+    record_fission_neutron(site);
   }
 
-  // If shared fission bank was full, and no fissions could be added,
-  // set the particle fission flag to false.
-  if (n_sites_stored == 0) {
+  // If shared fission bank was full, and no fission neutrons could be
+  // recorded, set the particle fission flag to false.
+  if (p.nu_bank().empty()) {
     p.fission() = false;
     return;
   }
 
-  // Set nu to the number of fission sites successfully stored. If the fission
-  // bank was not found to be full then these values are already equivalent.
-  nu = n_sites_stored;
+  // Set nu to the number of fission neutrons recorded, i.e., the sites stored
+  // in a bank plus those rejected by the time cutoff. If the fission bank was
+  // not found to be full then these values are already equivalent.
+  nu = p.nu_bank().size();
 
-  // Store the total weight banked for analog fission tallies
+  // Store the total weight recorded for analog fission tallies
   p.n_bank() = nu;
   p.wgt_bank() = nu / weight;
   for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {

src/physics_mg.cpp

@@ -129,11 +129,21 @@ void create_fission_sites(Particle& p)
   // or the secondary particle bank.
   bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
 
-  // Counter for the number of fission sites successfully stored to the shared
-  // fission bank or the secondary particle bank
-  int n_sites_stored;
+  // Record a fission neutron for analog fission tallies
+  auto record_fission_neutron = [&](const SourceSite& site) {
+    // Increment the number of neutrons born delayed
+    if (site.delayed_group > 0) {
+      nu_d[site.delayed_group - 1]++;
+    }
+
+    // Write fission particles to nuBank
+    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
+    nu_bank_entry.wgt = site.wgt;
+    nu_bank_entry.E = site.E;
+    nu_bank_entry.delayed_group = site.delayed_group;
+  };
 
-  for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
+  for (int i = 0; i < nu; ++i) {
     // Initialize fission site object with particle data
     SourceSite site;
     site.r = p.r();
@@ -171,9 +181,12 @@ void create_fission_sites(Particle& p)
         macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0);
       site.time -= std::log(prn(p.current_seed())) / decay_rate;
 
-      // Reject site if it exceeds time cutoff
+      // Reject site if it exceeds time cutoff. The neutron would be killed at
+      // birth, so it is not banked, but it was produced by this fission and
+      // is recorded for analog fission tallies like any other fission neutron.
       double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
       if (site.time > t_cutoff) {
+        record_fission_neutron(site);
         continue;
       }
     }
@@ -210,30 +223,22 @@ void create_fission_sites(Particle& p)
     // Set the delayed group on the particle as well
     p.delayed_group() = dg + 1;
 
-    // Increment the number of neutrons born delayed
-    if (p.delayed_group() > 0) {
-      nu_d[dg]++;
-    }
-
-    // Write fission particles to nuBank
-    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
-    nu_bank_entry.wgt = site.wgt;
-    nu_bank_entry.E = site.E;
-    nu_bank_entry.delayed_group = site.delayed_group;
+    record_fission_neutron(site);
   }
 
-  // If shared fission bank was full, and no fissions could be added,
-  // set the particle fission flag to false.
-  if (n_sites_stored == 0) {
+  // If shared fission bank was full, and no fission neutrons could be
+  // recorded, set the particle fission flag to false.
+  if (p.nu_bank().empty()) {
     p.fission() = false;
     return;
   }
 
-  // Set nu to the number of fission sites successfully stored. If the fission
-  // bank was not found to be full then these values are already equivalent.
-  nu = n_sites_stored;
+  // Set nu to the number of fission neutrons recorded, i.e., the sites stored
+  // in a bank plus those rejected by the time cutoff. If the fission bank was
+  // not found to be full then these values are already equivalent.
+  nu = p.nu_bank().size();
 
-  // Store the total weight banked for analog fission tallies
+  // Store the total weight recorded for analog fission tallies
   p.n_bank() = nu;
   p.wgt_bank() = nu / weight;
   for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {