Add delayed fission neutron emission

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

@@ -150,6 +150,7 @@ extern double source_rejection_fraction; //!< Minimum fraction of source sites
 
 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/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

src/physics_mg.cpp

@@ -139,8 +139,6 @@ void create_fission_sites(Particle& p)
     site.particle = ParticleType::neutron;
     site.time = p.time();
     site.wgt = 1. / weight;
-    site.parent_id = p.id();
-    site.progeny_id = p.n_progeny()++;
 
     // Sample the cosine of the angle, assuming fission neutrons are emitted
     // isotropically
@@ -165,6 +163,24 @@ void create_fission_sites(Particle& p)
     // of the code, 0 is prompt.
     site.delayed_group = dg + 1;
 
+    // If delayed product production, sample time of emission
+    if (dg != -1) {
+      auto& macro_xs = data::mg.macro_xs_[p.material()];
+      double decay_rate =
+        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
+      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > t_cutoff) {
+        continue;
+      }
+    }
+
+    // Set parent and progeny ID
+    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);

src/settings.cpp

@@ -114,6 +114,7 @@ int max_order {0};
 int n_log_bins {8000};
 int n_batches;
 int n_max_batches;
+int max_secondaries {10000};
 int max_history_splits {10'000'000};
 int max_tracks {1000};
 ResScatMethod res_scat_method {ResScatMethod::rvs};
@@ -1144,6 +1145,11 @@ void read_settings_xml(pugi::xml_node root)
     weight_windows_on = get_node_value_bool(root, "weight_windows_on");
   }
 
+  if (check_for_node(root, "max_secondaries")) {
+    settings::max_secondaries =
+      std::stoi(get_node_value(root, "max_secondaries"));
+  }
+
   if (check_for_node(root, "max_history_splits")) {
     settings::max_history_splits =
       std::stoi(get_node_value(root, "max_history_splits"));

src/tallies/tally_scoring.cpp

@@ -1624,8 +1624,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
       tally.estimator_ == TallyEstimator::COLLISION) {
     if (settings::survival_biasing) {
       // Determine weight that was absorbed
-      wgt_absorb = p.wgt_last() * p.neutron_xs(p.event_nuclide()).absorption /
-                   p.neutron_xs(p.event_nuclide()).total;
+      wgt_absorb = p.wgt_last() * p.macro_xs().absorption / p.macro_xs().total;
 
       // Then we either are alive and had a scatter (and so g changed),
       // or are dead and g did not change

tests/unit_tests/test_settings.py

@@ -67,6 +67,7 @@ def test_export_to_xml(run_in_tmpdir):
         )
     }
     s.max_particle_events = 100
+    s.max_secondaries = 1_000_000
     s.source_rejection_fraction = 0.01
 
     # Make sure exporting XML works
@@ -144,4 +145,5 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.random_ray['distance_active'] == 100.0
     assert s.random_ray['ray_source'].space.lower_left == [-1., -1., -1.]
     assert s.random_ray['ray_source'].space.upper_right == [1., 1., 1.]
+    assert s.max_secondaries == 1_000_000
     assert s.source_rejection_fraction == 0.01