Combing for fission site sampling, and delayed neutron emission time

docs/source/methods/neutron_physics.rst

@@ -290,7 +290,10 @@ create and store fission sites for the following generation. First, the average
 number of prompt and delayed neutrons must be determined to decide whether the
 secondary neutrons will be prompt or delayed. This is important because delayed
 neutrons have a markedly different spectrum from prompt neutrons, one that has a
-lower average energy of emission. The total number of neutrons emitted
+lower average energy of emission. Furthermore, in simulations where tracking
+time of neutrons is important, we need to consider the emission time delay of
+the secondary neutrons, which is dependent on the decay constant of the
+delayed neutron precursor. The total number of neutrons emitted
 :math:`\nu_t` is given as a function of incident energy in the ENDF format. Two
 representations exist for :math:`\nu_t`. The first is a polynomial of order
 :math:`N` with coefficients :math:`c_0,c_1,\dots,c_N`. If :math:`\nu_t` has this
@@ -306,8 +309,8 @@ interpolation law. The number of prompt neutrons released per fission event
 :math:`\nu_p` is also given as a function of incident energy and can be
 specified in a polynomial or tabular format. The number of delayed neutrons
 released per fission event :math:`\nu_d` can only be specified in a tabular
-format. In practice, we only need to determine :math:`nu_t` and
-:math:`nu_d`. Once these have been determined, we can calculated the delayed
+format. In practice, we only need to determine :math:`\nu_t` and
+:math:`\nu_d`. Once these have been determined, we can calculate the delayed
 neutron fraction
 
 .. math::
@@ -335,8 +338,14 @@ neutrons. Otherwise, we produce :math:`\lfloor \nu \rfloor + 1` neutrons. Then,
 for each fission site produced, we sample the outgoing angle and energy
 according to the algorithms given in :ref:`sample-angle` and
 :ref:`sample-energy` respectively. If the neutron is to be born delayed, then
-there is an extra step of sampling a delayed neutron precursor group since they
-each have an associated secondary energy distribution.
+there is an extra step of sampling a delayed neutron precursor group to get the
+associated secondary energy distribution and the decay constant
+:math:`\lambda`, which is needed to sample the emission delay time :math:`t_d`:
+
+.. math::
+    :label: sample-delay-time
+
+    t_d = -\frac{\ln \xi}{\lambda}.
 
 The sampled outgoing angle and energy of fission neutrons along with the
 position of the collision site are stored in an array called the fission

examples/pincell_pulsed/run_pulse.py

@@ -0,0 +1,101 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import openmc
+
+###############################################################################
+# Create materials for the problem
+
+uo2 = openmc.Material(name="UO2 fuel at 2.4% wt enrichment")
+uo2.set_density("g/cm3", 10.29769)
+uo2.add_element("U", 1.0, enrichment=2.4)
+uo2.add_element("O", 2.0)
+
+helium = openmc.Material(name="Helium for gap")
+helium.set_density("g/cm3", 0.001598)
+helium.add_element("He", 2.4044e-4)
+
+zircaloy = openmc.Material(name="Zircaloy 4")
+zircaloy.set_density("g/cm3", 6.55)
+zircaloy.add_element("Sn", 0.014, "wo")
+zircaloy.add_element("Fe", 0.00165, "wo")
+zircaloy.add_element("Cr", 0.001, "wo")
+zircaloy.add_element("Zr", 0.98335, "wo")
+
+borated_water = openmc.Material(name="Borated water")
+borated_water.set_density("g/cm3", 0.740582)
+borated_water.add_element("B", 2.0e-4)  # 3x the original pincell
+borated_water.add_element("H", 5.0e-2)
+borated_water.add_element("O", 2.4e-2)
+borated_water.add_s_alpha_beta("c_H_in_H2O")
+
+###############################################################################
+# Define problem geometry
+
+# Create cylindrical surfaces
+fuel_or = openmc.ZCylinder(r=0.39218, name="Fuel OR")
+clad_ir = openmc.ZCylinder(r=0.40005, name="Clad IR")
+clad_or = openmc.ZCylinder(r=0.45720, name="Clad OR")
+
+# Create a region represented as the inside of a rectangular prism
+pitch = 1.25984
+box = openmc.model.RectangularPrism(pitch, pitch, boundary_type="reflective")
+
+# Create cells, mapping materials to regions
+fuel = openmc.Cell(fill=uo2, region=-fuel_or)
+gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
+clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
+water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
+
+# Create a model and assign geometry
+model = openmc.Model()
+model.geometry = openmc.Geometry([fuel, gap, clad, water])
+
+###############################################################################
+# Define problem settings
+
+# Set the mode
+model.settings.run_mode = "fixed source"
+
+# Indicate how many batches and particles to run
+model.settings.batches = 10
+model.settings.particles = 10000
+
+# Set time cutoff (we only care about t < 100 seconds, see tally below)
+model.settings.cutoff = {"time_neutron": 100}
+
+# Create the neutron pulse source (by default, isotropic direction, t=0)
+space = openmc.stats.Point()  # At the origin (0, 0, 0)
+energy = openmc.stats.delta_function(14.1e6)  # At 14.1 MeV
+model.settings.source = openmc.IndependentSource(space=space, energy=energy)
+
+###############################################################################
+# Define tallies
+
+# Create time filter
+t_grid = np.insert(np.logspace(-6, 2, 100), 0, 0.0)
+time_filter = openmc.TimeFilter(t_grid)
+
+# Tally for total neutron density in time
+density_tally = openmc.Tally(name="Density")
+density_tally.filters = [time_filter]
+density_tally.scores = ["inverse-velocity"]
+
+# Add tallies to model
+model.tallies = openmc.Tallies([density_tally])
+
+
+# Run the model
+model.run(apply_tally_results=True)
+
+# Bin-averaged result
+density_mean = density_tally.mean.ravel() / np.diff(t_grid)
+
+# Plot particle density versus time
+fig, ax = plt.subplots()
+ax.stairs(density_mean, t_grid)
+ax.set_xscale("log")
+ax.set_yscale("log")
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("Total density")
+ax.grid()
+plt.show()

include/openmc/settings.h

@@ -1,3 +1,3 @@
 #ifndef OPENMC_SETTINGS_H
 #define OPENMC_SETTINGS_H
 
@@ -150,6 +150,7 @@
 
 extern int
   max_history_splits; //!< maximum number of particle splits for weight windows
+extern int max_secondaries; //!< maximum number of secondaries in the bank
 extern int64_t ssw_max_particles; //!< maximum number of particles to be
                                   //!< banked on surfaces per process
 extern int64_t ssw_max_files;     //!< maximum number of surface source files

openmc/settings.py

@@ -128,6 +128,10 @@ class Settings:
         Maximum number of times a particle can split during a history
 
         .. versionadded:: 0.13
+    max_secondaries : int
+        Maximum secondary bank size
+
+        .. versionadded:: 0.15.3
     max_tracks : int
         Maximum number of tracks written to a track file (per MPI process).
 
@@ -1137,6 +1141,16 @@ def max_history_splits(self, value: int):
         cv.check_greater_than('max particle splits', value, 0)
         self._max_history_splits = value
 
+    @property
+    def max_secondaries(self) -> int:
+        return self._max_secondaries
+
+    @max_secondaries.setter
+    def max_secondaries(self, value: int):
+        cv.check_type('maximum secondary bank size', value, Integral)
+        cv.check_greater_than('max secondary bank size', value, 0)
+        self._max_secondaries = value
+
     @property
     def max_tracks(self) -> int:
         return self._max_tracks
@@ -1673,6 +1687,11 @@ def _create_max_history_splits_subelement(self, root):
             elem = ET.SubElement(root, "max_history_splits")
             elem.text = str(self._max_history_splits)
 
+    def _create_max_secondaries_subelement(self, root):
+        if self._max_secondaries is not None:
+            elem = ET.SubElement(root, "max_secondaries")
+            elem.text = str(self._max_secondaries)
+
     def _create_max_tracks_subelement(self, root):
         if self._max_tracks is not None:
             elem = ET.SubElement(root, "max_tracks")
@@ -2073,6 +2092,11 @@ def _max_history_splits_from_xml_element(self, root):
         if text is not None:
             self.max_history_splits = int(text)
 
+    def _max_secondaries_from_xml_element(self, root):
+        text = get_text(root, 'max_secondaries')
+        if text is not None:
+            self.max_secondaries = int(text)
+
     def _max_tracks_from_xml_element(self, root):
         text = get_text(root, 'max_tracks')
         if text is not None:

src/eigenvalue.cpp

@@ -127,30 +127,8 @@ void synchronize_bank()
       "No fission sites banked on MPI rank " + std::to_string(mpi::rank));
   }
 
-  // Make sure all processors start at the same point for random sampling. Then
-  // skip ahead in the sequence using the starting index in the 'global'
-  // fission bank for each processor.
-
-  int64_t id = simulation::total_gen + overall_generation();
-  uint64_t seed = init_seed(id, STREAM_TRACKING);
-  advance_prn_seed(start, &seed);
-
-  // Determine how many fission sites we need to sample from the source bank
-  // and the probability for selecting a site.
-
-  int64_t sites_needed;
-  if (total < settings::n_particles) {
-    sites_needed = settings::n_particles % total;
-  } else {
-    sites_needed = settings::n_particles;
-  }
-  double p_sample = static_cast<double>(sites_needed) / total;
-
   simulation::time_bank_sample.start();
 
-  // ==========================================================================
-  // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
-
   // Allocate temporary source bank -- we don't really know how many fission
   // sites were created, so overallocate by a factor of 3
   int64_t index_temp = 0;
@@ -165,33 +143,38 @@ void synchronize_bank()
       temp_delayed_groups, temp_lifetimes, 3 * simulation::work_per_rank);
   }
 
-  for (int64_t i = 0; i < simulation::fission_bank.size(); i++) {
-    const auto& site = simulation::fission_bank[i];
+  // ==========================================================================
+  // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
 
-    // If there are less than n_particles particles banked, automatically add
-    // int(n_particles/total) sites to temp_sites. For example, if you need
-    // 1000 and 300 were banked, this would add 3 source sites per banked site
-    // and the remaining 100 would be randomly sampled.
-    if (total < settings::n_particles) {
-      for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
-        temp_sites[index_temp] = site;
-        if (settings::ifp_on) {
-          copy_ifp_data_from_fission_banks(
-            i, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]);
-        }
-        ++index_temp;
-      }
-    }
+  // We use Uniform Combing method to exactly get the targeted particle size
+  // [https://doi.org/10.1080/00295639.2022.2091906]
 
-    // Randomly sample sites needed
-    if (prn(&seed) < p_sample) {
-      temp_sites[index_temp] = site;
-      if (settings::ifp_on) {
-        copy_ifp_data_from_fission_banks(
-          i, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]);
-      }
-      ++index_temp;
+  // Make sure all processors use the same random number seed.
+  int64_t id = simulation::total_gen + overall_generation();
+  uint64_t seed = init_seed(id, STREAM_TRACKING);
+
+  // Comb specification
+  double teeth_distance = static_cast<double>(total) / settings::n_particles;
+  double teeth_offset = prn(&seed) * teeth_distance;
+
+  // First and last hitting tooth
+  int64_t end = start + simulation::fission_bank.size();
+  int64_t tooth_start = std::ceil((start - teeth_offset) / teeth_distance);
+  int64_t tooth_end = std::floor((end - teeth_offset) / teeth_distance) + 1;
+
+  // Locally comb particles in fission_bank
+  double tooth = tooth_start * teeth_distance + teeth_offset;
+  for (int64_t i = tooth_start; i < tooth_end; i++) {
+    int64_t idx = std::floor(tooth) - start;
+    temp_sites[index_temp] = simulation::fission_bank[idx];
+    if (settings::ifp_on) {
+      copy_ifp_data_from_fission_banks(
+        idx, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]);
     }
+    ++index_temp;
+
+    // Next tooth
+    tooth += teeth_distance;
   }
 
   // At this point, the sampling of source sites is done and now we need to
@@ -217,37 +200,6 @@ void synchronize_bank()
   finish = index_temp;
 #endif
 
-  // Now that the sampling is complete, we need to ensure that we have exactly
-  // n_particles source sites. The way this is done in a reproducible manner is
-  // to adjust only the source sites on the last processor.
-
-  if (mpi::rank == mpi::n_procs - 1) {
-    if (finish > settings::n_particles) {
-      // If we have extra sites sampled, we will simply discard the extra
-      // ones on the last processor
-      index_temp = settings::n_particles - start;
-
-    } else if (finish < settings::n_particles) {
-      // If we have too few sites, repeat sites from the very end of the
-      // fission bank
-      sites_needed = settings::n_particles - finish;
-      // TODO: sites_needed > simulation::fission_bank.size() or other test to
-      // make sure we don't need info from other proc
-      for (int i = 0; i < sites_needed; ++i) {
-        int i_bank = simulation::fission_bank.size() - sites_needed + i;
-        temp_sites[index_temp] = simulation::fission_bank[i_bank];
-        if (settings::ifp_on) {
-          copy_ifp_data_from_fission_banks(i_bank,
-            temp_delayed_groups[index_temp], temp_lifetimes[index_temp]);
-        }
-        ++index_temp;
-      }
-    }
-
-    // the last processor should not be sending sites to right
-    finish = simulation::work_index[mpi::rank + 1];
-  }
-
   simulation::time_bank_sample.stop();
   simulation::time_bank_sendrecv.start();
 

src/finalize.cpp

@@ -92,6 +92,7 @@ int openmc_finalize()
   settings::max_lost_particles = 10;
   settings::max_order = 0;
   settings::max_particles_in_flight = 100000;
+  settings::max_secondaries = 10000;
   settings::max_particle_events = 1'000'000;
   settings::max_history_splits = 10'000'000;
   settings::max_tracks = 1000;

src/physics.cpp

@@ -115,7 +115,7 @@ void sample_neutron_reaction(Particle& p)
 
       // Make sure particle population doesn't grow out of control for
       // subcritical multiplication problems.
-      if (p.secondary_bank().size() >= 10000) {
+      if (p.secondary_bank().size() >= settings::max_secondaries) {
         fatal_error(
           "The secondary particle bank appears to be growing without "
           "bound. You are likely running a subcritical multiplication problem "
@@ -210,13 +210,23 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     site.particle = ParticleType::neutron;
     site.time = p.time();
     site.wgt = 1. / weight;
-    site.parent_id = p.id();
-    site.progeny_id = p.n_progeny()++;
     site.surf_id = 0;
 
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, &site, p);
 
+    // Reject site if it exceeds time cutoff
+    if (site.delayed_group > 0) {
+      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > t_cutoff) {
+        continue;
+      }
+    }
+
+    // Set parent and progeny IDs
+    site.parent_id = p.id();
+    site.progeny_id = p.n_progeny()++;
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);
@@ -1069,6 +1079,10 @@ void sample_fission_neutron(
     // set the delayed group for the particle born from fission
     site->delayed_group = group;
 
+    // Sample time of emission based on decay constant of precursor
+    double decay_rate = rx.products_[site->delayed_group].decay_rate_;
+    site->time -= std::log(prn(p.current_seed())) / decay_rate;
+
   } else {
     // ====================================================================
     // PROMPT NEUTRON SAMPLED