From e5cf31ff9b9860f35def79965a6f4da7fc961092 Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Tue, 27 Feb 2024 23:15:40 -0800
Subject: [PATCH 01/13] add mg speed access, adjust fission emission time

---
 src/particle.cpp   |  9 +++++++++
 src/physics_mg.cpp | 16 ++++++++++++++++
 2 files changed, 25 insertions(+)

diff --git a/src/particle.cpp b/src/particle.cpp
index 549de5cff65..f69f2912536 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -44,6 +44,15 @@ namespace openmc {
 
 double Particle::speed() const
 {
+  // Multigroup speed?
+  if (!settings::run_CE) {
+    auto& macro_xs = data::mg.macro_xs_[material()];
+    int macro_t = mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(u());
+    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, g(), nullptr,
+                   nullptr, nullptr, macro_t, macro_a);
+  }
+
   // Determine mass in eV/c^2
   double mass;
   switch (this->type()) {
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 361cf5affcf..0be902ec813 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -130,6 +130,7 @@ void create_fission_sites(Particle& p)
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;
+    site.time = p.time();
 
     // Sample the cosine of the angle, assuming fission neutrons are emitted
     // isotropically
@@ -154,6 +155,21 @@ 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 time_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > time_cutoff) {
+        p.n_progeny()--;
+        continue;
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);

From ef0869ba9f000ce31c5b7bf571ca9a8776176daf Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Tue, 27 Feb 2024 23:16:32 -0800
Subject: [PATCH 02/13] clang-format

---
 src/physics_mg.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 0be902ec813..f950cb6c689 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -158,12 +158,13 @@ void create_fission_sites(Particle& p)
     // 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);
+      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 time_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      double time_cutoff =
+        settings::time_cutoff[static_cast<int>(site.particle)];
       if (site.time > time_cutoff) {
         p.n_progeny()--;
         continue;

From 508fa13f2b4d3213d1e60dabd3ec3298129d96de Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Wed, 28 Feb 2024 10:31:21 -0800
Subject: [PATCH 03/13] fix traveled distance at time futoff

---
 src/particle.cpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/particle.cpp b/src/particle.cpp
index f69f2912536..6740dd039d4 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -246,6 +246,9 @@ void Particle::event_advance()
     double push_back_distance = speed() * dt;
     this->move_distance(-push_back_distance);
     hit_time_boundary = true;
+
+    // Reduce the distance traveled for tallying
+    distance -= push_back_distance;
   }
 
   // Score track-length tallies

From ad07c3d9a6a7c28d32e6f41edfb52a5ed878f0da Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Wed, 28 Feb 2024 13:45:38 -0800
Subject: [PATCH 04/13] add delayed fission time

---
 src/physics.cpp    | 14 ++++++++++++++
 src/physics_mg.cpp |  2 +-
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/src/physics.cpp b/src/physics.cpp
index 69e74f2eafd..65ff5d4aa59 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -209,6 +209,20 @@ 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);
 
+    // If delayed product production, sample time of emission
+    if (site.delayed_group > 0) {
+      double decay_rate = rx.products_[site.delayed_group].decay_rate_;
+      site.time -= std::log(prn(p.current_seed())) / decay_rate;
+
+      // Reject site if it exceeds time cutoff
+      double time_cutoff =
+        settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > time_cutoff) {
+        p.n_progeny()--;
+        continue;
+      }
+    }
+
     // Store fission site in bank
     if (use_fission_bank) {
       int64_t idx = simulation::fission_bank.thread_safe_append(site);
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index f950cb6c689..7eaec98236b 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -127,10 +127,10 @@ void create_fission_sites(Particle& p)
     SourceSite site;
     site.r = p.r();
     site.particle = ParticleType::neutron;
+    site.time = p.time();
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;
-    site.time = p.time();
 
     // Sample the cosine of the angle, assuming fission neutrons are emitted
     // isotropically

From c72f753b16009c28b40012dc3433c5da0c2d4d97 Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Wed, 28 Feb 2024 16:59:48 -0800
Subject: [PATCH 05/13] update docs: sampling delay time

---
 docs/source/methods/neutron_physics.rst | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst
index 70ace4a3532..1117de83779 100644
--- a/docs/source/methods/neutron_physics.rst
+++ b/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 calculated 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

From 08c883e6b39ff00c0b79092b973fbe07c9d6a688 Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Wed, 28 Feb 2024 19:57:39 -0800
Subject: [PATCH 06/13] add settings:max_secondaries

---
 include/openmc/settings.h |  1 +
 openmc/settings.py        | 27 +++++++++++++++++++++++++++
 src/finalize.cpp          |  1 +
 src/physics.cpp           |  2 +-
 src/settings.cpp          |  5 +++++
 5 files changed, 35 insertions(+), 1 deletion(-)

diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 69a8d7d13be..2921ca23c19 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -123,6 +123,7 @@ extern std::unordered_set<int>
 extern int max_splits; //!< maximum number of particle splits for weight windows
 extern int64_t max_surface_particles; //!< maximum number of particles to be
                                       //!< banked on surfaces per process
+extern int max_secondaries; //!< maximum number of secondaries in the bank
 extern TemperatureMethod
   temperature_method; //!< method for choosing temperatures
 extern double
diff --git a/openmc/settings.py b/openmc/settings.py
index a6d273901e0..4cf8d243ff3 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -114,6 +114,10 @@ class Settings:
     max_splits : int
         Maximum number of times a particle can split during a history
 
+        .. versionadded:: 0.15
+    max_secondaries : int
+        Maximum secondary bank size
+
         .. versionadded:: 0.13
     max_tracks : int
         Maximum number of tracks written to a track file (per MPI process).
@@ -341,6 +345,7 @@ def __init__(self, **kwargs):
         self._weight_windows_file = None
         self._weight_window_checkpoints = {}
         self._max_splits = None
+        self._max_secondaries = None
         self._max_tracks = None
 
         for key, value in kwargs.items():
@@ -986,6 +991,16 @@ def max_splits(self, value: int):
         cv.check_greater_than('max particle splits', value, 0)
         self._max_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
@@ -1417,6 +1432,11 @@ def _create_max_splits_subelement(self, root):
             elem = ET.SubElement(root, "max_splits")
             elem.text = str(self._max_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")
@@ -1763,6 +1783,11 @@ def _max_splits_from_xml_element(self, root):
         if text is not None:
             self.max_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:
@@ -1828,6 +1853,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_weight_windows_file_element(element)
         self._create_weight_window_checkpoints_subelement(element)
         self._create_max_splits_subelement(element)
+        self._create_max_secondaries_subelement(element)
         self._create_max_tracks_subelement(element)
 
         # Clean the indentation in the file to be user-readable
@@ -1930,6 +1956,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._weight_window_generators_from_xml_element(elem, meshes)
         settings._weight_window_checkpoints_from_xml_element(elem)
         settings._max_splits_from_xml_element(elem)
+        settings._max_secondaries_from_xml_element(elem)
         settings._max_tracks_from_xml_element(elem)
 
         # TODO: Get volume calculations
diff --git a/src/finalize.cpp b/src/finalize.cpp
index a41c2f783a5..c6a5cfcb21c 100644
--- a/src/finalize.cpp
+++ b/src/finalize.cpp
@@ -91,6 +91,7 @@ int openmc_finalize()
   settings::max_lost_particles = 10;
   settings::max_order = 0;
   settings::max_particles_in_flight = 100000;
+  settings::max_secondaries = 100000;
   settings::max_splits = 1000;
   settings::max_tracks = 1000;
   settings::max_write_lost_particles = -1;
diff --git a/src/physics.cpp b/src/physics.cpp
index 65ff5d4aa59..ebea5a22515 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -113,7 +113,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 "
diff --git a/src/settings.cpp b/src/settings.cpp
index a5256c9c26d..77d973112d7 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -105,6 +105,7 @@ int max_order {0};
 int n_log_bins {8000};
 int n_batches;
 int n_max_batches;
+int max_secondaries {10000};
 int max_splits {1000};
 int max_tracks {1000};
 ResScatMethod res_scat_method {ResScatMethod::rvs};
@@ -930,6 +931,10 @@ 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_splits")) {
     settings::max_splits = std::stoi(get_node_value(root, "max_splits"));
   }

From 14f96c962e814217089ad4cfcf770b6de546c6cf Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Thu, 29 Feb 2024 10:47:00 -0800
Subject: [PATCH 07/13] update default max_secondaries at in finalize.cpp

---
 src/finalize.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/finalize.cpp b/src/finalize.cpp
index c6a5cfcb21c..70722add4ba 100644
--- a/src/finalize.cpp
+++ b/src/finalize.cpp
@@ -91,7 +91,7 @@ int openmc_finalize()
   settings::max_lost_particles = 10;
   settings::max_order = 0;
   settings::max_particles_in_flight = 100000;
-  settings::max_secondaries = 100000;
+  settings::max_secondaries = 10000;
   settings::max_splits = 1000;
   settings::max_tracks = 1000;
   settings::max_write_lost_particles = -1;

From e388596fa44f09d820dd85eea9ec4fb66c77dbdc Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansi@oregonstate.edu>
Date: Mon, 11 Mar 2024 09:09:10 -0700
Subject: [PATCH 08/13] Update docs/source/methods/neutron_physics.rst

Co-authored-by: Gavin Ridley <gavin.keith.ridley@gmail.com>
---
 docs/source/methods/neutron_physics.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst
index 1117de83779..9515055692f 100644
--- a/docs/source/methods/neutron_physics.rst
+++ b/docs/source/methods/neutron_physics.rst
@@ -310,7 +310,7 @@ interpolation law. The number of prompt neutrons released per fission event
 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
+:math:`\nu_d`. Once these have been determined, we can calculate the delayed
 neutron fraction
 
 .. math::

From f729e6a2ecf6409f2497325bffb801298f12feb2 Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Wed, 13 Mar 2024 13:19:04 -0700
Subject: [PATCH 09/13] add pulsed pincell example

---
 examples/pincell_pulsed/build_xml.py   | 96 ++++++++++++++++++++++++++
 examples/pincell_pulsed/plot_result.py | 23 ++++++
 2 files changed, 119 insertions(+)
 create mode 100644 examples/pincell_pulsed/build_xml.py
 create mode 100644 examples/pincell_pulsed/plot_result.py

diff --git a/examples/pincell_pulsed/build_xml.py b/examples/pincell_pulsed/build_xml.py
new file mode 100644
index 00000000000..d8f5b943ee8
--- /dev/null
+++ b/examples/pincell_pulsed/build_xml.py
@@ -0,0 +1,96 @@
+from math import log10
+
+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")
+
+# Collect the materials together and export to XML
+materials = openmc.Materials([uo2, helium, zircaloy, borated_water])
+materials.export_to_xml()
+
+###############################################################################
+# 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 geometry and export to XML
+geometry = openmc.Geometry([fuel, gap, clad, water])
+geometry.export_to_xml()
+
+###############################################################################
+# Define problem settings
+
+# Set the mode
+settings = openmc.Settings()
+settings.run_mode = "fixed source"
+
+# Indicate how many batches and particles to run
+settings.batches = 10
+settings.particles = 10000
+
+# Set time cutoff
+# (because we only care about solutions in t < 100 seconds, see tally below)
+settings.cutoff = {"time_neutron": 100}
+
+# Create the neutron pulse source
+# (by default, isotropic direction, at t = 0 s)
+space = openmc.stats.Point()  # At the origin (0, 0, 0)
+energy = openmc.stats.Discrete([1.41e7], [1.0])  # At 14.1 MeV
+settings.source = openmc.IndependentSource(space=space, energy=energy)
+
+# Export to XML
+settings.export_to_xml()
+
+###############################################################################
+# 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"]
+
+# Instantiate a Tallies collection and export to XML
+tallies = openmc.Tallies([density_tally])
+tallies.export_to_xml()
diff --git a/examples/pincell_pulsed/plot_result.py b/examples/pincell_pulsed/plot_result.py
new file mode 100644
index 00000000000..6d4aae811f4
--- /dev/null
+++ b/examples/pincell_pulsed/plot_result.py
@@ -0,0 +1,23 @@
+import matplotlib.pyplot as plt
+import openmc
+
+
+# Get results from statepoint
+with openmc.StatePoint("statepoint.10.h5") as sp:
+    tally = sp.get_tally(name="Density")
+
+# Get the time grid
+t = tally.filters[0].bins
+dt = t[:, 1] - t[:, 0]
+t_mid = 0.5 * (t[:, 0] + t[:, 1])
+
+# Bin-averaged result
+density_mean = tally.get_values(value="mean").ravel() / dt
+
+# Plot flux spectrum
+fig, ax = plt.subplots()
+ax.loglog(t_mid, density_mean, "ok", fillstyle="none")
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("Total density")
+ax.grid()
+plt.show()

From 6ca1ae09547dd78684b31a1bfe155d6eed141b2b Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansyah.ilham@gmail.com>
Date: Wed, 3 Apr 2024 10:08:54 -0700
Subject: [PATCH 10/13] debug mg_mode + collision_estimator + survival_biasing

---
 src/tallies/tally_scoring.cpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index e798161ec2f..fc3236b8aee 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -1558,8 +1558,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

From 718356f32bf9b83d243cd3ae857df8de2be8e406 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 24 Apr 2024 11:56:40 -0500
Subject: [PATCH 11/13] Small fixes/updates

---
 examples/pincell_pulsed/build_xml.py   | 32 +++++++++-----------------
 examples/pincell_pulsed/plot_result.py |  9 ++++----
 openmc/settings.py                     |  4 ++--
 src/particle.cpp                       |  8 +++----
 src/physics.cpp                        | 19 ++++++---------
 src/physics_mg.cpp                     | 17 +++++---------
 6 files changed, 35 insertions(+), 54 deletions(-)

diff --git a/examples/pincell_pulsed/build_xml.py b/examples/pincell_pulsed/build_xml.py
index d8f5b943ee8..7cd1617dfeb 100644
--- a/examples/pincell_pulsed/build_xml.py
+++ b/examples/pincell_pulsed/build_xml.py
@@ -1,5 +1,3 @@
-from math import log10
-
 import numpy as np
 import openmc
 
@@ -29,10 +27,6 @@
 borated_water.add_element("O", 2.4e-2)
 borated_water.add_s_alpha_beta("c_H_in_H2O")
 
-# Collect the materials together and export to XML
-materials = openmc.Materials([uo2, helium, zircaloy, borated_water])
-materials.export_to_xml()
-
 ###############################################################################
 # Define problem geometry
 
@@ -51,33 +45,29 @@
 clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
 water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
 
-# Create a geometry and export to XML
-geometry = openmc.Geometry([fuel, gap, clad, water])
-geometry.export_to_xml()
+# Create a model and assign geometry
+model = openmc.Model()
+model.geometry = openmc.Geometry([fuel, gap, clad, water])
 
 ###############################################################################
 # Define problem settings
 
 # Set the mode
-settings = openmc.Settings()
-settings.run_mode = "fixed source"
+model.settings.run_mode = "fixed source"
 
 # Indicate how many batches and particles to run
-settings.batches = 10
-settings.particles = 10000
+model.settings.batches = 10
+model.settings.particles = 10000
 
 # Set time cutoff
 # (because we only care about solutions in t < 100 seconds, see tally below)
-settings.cutoff = {"time_neutron": 100}
+model.settings.cutoff = {"time_neutron": 100}
 
 # Create the neutron pulse source
 # (by default, isotropic direction, at t = 0 s)
 space = openmc.stats.Point()  # At the origin (0, 0, 0)
 energy = openmc.stats.Discrete([1.41e7], [1.0])  # At 14.1 MeV
-settings.source = openmc.IndependentSource(space=space, energy=energy)
-
-# Export to XML
-settings.export_to_xml()
+model.settings.source = openmc.IndependentSource(space=space, energy=energy)
 
 ###############################################################################
 # Define tallies
@@ -91,6 +81,6 @@
 density_tally.filters = [time_filter]
 density_tally.scores = ["inverse-velocity"]
 
-# Instantiate a Tallies collection and export to XML
-tallies = openmc.Tallies([density_tally])
-tallies.export_to_xml()
+# Add tallies to model
+model.tallies = openmc.Tallies([density_tally])
+model.export_to_model_xml()
diff --git a/examples/pincell_pulsed/plot_result.py b/examples/pincell_pulsed/plot_result.py
index 6d4aae811f4..29dae470c5f 100644
--- a/examples/pincell_pulsed/plot_result.py
+++ b/examples/pincell_pulsed/plot_result.py
@@ -1,4 +1,5 @@
 import matplotlib.pyplot as plt
+import numpy as np
 import openmc
 
 
@@ -7,12 +8,12 @@
     tally = sp.get_tally(name="Density")
 
 # Get the time grid
-t = tally.filters[0].bins
-dt = t[:, 1] - t[:, 0]
-t_mid = 0.5 * (t[:, 0] + t[:, 1])
+t = tally.filters[0].values
+dt = np.diff(t)
+t_mid = 0.5 * (t[:-1] + t[1:])
 
 # Bin-averaged result
-density_mean = tally.get_values(value="mean").ravel() / dt
+density_mean = tally.mean.ravel() / dt
 
 # Plot flux spectrum
 fig, ax = plt.subplots()
diff --git a/openmc/settings.py b/openmc/settings.py
index 4cf8d243ff3..35ccd285920 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -114,11 +114,11 @@ class Settings:
     max_splits : int
         Maximum number of times a particle can split during a history
 
-        .. versionadded:: 0.15
+        .. versionadded:: 0.13
     max_secondaries : int
         Maximum secondary bank size
 
-        .. versionadded:: 0.13
+        .. versionadded:: 0.14.1
     max_tracks : int
         Maximum number of tracks written to a track file (per MPI process).
 
diff --git a/src/particle.cpp b/src/particle.cpp
index 6740dd039d4..fecb411eff1 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -46,10 +46,10 @@ double Particle::speed() const
 {
   // Multigroup speed?
   if (!settings::run_CE) {
-    auto& macro_xs = data::mg.macro_xs_[material()];
-    int macro_t = mg_xs_cache().t;
-    int macro_a = macro_xs.get_angle_index(u());
-    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, g(), nullptr,
+    auto& macro_xs = data::mg.macro_xs_[this->material()];
+    int macro_t = this->mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(this->u());
+    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
                    nullptr, nullptr, macro_t, macro_a);
   }
 
diff --git a/src/physics.cpp b/src/physics.cpp
index ebea5a22515..683293bfb49 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -202,23 +202,19 @@ 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);
 
-    // If delayed product production, sample time of emission
+    // If neutron is delayed, sample time of emission
     if (site.delayed_group > 0) {
       double decay_rate = rx.products_[site.delayed_group].decay_rate_;
       site.time -= std::log(prn(p.current_seed())) / decay_rate;
 
       // Reject site if it exceeds time cutoff
-      double time_cutoff =
-        settings::time_cutoff[static_cast<int>(site.particle)];
-      if (site.time > time_cutoff) {
-        p.n_progeny()--;
+      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > t_cutoff) {
         continue;
       }
     }
@@ -232,11 +228,6 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
           "in this generation will not be banked. Results may be "
           "non-deterministic.");
 
-        // Decrement number of particle progeny as storage was unsuccessful.
-        // This step is needed so that the sum of all progeny is equal to the
-        // size of the shared fission bank.
-        p.n_progeny()--;
-
         // Break out of loop as no more sites can be added to fission bank
         break;
       }
@@ -244,6 +235,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
       p.secondary_bank().push_back(site);
     }
 
+    // Set parent and progeny ID
+    site.parent_id = p.id();
+    site.progeny_id = p.n_progeny()++;
+
     // Set the delayed group on the particle as well
     p.delayed_group() = site.delayed_group;
 
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 7eaec98236b..cbe954cc9a4 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -129,8 +129,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
@@ -163,10 +161,8 @@ void create_fission_sites(Particle& p)
       site.time -= std::log(prn(p.current_seed())) / decay_rate;
 
       // Reject site if it exceeds time cutoff
-      double time_cutoff =
-        settings::time_cutoff[static_cast<int>(site.particle)];
-      if (site.time > time_cutoff) {
-        p.n_progeny()--;
+      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      if (site.time > t_cutoff) {
         continue;
       }
     }
@@ -180,11 +176,6 @@ void create_fission_sites(Particle& p)
           "in this generation will not be banked. Results may be "
           "non-deterministic.");
 
-        // Decrement number of particle progeny as storage was unsuccessful.
-        // This step is needed so that the sum of all progeny is equal to the
-        // size of the shared fission bank.
-        p.n_progeny()--;
-
         // Break out of loop as no more sites can be added to fission bank
         break;
       }
@@ -192,6 +183,10 @@ void create_fission_sites(Particle& p)
       p.secondary_bank().push_back(site);
     }
 
+    // Set parent and progeny ID
+    site.parent_id = p.id();
+    site.progeny_id = p.n_progeny()++;
+
     // Set the delayed group on the particle as well
     p.delayed_group() = dg + 1;
 

From 6b2a92ad3f37073ec60d82afa8e48ee4e31945c2 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 24 Apr 2024 12:04:40 -0500
Subject: [PATCH 12/13] Add max_secondaries on settings unit test

---
 tests/unit_tests/test_settings.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py
index 650bfd18680..79bbf49997d 100644
--- a/tests/unit_tests/test_settings.py
+++ b/tests/unit_tests/test_settings.py
@@ -65,8 +65,8 @@ def test_export_to_xml(run_in_tmpdir):
             space=openmc.stats.Box((-1., -1., -1.), (1., 1., 1.))
         )
     }
-
     s.max_particle_events = 100
+    s.max_secondaries = 1_000_000
 
     # Make sure exporting XML works
     s.export_to_xml()
@@ -142,3 +142,4 @@ 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

From 74c82c2a46caf4dac155fc196737d9ca606393e0 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 24 Apr 2024 12:07:10 -0500
Subject: [PATCH 13/13] Fix C++ formatting on settings.cpp

---
 src/settings.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/settings.cpp b/src/settings.cpp
index 90034ddf79d..b1f5da39a34 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -981,7 +981,8 @@ void read_settings_xml(pugi::xml_node root)
   }
 
   if (check_for_node(root, "max_secondaries")) {
-    settings::max_secondaries = std::stoi(get_node_value(root, "max_secondaries"));
+    settings::max_secondaries =
+      std::stoi(get_node_value(root, "max_secondaries"));
   }
 
   if (check_for_node(root, "max_splits")) {