From 868ef8899de67783a7b179c40977c4eeaaa38cd0 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Thu, 7 Aug 2025 12:45:47 +0300
Subject: [PATCH 01/10] fix for time_filter & mesh_filter & tracklength filter

---
 include/openmc/tallies/tally.h |   4 ++
 src/particle.cpp               |  56 +++++++++-------
 src/tallies/tally.cpp          |  58 +++++++++++++----
 tests/unit_tests/test_mesh.py  | 115 +++++++++++++++++++++++++++++++++
 4 files changed, 200 insertions(+), 33 deletions(-)

diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 8c088c46089..29337ae51f7 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -208,6 +208,10 @@ extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
 extern vector<int> active_pulse_height_tallies;
 extern vector<int> pulse_height_cells;
+extern vector<double> time_grid;
+
+double distance_to_time_boundary(double time, double speed);
+
 } // namespace model
 
 namespace simulation {
diff --git a/src/particle.cpp b/src/particle.cpp
index a0f1d66bd5f..115c36759d4 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -243,39 +243,51 @@ void Particle::event_advance()
     (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
 
   // Select smaller of the three distances
-  double distance =
+  double total_distance =
     std::min({boundary().distance(), collision_distance(), distance_cutoff});
 
-  // Advance particle in space and time
-  // Short-term solution until the surface source is revised and we can use
-  // this->move_distance(distance)
-  for (int j = 0; j < n_coord(); ++j) {
-    coord(j).r() += distance * coord(j).u();
-  }
-  double dt = distance / speed;
-  this->time() += dt;
-  this->lifetime() += dt;
+  double distance_traveled = 0.0;
+  while (distance_traveled < total_distance) {
+
+    double distance = std::min(model::distance_to_time_boundary(time(), speed),
+      total_distance - distance_traveled);
+    double dt = distance / speed;
+
+    // Save particle last state for tracklength tallies
+    this->time_last() = time();
+    this->r_last() = r();
+
+    // Advance particle in space and time
+    // Short-term solution until the surface source is revised and we can use
+    // this->move_distance(distance)
+    for (int j = 0; j < n_coord(); ++j) {
+      coord(j).r() += distance * coord(j).u();
+    }
+    this->time() += dt;
+    this->lifetime() += dt;
+    distance_traveled += distance;
 
-  // Score track-length tallies
-  if (!model::active_tracklength_tallies.empty()) {
-    score_tracklength_tally(*this, distance);
+    // Score track-length tallies
+    if (!model::active_tracklength_tallies.empty()) {
+      score_tracklength_tally(*this, distance);
+    }
+
+    // Score flux derivative accumulators for differential tallies.
+    if (!model::active_tallies.empty()) {
+      score_track_derivative(*this, distance);
+    }
+>>>>>>> 31a712bb0 (fix for time_filter & mesh_filter & tracklength filter)
   }
 
   // Score track-length estimate of k-eff
   if (settings::run_mode == RunMode::EIGENVALUE &&
       type() == ParticleType::neutron) {
-    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
+    keff_tally_tracklength() += wgt() * total_distance * macro_xs().nu_fission;
   }
 
-  // Score flux derivative accumulators for differential tallies.
-  if (!model::active_tallies.empty()) {
-    score_track_derivative(*this, distance);
-  }
-
-  // Set particle weight to zero if it hit the time boundary
-  if (distance == distance_cutoff) {
+  // Set particle weight to zero if it hit the time cutoff boundary
+  if (total_distance == distance_cutoff)
     wgt() = 0.0;
-  }
 }
 
 void Particle::event_cross_surface()
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index ba183860269..11f7c13c8f8 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -27,10 +27,12 @@
 #include "openmc/tallies/filter_legendre.h"
 #include "openmc/tallies/filter_mesh.h"
 #include "openmc/tallies/filter_meshborn.h"
+#include "openmc/tallies/filter_meshmaterial.h"
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_particle.h"
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_surface.h"
+#include "openmc/tallies/filter_time.h"
 #include "openmc/xml_interface.h"
 
 #include "xtensor/xadapt.hpp"
@@ -61,6 +63,30 @@ vector<int> active_meshsurf_tallies;
 vector<int> active_surface_tallies;
 vector<int> active_pulse_height_tallies;
 vector<int> pulse_height_cells;
+vector<double> time_grid;
+
+double distance_to_time_boundary(double time, double speed)
+{
+  if (time_grid.size() == 0) {
+    return INFTY;
+  } else if (time >= time_grid[time_grid.size() - 1]) {
+    return INFTY;
+  } else {
+    return (*std::upper_bound(time_grid.begin(), time_grid.end(), time) -
+             time) *
+           speed;
+  }
+}
+
+void add_to_time_grid(vector<double> grid)
+{
+  auto temp = time_grid;
+  time_grid.resize(time_grid.size() + grid.size());
+  std::merge(temp.begin(), temp.end(), grid.begin(), grid.end(),
+    std::back_inserter(time_grid));
+  auto last_unique = std::unique(time_grid.begin(), time_grid.end());
+  time_grid.erase(last_unique, time_grid.end());
+}
 } // namespace model
 
 namespace simulation {
@@ -243,8 +269,8 @@ Tally::Tally(pugi::xml_node node)
     for (int score : scores_) {
       switch (score) {
       case SCORE_PULSE_HEIGHT:
-        fatal_error(
-          "For pulse-height tallies, photon transport needs to be activated.");
+        fatal_error("For pulse-height tallies, photon transport needs to be "
+                    "activated.");
         break;
       }
     }
@@ -318,7 +344,8 @@ Tally::Tally(pugi::xml_node node)
         if (has_energyout && i_nuc == -1) {
           fatal_error(fmt::format(
             "Error on tally {}: Cannot use a "
-            "'nuclide_density' or 'temperature' derivative on a tally with an "
+            "'nuclide_density' or 'temperature' derivative on a tally with "
+            "an "
             "outgoing energy filter and 'total' nuclide rate. Instead, tally "
             "each nuclide in the material individually.",
             id_));
@@ -493,9 +520,9 @@ void Tally::add_filter(Filter* filter)
 
 void Tally::set_strides()
 {
-  // Set the strides.  Filters are traversed in reverse so that the last filter
-  // has the shortest stride in memory and the first filter has the longest
-  // stride.
+  // Set the strides.  Filters are traversed in reverse so that the last
+  // filter has the shortest stride in memory and the first filter has the
+  // longest stride.
   auto n = filters_.size();
   strides_.resize(n, 0);
   int stride = 1;
@@ -551,7 +578,8 @@ void Tally::set_scores(const vector<std::string>& scores)
 
   // Iterate over the given scores.
   for (auto score_str : scores) {
-    // Make sure a delayed group filter wasn't used with an incompatible score.
+    // Make sure a delayed group filter wasn't used with an incompatible
+    // score.
     if (delayedgroup_filter_ != C_NONE) {
       if (score_str != "delayed-nu-fission" && score_str != "decay-rate")
         fatal_error("Cannot tally " + score_str + "with a delayedgroup filter");
@@ -984,8 +1012,8 @@ void reduce_tally_results()
     }
   }
 
-  // Note that global tallies are *always* reduced even when no_reduce option is
-  // on.
+  // Note that global tallies are *always* reduced even when no_reduce option
+  // is on.
 
   // Get view of global tally values
   auto& gt = simulation::global_tallies;
@@ -1073,6 +1101,7 @@ void setup_active_tallies()
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
+  model::time_grid.clear();
 
   for (auto i = 0; i < model::tallies.size(); ++i) {
     const auto& tally {*model::tallies[i]};
@@ -1088,6 +1117,12 @@ void setup_active_tallies()
           break;
         case TallyEstimator::TRACKLENGTH:
           model::active_tracklength_tallies.push_back(i);
+          if (auto time_filter = tally.get_filter<TimeFilter>()) {
+            if ((tally.get_filter<MeshFilter>() != nullptr) ||
+                (tally.get_filter<MeshMaterialFilter>() != nullptr)) {
+              model::add_to_time_grid(time_filter->bins());
+            }
+          }
           break;
         case TallyEstimator::COLLISION:
           model::active_collision_tallies.push_back(i);
@@ -1127,6 +1162,7 @@ void free_memory_tally()
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
+  model::time_grid.clear();
 
   model::tally_map.clear();
 }
@@ -1465,8 +1501,8 @@ extern "C" int openmc_tally_get_n_realizations(int32_t index, int32_t* n)
   return 0;
 }
 
-//! \brief Returns a pointer to a tally results array along with its shape. This
-//! allows a user to obtain in-memory tally results from Python directly.
+//! \brief Returns a pointer to a tally results array along with its shape.
+//! This allows a user to obtain in-memory tally results from Python directly.
 extern "C" int openmc_tally_results(
   int32_t index, double** results, size_t* shape)
 {
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 67ca4028e06..e84c3d5caec 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -7,6 +7,7 @@
 import openmc
 import openmc.lib
 from openmc.utility_funcs import change_directory
+from uncertainties.unumpy import uarray, nominal_values, std_devs
 
 
 @pytest.mark.parametrize("val_left,val_right", [(0, 0), (-1., -1.), (2.0, 2)])
@@ -651,3 +652,117 @@ def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
 
     # Run the model; this should not cause an infinite loop
     model.run()
+
+
+def test_filter_time_mesh(run_in_tmpdir):
+    """
+    Test that TimeFilter+MeshFilter with collision estimator agree with tracklength estimator.
+    """
+
+    # ===========================================================================
+    # Set Material
+    # ===========================================================================
+
+    mat = openmc.Material()
+    mat.add_nuclide('Fe56', 1.0)
+    mat.set_density('g/cm3', 7.8)
+
+    # ===========================================================================
+    # Set geometry
+    # ===========================================================================
+
+    # Instantiate ZCylinder surfaces
+    surf_Z1 = openmc.XPlane(surface_id=1, x0=-1e10, boundary_type="reflective")
+    surf_Z2 = openmc.XPlane(surface_id=2, x0=1e10, boundary_type="reflective")
+
+    # Instantiate Cells
+    cell_F = openmc.Cell(cell_id=1, name="F")
+
+    # Use surface half-spaces to define regions
+    cell_F.region = +surf_Z1 & -surf_Z2
+
+    # Register Materials with Cells
+    cell_F.fill = mat
+
+    # Instantiate Universes
+    root = openmc.Universe(universe_id=0, name="root universe", cells=[cell_F])
+
+    # Instantiate a Geometry, register the root Universe, and export to XML
+    geometry = openmc.Geometry(root)
+
+    # ===========================================================================
+    # Settings
+    # ===========================================================================
+
+    # Instantiate a Settings object, set all runtime parameters, and export to XML
+    settings_file = openmc.Settings()
+    settings_file.run_mode = "fixed source"
+    settings_file.particles = 10000
+    settings_file.batches = 20
+    settings_file.output = {"tallies": False}
+    settings_file.cutoff = {"time_neutron": 1E-7}
+
+    # Create an initial uniform spatial source distribution over fissionable zones
+    delta_dist = openmc.stats.Point()
+    isotropic = openmc.stats.Isotropic()
+    settings_file.source = openmc.IndependentSource(space=delta_dist, angle=isotropic)
+
+    # ===========================================================================
+    # Set tallies
+    # ===========================================================================
+
+    # Create a mesh filter that can be used in a tally
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (21, 1, 1)
+    mesh.lower_left = (-20.5, -1e10, -1e10)
+    mesh.upper_right = (20.5, 1e10, 1e10)
+    time_grid = np.linspace(0.0, 1E-7, 21)
+
+    mesh_filter = openmc.MeshFilter(mesh)
+    time_filter = openmc.TimeFilter(time_grid)
+
+    # Now use the mesh filter in a tally and indicate what scores are desired
+    tally1 = openmc.Tally(name="collision")
+    tally1.estimator = "collision"
+    tally1.filters = [time_filter, mesh_filter]
+    tally1.scores = ["flux"]
+
+    tally2 = openmc.Tally(name="tracklength")
+    tally2.estimator = "tracklength"
+    tally2.filters = [time_filter, mesh_filter]
+    tally2.scores = ["flux"]
+
+    # Instantiate a tallies collection
+    tallies = openmc.Tallies([tally1, tally2])
+    
+    # ===========================================================================
+    # Set the model
+    # ===========================================================================
+
+    model = openmc.Model()
+    model.geometry = geometry
+    model.settings = settings_file
+    model.tallies = tallies
+
+    # ===========================================================================
+    # Run and post-process
+    # ===========================================================================
+
+    sp_filename = model.run()
+
+    # Get radial flux distribution
+    with openmc.StatePoint(sp_filename) as sp:
+        flux_collision = sp.tallies[tally1.id].mean.ravel()
+        flux_collision_unc = sp.tallies[tally1.id].std_dev.ravel()
+        flux_tracklength = sp.tallies[tally2.id].mean.ravel()
+        flux_tracklength_unc = sp.tallies[tally2.id].std_dev.ravel()
+
+    # Construct arrays with uncertainties
+    collision = uarray(flux_collision, flux_collision_unc)
+    tracklength = uarray(flux_tracklength, flux_tracklength_unc)
+    delta = collision - tracklength
+
+    # Check that difference is within uncertainty
+    diff = nominal_values(delta)
+    std_dev = std_devs(delta)
+    assert np.all(diff < 4*std_dev)

From c8f2d25f4ab97c46e0b7c83e20f9ce0a8b0264a8 Mon Sep 17 00:00:00 2001
From: GuySten <62616591+GuySten@users.noreply.github.com>
Date: Mon, 18 Aug 2025 20:58:44 +0200
Subject: [PATCH 02/10] Update particle.cpp

---
 src/particle.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/particle.cpp b/src/particle.cpp
index 115c36759d4..79f5cbaec0c 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -276,7 +276,6 @@ void Particle::event_advance()
     if (!model::active_tallies.empty()) {
       score_track_derivative(*this, distance);
     }
->>>>>>> 31a712bb0 (fix for time_filter & mesh_filter & tracklength filter)
   }
 
   // Score track-length estimate of k-eff

From ae31b241bc606f299f2da0459382b9a9afe8f479 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 19 Aug 2025 04:34:11 +0200
Subject: [PATCH 03/10] change approach of code to leave transport logic as is
 and virtually move particle between time boundaries in the tally scoring
 section

---
 include/openmc/tallies/tally.h         |  1 +
 include/openmc/tallies/tally_scoring.h | 12 +++++
 src/particle.cpp                       | 55 ++++++++------------
 src/tallies/tally.cpp                  | 17 ++++---
 src/tallies/tally_scoring.cpp          | 61 +++++++++++++++++++++--
 tests/unit_tests/test_mesh.py          | 69 +++++++++++++-------------
 6 files changed, 138 insertions(+), 77 deletions(-)

diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 29337ae51f7..96687efcc3a 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -203,6 +203,7 @@ extern vector<unique_ptr<Tally>> tallies;
 extern vector<int> active_tallies;
 extern vector<int> active_analog_tallies;
 extern vector<int> active_tracklength_tallies;
+extern vector<int> active_timed_tracklength_tallies;
 extern vector<int> active_collision_tallies;
 extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
index 28f1f16223d..454a7db778f 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -91,6 +91,18 @@ void score_analog_tally_mg(Particle& p);
 //! \param distance The distance in [cm] traveled by the particle
 void score_tracklength_tally(Particle& p, double distance);
 
+//! Score tallies using a tracklength estimate of the flux.
+//
+//! This is triggered at every event (surface crossing, lattice crossing, or
+//! collision) and thus cannot be done for tallies that require post-collision
+//! information.
+//
+//! \param p The particle being tracked
+//! \param flux The tracklength flux estimate in [cm].
+//! \param tallies A vector of the indices of the tallies to score to
+void score_tracklength_tally(
+  Particle& p, double flux, const vector<int>& tallies);
+
 //! Score surface or mesh-surface tallies for particle currents.
 //
 //! \param p The particle being tracked
diff --git a/src/particle.cpp b/src/particle.cpp
index 79f5cbaec0c..a0f1d66bd5f 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -243,50 +243,39 @@ void Particle::event_advance()
     (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
 
   // Select smaller of the three distances
-  double total_distance =
+  double distance =
     std::min({boundary().distance(), collision_distance(), distance_cutoff});
 
-  double distance_traveled = 0.0;
-  while (distance_traveled < total_distance) {
-
-    double distance = std::min(model::distance_to_time_boundary(time(), speed),
-      total_distance - distance_traveled);
-    double dt = distance / speed;
-
-    // Save particle last state for tracklength tallies
-    this->time_last() = time();
-    this->r_last() = r();
-
-    // Advance particle in space and time
-    // Short-term solution until the surface source is revised and we can use
-    // this->move_distance(distance)
-    for (int j = 0; j < n_coord(); ++j) {
-      coord(j).r() += distance * coord(j).u();
-    }
-    this->time() += dt;
-    this->lifetime() += dt;
-    distance_traveled += distance;
-
-    // Score track-length tallies
-    if (!model::active_tracklength_tallies.empty()) {
-      score_tracklength_tally(*this, distance);
-    }
+  // Advance particle in space and time
+  // Short-term solution until the surface source is revised and we can use
+  // this->move_distance(distance)
+  for (int j = 0; j < n_coord(); ++j) {
+    coord(j).r() += distance * coord(j).u();
+  }
+  double dt = distance / speed;
+  this->time() += dt;
+  this->lifetime() += dt;
 
-    // Score flux derivative accumulators for differential tallies.
-    if (!model::active_tallies.empty()) {
-      score_track_derivative(*this, distance);
-    }
+  // Score track-length tallies
+  if (!model::active_tracklength_tallies.empty()) {
+    score_tracklength_tally(*this, distance);
   }
 
   // Score track-length estimate of k-eff
   if (settings::run_mode == RunMode::EIGENVALUE &&
       type() == ParticleType::neutron) {
-    keff_tally_tracklength() += wgt() * total_distance * macro_xs().nu_fission;
+    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
   }
 
-  // Set particle weight to zero if it hit the time cutoff boundary
-  if (total_distance == distance_cutoff)
+  // Score flux derivative accumulators for differential tallies.
+  if (!model::active_tallies.empty()) {
+    score_track_derivative(*this, distance);
+  }
+
+  // Set particle weight to zero if it hit the time boundary
+  if (distance == distance_cutoff) {
     wgt() = 0.0;
+  }
 }
 
 void Particle::event_cross_surface()
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 11f7c13c8f8..a96afedc725 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -58,6 +58,7 @@ vector<unique_ptr<Tally>> tallies;
 vector<int> active_tallies;
 vector<int> active_analog_tallies;
 vector<int> active_tracklength_tallies;
+vector<int> active_timed_tracklength_tallies;
 vector<int> active_collision_tallies;
 vector<int> active_meshsurf_tallies;
 vector<int> active_surface_tallies;
@@ -1097,6 +1098,7 @@ void setup_active_tallies()
   model::active_tallies.clear();
   model::active_analog_tallies.clear();
   model::active_tracklength_tallies.clear();
+  model::active_timed_tracklength_tallies.clear();
   model::active_collision_tallies.clear();
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
@@ -1108,6 +1110,9 @@ void setup_active_tallies()
 
     if (tally.active_) {
       model::active_tallies.push_back(i);
+      bool mesh_present = ((tally.get_filter<MeshFilter>() != nullptr) ||
+                           (tally.get_filter<MeshMaterialFilter>() != nullptr));
+      auto time_filter = tally.get_filter<TimeFilter>();
       switch (tally.type_) {
 
       case TallyType::VOLUME:
@@ -1116,12 +1121,11 @@ void setup_active_tallies()
           model::active_analog_tallies.push_back(i);
           break;
         case TallyEstimator::TRACKLENGTH:
-          model::active_tracklength_tallies.push_back(i);
-          if (auto time_filter = tally.get_filter<TimeFilter>()) {
-            if ((tally.get_filter<MeshFilter>() != nullptr) ||
-                (tally.get_filter<MeshMaterialFilter>() != nullptr)) {
-              model::add_to_time_grid(time_filter->bins());
-            }
+          if ((time_filter != nullptr) && mesh_present) {
+            model::active_timed_tracklength_tallies.push_back(i);
+            model::add_to_time_grid(time_filter->bins());
+          } else {
+            model::active_tracklength_tallies.push_back(i);
           }
           break;
         case TallyEstimator::COLLISION:
@@ -1158,6 +1162,7 @@ void free_memory_tally()
   model::active_tallies.clear();
   model::active_analog_tallies.clear();
   model::active_tracklength_tallies.clear();
+  model::active_timed_tracklength_tallies.clear();
   model::active_collision_tallies.clear();
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index e73fb90f311..859eec855f6 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2404,15 +2404,70 @@ void score_analog_tally_mg(Particle& p)
     match.bins_present_ = false;
 }
 
-void score_tracklength_tally(Particle& p, double distance)
+void score_tracklength_tally(Particle& p, double total_distance)
 {
+
+  if (!model::active_timed_tracklength_tallies.empty()) {
+    double speed = p.speed();
+    double total_dt = total_distance / speed;
+
+    // save particle last state
+    auto time_last = p.time_last();
+    auto r_last = p.r_last();
+
+    // move particle back
+    p.time() -= total_dt;
+    p.lifetime() -= total_dt;
+
+    for (int j = 0; j < p.n_coord(); ++j) {
+      p.coord(j).r() -= total_distance * p.coord(j).u();
+    }
+
+    double distance_traveled = 0.0;
+    while (distance_traveled < total_distance) {
+
+      double distance =
+        std::min(model::distance_to_time_boundary(p.time(), speed),
+          total_distance - distance_traveled);
+      double dt = distance / speed;
+
+      // Save particle last state for tracklength tallies
+      p.time_last() = p.time();
+      p.r_last() = p.r();
+
+      // Advance particle in space and time
+      // Short-term solution until the surface source is revised and we can use
+      // this->move_distance(distance)
+      for (int j = 0; j < p.n_coord(); ++j) {
+        p.coord(j).r() += distance * p.coord(j).u();
+      }
+      p.time() += dt;
+      p.lifetime() += dt;
+
+      // Determine the tracklength estimate of the flux
+      double flux = p.wgt() * distance;
+
+      score_tracklength_tally(p, flux, model::active_timed_tracklength_tallies);
+      distance_traveled += distance;
+    }
+
+    p.time_last() = time_last;
+    p.r_last() = r_last;
+  }
+
   // Determine the tracklength estimate of the flux
-  double flux = p.wgt() * distance;
+  double total_flux = p.wgt() * total_distance;
+
+  score_tracklength_tally(p, total_flux, model::active_tracklength_tallies);
+}
 
+void score_tracklength_tally(
+  Particle& p, double flux, const vector<int>& tallies)
+{
   // Set 'none' value for log union grid index
   int i_log_union = C_NONE;
 
-  for (auto i_tally : model::active_tracklength_tallies) {
+  for (auto i_tally : tallies) {
     const Tally& tally {*model::tallies[i_tally]};
 
     // Initialize an iterator over valid filter bin combinations.  If there are
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index e84c3d5caec..be4bc293485 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -653,48 +653,48 @@ def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Run the model; this should not cause an infinite loop
     model.run()
 
-
+    
 def test_filter_time_mesh(run_in_tmpdir):
     """
     Test that TimeFilter+MeshFilter with collision estimator agree with tracklength estimator.
     """
 
-    # ===========================================================================
-    # Set Material
-    # ===========================================================================
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    #Set Material
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
 
     mat = openmc.Material()
     mat.add_nuclide('Fe56', 1.0)
     mat.set_density('g/cm3', 7.8)
 
-    # ===========================================================================
-    # Set geometry
-    # ===========================================================================
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    #Set geometry
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
 
-    # Instantiate ZCylinder surfaces
+    #Instantiate ZCylinder surfaces
     surf_Z1 = openmc.XPlane(surface_id=1, x0=-1e10, boundary_type="reflective")
     surf_Z2 = openmc.XPlane(surface_id=2, x0=1e10, boundary_type="reflective")
 
-    # Instantiate Cells
+    #Instantiate Cells
     cell_F = openmc.Cell(cell_id=1, name="F")
 
-    # Use surface half-spaces to define regions
+    #Use surface half - spaces to define regions
     cell_F.region = +surf_Z1 & -surf_Z2
 
-    # Register Materials with Cells
+    #Register Materials with Cells
     cell_F.fill = mat
 
-    # Instantiate Universes
+    #Instantiate Universes
     root = openmc.Universe(universe_id=0, name="root universe", cells=[cell_F])
 
-    # Instantiate a Geometry, register the root Universe, and export to XML
+    #Instantiate a Geometry, register the root Universe, and export to XML
     geometry = openmc.Geometry(root)
 
-    # ===========================================================================
-    # Settings
-    # ===========================================================================
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    #Settings
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
 
-    # Instantiate a Settings object, set all runtime parameters, and export to XML
+    #Instantiate a Settings object, set all runtime parameters, and export to XML
     settings_file = openmc.Settings()
     settings_file.run_mode = "fixed source"
     settings_file.particles = 10000
@@ -702,16 +702,16 @@ def test_filter_time_mesh(run_in_tmpdir):
     settings_file.output = {"tallies": False}
     settings_file.cutoff = {"time_neutron": 1E-7}
 
-    # Create an initial uniform spatial source distribution over fissionable zones
+    #Create an initial uniform spatial source distribution over fissionable zones
     delta_dist = openmc.stats.Point()
     isotropic = openmc.stats.Isotropic()
     settings_file.source = openmc.IndependentSource(space=delta_dist, angle=isotropic)
 
-    # ===========================================================================
-    # Set tallies
-    # ===========================================================================
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    #Set tallies
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
 
-    # Create a mesh filter that can be used in a tally
+    #Create a mesh filter that can be used in a tally
     mesh = openmc.RegularMesh()
     mesh.dimension = (21, 1, 1)
     mesh.lower_left = (-20.5, -1e10, -1e10)
@@ -721,7 +721,7 @@ def test_filter_time_mesh(run_in_tmpdir):
     mesh_filter = openmc.MeshFilter(mesh)
     time_filter = openmc.TimeFilter(time_grid)
 
-    # Now use the mesh filter in a tally and indicate what scores are desired
+    #Now use the mesh filter in a tally and indicate what scores are desired
     tally1 = openmc.Tally(name="collision")
     tally1.estimator = "collision"
     tally1.filters = [time_filter, mesh_filter]
@@ -732,37 +732,36 @@ def test_filter_time_mesh(run_in_tmpdir):
     tally2.filters = [time_filter, mesh_filter]
     tally2.scores = ["flux"]
 
-    # Instantiate a tallies collection
+    #Instantiate a tallies collection
     tallies = openmc.Tallies([tally1, tally2])
-    
-    # ===========================================================================
-    # Set the model
-    # ===========================================================================
+
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    #Set the model
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
 
     model = openmc.Model()
     model.geometry = geometry
     model.settings = settings_file
     model.tallies = tallies
 
-    # ===========================================================================
-    # Run and post-process
-    # ===========================================================================
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    #Run and post - process
+    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
 
     sp_filename = model.run()
 
-    # Get radial flux distribution
+    #Get radial flux distribution
     with openmc.StatePoint(sp_filename) as sp:
         flux_collision = sp.tallies[tally1.id].mean.ravel()
         flux_collision_unc = sp.tallies[tally1.id].std_dev.ravel()
         flux_tracklength = sp.tallies[tally2.id].mean.ravel()
         flux_tracklength_unc = sp.tallies[tally2.id].std_dev.ravel()
 
-    # Construct arrays with uncertainties
+    #Construct arrays with uncertainties
     collision = uarray(flux_collision, flux_collision_unc)
     tracklength = uarray(flux_tracklength, flux_tracklength_unc)
     delta = collision - tracklength
 
-    # Check that difference is within uncertainty
+    #Check that difference is within uncertainty
     diff = nominal_values(delta)
     std_dev = std_devs(delta)
-    assert np.all(diff < 4*std_dev)

From 48ebdb78256677d10308cd3ba89a92ad693fa895 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Fri, 22 Aug 2025 17:14:20 +0200
Subject: [PATCH 04/10] simplify code by using move_distance

---
 src/tallies/tally_scoring.cpp | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 859eec855f6..30d5a34263d 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2416,13 +2416,10 @@ void score_tracklength_tally(Particle& p, double total_distance)
     auto r_last = p.r_last();
 
     // move particle back
+    p.move_distance(-total_distance);
     p.time() -= total_dt;
     p.lifetime() -= total_dt;
 
-    for (int j = 0; j < p.n_coord(); ++j) {
-      p.coord(j).r() -= total_distance * p.coord(j).u();
-    }
-
     double distance_traveled = 0.0;
     while (distance_traveled < total_distance) {
 
@@ -2436,11 +2433,7 @@ void score_tracklength_tally(Particle& p, double total_distance)
       p.r_last() = p.r();
 
       // Advance particle in space and time
-      // Short-term solution until the surface source is revised and we can use
-      // this->move_distance(distance)
-      for (int j = 0; j < p.n_coord(); ++j) {
-        p.coord(j).r() += distance * p.coord(j).u();
-      }
+      p.move_distance(distance);
       p.time() += dt;
       p.lifetime() += dt;
 

From fa2a7cc618eb9bac33135fc3ab2bb4e7978ad2e6 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Mon, 25 Aug 2025 13:32:07 +0200
Subject: [PATCH 05/10] fix when only timed tracklength tallies appear

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

diff --git a/src/particle.cpp b/src/particle.cpp
index fb41d82e950..da2d1e6855d 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -253,7 +253,8 @@ void Particle::event_advance()
   this->lifetime() += dt;
 
   // Score track-length tallies
-  if (!model::active_tracklength_tallies.empty()) {
+  if (!model::active_tracklength_tallies.empty() ||
+      !model::active_timed_tracklength_tallies.empty()) {
     score_tracklength_tally(*this, distance);
   }
 

From 987a1ece335bd052e525591b173cc0b8f583748d Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Mon, 25 Aug 2025 17:59:45 +0200
Subject: [PATCH 06/10] fix no assert in test

---
 tests/unit_tests/test_mesh.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index be4bc293485..f8748280765 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -765,3 +765,4 @@ def test_filter_time_mesh(run_in_tmpdir):
     #Check that difference is within uncertainty
     diff = nominal_values(delta)
     std_dev = std_devs(delta)
+    assert np.all(diff < 3*std_dev)

From 9174e53d742de9d4afb62d589509412d93f87ca2 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 26 Aug 2025 02:36:28 -0700
Subject: [PATCH 07/10] cosmetic fixes

---
 include/openmc/tallies/tally_scoring.h |   8 +-
 src/particle.cpp                       |   8 +-
 src/tallies/tally_scoring.cpp          | 104 +++++++++++++------------
 3 files changed, 62 insertions(+), 58 deletions(-)

diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
index 454a7db778f..4927b5227f9 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -91,17 +91,15 @@ void score_analog_tally_mg(Particle& p);
 //! \param distance The distance in [cm] traveled by the particle
 void score_tracklength_tally(Particle& p, double distance);
 
-//! Score tallies using a tracklength estimate of the flux.
+//! Score time filtered tallies using a tracklength estimate of the flux.
 //
 //! This is triggered at every event (surface crossing, lattice crossing, or
 //! collision) and thus cannot be done for tallies that require post-collision
 //! information.
 //
 //! \param p The particle being tracked
-//! \param flux The tracklength flux estimate in [cm].
-//! \param tallies A vector of the indices of the tallies to score to
-void score_tracklength_tally(
-  Particle& p, double flux, const vector<int>& tallies);
+//! \param distance The distance in [cm] traveled by the particle
+void score_timed_tracklength_tally(Particle& p, double distance);
 
 //! Score surface or mesh-surface tallies for particle currents.
 //
diff --git a/src/particle.cpp b/src/particle.cpp
index da2d1e6855d..9f273a2dd67 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -252,9 +252,13 @@ void Particle::event_advance()
   this->time() += dt;
   this->lifetime() += dt;
 
+  // Score timed track-length tallies
+  if (!model::active_timed_tracklength_tallies.empty()) {
+    score_timed_tracklength_tally(*this, distance);
+  }
+
   // Score track-length tallies
-  if (!model::active_tracklength_tallies.empty() ||
-      !model::active_timed_tracklength_tallies.empty()) {
+  if (!model::active_tracklength_tallies.empty()) {
     score_tracklength_tally(*this, distance);
   }
 
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 30d5a34263d..e8abd1b79a3 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2404,57 +2404,7 @@ void score_analog_tally_mg(Particle& p)
     match.bins_present_ = false;
 }
 
-void score_tracklength_tally(Particle& p, double total_distance)
-{
-
-  if (!model::active_timed_tracklength_tallies.empty()) {
-    double speed = p.speed();
-    double total_dt = total_distance / speed;
-
-    // save particle last state
-    auto time_last = p.time_last();
-    auto r_last = p.r_last();
-
-    // move particle back
-    p.move_distance(-total_distance);
-    p.time() -= total_dt;
-    p.lifetime() -= total_dt;
-
-    double distance_traveled = 0.0;
-    while (distance_traveled < total_distance) {
-
-      double distance =
-        std::min(model::distance_to_time_boundary(p.time(), speed),
-          total_distance - distance_traveled);
-      double dt = distance / speed;
-
-      // Save particle last state for tracklength tallies
-      p.time_last() = p.time();
-      p.r_last() = p.r();
-
-      // Advance particle in space and time
-      p.move_distance(distance);
-      p.time() += dt;
-      p.lifetime() += dt;
-
-      // Determine the tracklength estimate of the flux
-      double flux = p.wgt() * distance;
-
-      score_tracklength_tally(p, flux, model::active_timed_tracklength_tallies);
-      distance_traveled += distance;
-    }
-
-    p.time_last() = time_last;
-    p.r_last() = r_last;
-  }
-
-  // Determine the tracklength estimate of the flux
-  double total_flux = p.wgt() * total_distance;
-
-  score_tracklength_tally(p, total_flux, model::active_tracklength_tallies);
-}
-
-void score_tracklength_tally(
+void score_tracklength_tally_general(
   Particle& p, double flux, const vector<int>& tallies)
 {
   // Set 'none' value for log union grid index
@@ -2529,6 +2479,58 @@ void score_tracklength_tally(
     match.bins_present_ = false;
 }
 
+void score_timed_tracklength_tally(Particle& p, double total_distance)
+{
+  double speed = p.speed();
+  double total_dt = total_distance / speed;
+
+  // save particle last state
+  auto time_last = p.time_last();
+  auto r_last = p.r_last();
+
+  // move particle back
+  p.move_distance(-total_distance);
+  p.time() -= total_dt;
+  p.lifetime() -= total_dt;
+
+  double distance_traveled = 0.0;
+  while (distance_traveled < total_distance) {
+
+    double distance =
+      std::min(model::distance_to_time_boundary(p.time(), speed),
+        total_distance - distance_traveled);
+    double dt = distance / speed;
+
+    // Save particle last state for tracklength tallies
+    p.time_last() = p.time();
+    p.r_last() = p.r();
+
+    // Advance particle in space and time
+    p.move_distance(distance);
+    p.time() += dt;
+    p.lifetime() += dt;
+
+    // Determine the tracklength estimate of the flux
+    double flux = p.wgt() * distance;
+
+    score_tracklength_tally_general(p, flux, model::active_timed_tracklength_tallies);
+    distance_traveled += distance;
+  }
+
+  p.time_last() = time_last;
+  p.r_last() = r_last;
+
+}
+
+void score_tracklength_tally(Particle& p, double distance)
+{
+
+  // Determine the tracklength estimate of the flux
+  double flux = p.wgt() * distance;
+
+  score_tracklength_tally_general(p, flux, model::active_tracklength_tallies);
+}
+
 void score_collision_tally(Particle& p)
 {
   // Determine the collision estimate of the flux

From fa06444b169cc658603b68c2be8f3f363c48d183 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 26 Aug 2025 02:38:58 -0700
Subject: [PATCH 08/10] more cosmetic fixes

---
 include/openmc/tallies/tally_scoring.h | 4 ++--
 src/tallies/tally_scoring.cpp          | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
index 4927b5227f9..c3ab779e6a1 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -98,8 +98,8 @@ void score_tracklength_tally(Particle& p, double distance);
 //! information.
 //
 //! \param p The particle being tracked
-//! \param distance The distance in [cm] traveled by the particle
-void score_timed_tracklength_tally(Particle& p, double distance);
+//! \param total_distance The distance in [cm] traveled by the particle
+void score_timed_tracklength_tally(Particle& p, double total_distance);
 
 //! Score surface or mesh-surface tallies for particle currents.
 //
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index e8abd1b79a3..feb0235f2e3 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2513,13 +2513,13 @@ void score_timed_tracklength_tally(Particle& p, double total_distance)
     // Determine the tracklength estimate of the flux
     double flux = p.wgt() * distance;
 
-    score_tracklength_tally_general(p, flux, model::active_timed_tracklength_tallies);
+    score_tracklength_tally_general(
+      p, flux, model::active_timed_tracklength_tallies);
     distance_traveled += distance;
   }
 
   p.time_last() = time_last;
   p.r_last() = r_last;
-
 }
 
 void score_tracklength_tally(Particle& p, double distance)

From 3177de3591b8dd7fa96ca88876a96d97fc884dd9 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 12 Sep 2025 13:04:43 -0500
Subject: [PATCH 09/10] Update building of time_grid, style changes

---
 include/openmc/tallies/tally.h |   9 ++-
 src/tallies/tally.cpp          |  69 ++++++++++++---------
 src/tallies/tally_scoring.cpp  |   5 +-
 tests/unit_tests/test_mesh.py  | 107 +++++++++------------------------
 4 files changed, 78 insertions(+), 112 deletions(-)

diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 96687efcc3a..3beeb9d5ac5 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -211,8 +211,6 @@ extern vector<int> active_pulse_height_tallies;
 extern vector<int> pulse_height_cells;
 extern vector<double> time_grid;
 
-double distance_to_time_boundary(double time, double speed);
-
 } // namespace model
 
 namespace simulation {
@@ -244,6 +242,13 @@ void read_tallies_xml(pugi::xml_node root);
 //! batch to a new random variable
 void accumulate_tallies();
 
+//! Determine distance to next time boundary
+//
+//! \param time Current time of particle
+//! \param speed Speed of particle
+//! \return Distance to next time boundary (or INFTY if none)
+double distance_to_time_boundary(double time, double speed);
+
 //! Determine which tallies should be active
 void setup_active_tallies();
 
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index a96afedc725..ae0bffe6eb0 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -40,9 +40,10 @@
 #include "xtensor/xview.hpp"
 #include <fmt/core.h>
 
-#include <algorithm> // for max
+#include <algorithm> // for max, set_union
 #include <cassert>
-#include <cstddef> // for size_t
+#include <cstddef>  // for size_t
+#include <iterator> // for back_inserter
 #include <string>
 
 namespace openmc {
@@ -65,29 +66,6 @@ vector<int> active_surface_tallies;
 vector<int> active_pulse_height_tallies;
 vector<int> pulse_height_cells;
 vector<double> time_grid;
-
-double distance_to_time_boundary(double time, double speed)
-{
-  if (time_grid.size() == 0) {
-    return INFTY;
-  } else if (time >= time_grid[time_grid.size() - 1]) {
-    return INFTY;
-  } else {
-    return (*std::upper_bound(time_grid.begin(), time_grid.end(), time) -
-             time) *
-           speed;
-  }
-}
-
-void add_to_time_grid(vector<double> grid)
-{
-  auto temp = time_grid;
-  time_grid.resize(time_grid.size() + grid.size());
-  std::merge(temp.begin(), temp.end(), grid.begin(), grid.end(),
-    std::back_inserter(time_grid));
-  auto last_unique = std::unique(time_grid.begin(), time_grid.end());
-  time_grid.erase(last_unique, time_grid.end());
-}
 } // namespace model
 
 namespace simulation {
@@ -1093,6 +1071,39 @@ void accumulate_tallies()
   }
 }
 
+double distance_to_time_boundary(double time, double speed)
+{
+  if (model::time_grid.empty()) {
+    return INFTY;
+  } else if (time >= model::time_grid.back()) {
+    return INFTY;
+  } else {
+    double next_time =
+      *std::upper_bound(model::time_grid.begin(), model::time_grid.end(), time);
+    return (next_time - time) * speed;
+  }
+}
+
+//! Add new points to the global time grid
+//
+//! \param grid Vector of new time points to add
+void add_to_time_grid(vector<double> grid)
+{
+  if (grid.empty())
+    return;
+
+  // Create new vector with enough space to hold old and new grid points
+  vector<double> merged;
+  merged.reserve(model::time_grid.size() + grid.size());
+
+  // Merge and remove duplicates
+  std::set_union(model::time_grid.begin(), model::time_grid.end(), grid.begin(),
+    grid.end(), std::back_inserter(merged));
+
+  // Swap in the new grid
+  model::time_grid.swap(merged);
+}
+
 void setup_active_tallies()
 {
   model::active_tallies.clear();
@@ -1110,8 +1121,8 @@ void setup_active_tallies()
 
     if (tally.active_) {
       model::active_tallies.push_back(i);
-      bool mesh_present = ((tally.get_filter<MeshFilter>() != nullptr) ||
-                           (tally.get_filter<MeshMaterialFilter>() != nullptr));
+      bool mesh_present = (tally.get_filter<MeshFilter>() ||
+                           tally.get_filter<MeshMaterialFilter>());
       auto time_filter = tally.get_filter<TimeFilter>();
       switch (tally.type_) {
 
@@ -1121,9 +1132,9 @@ void setup_active_tallies()
           model::active_analog_tallies.push_back(i);
           break;
         case TallyEstimator::TRACKLENGTH:
-          if ((time_filter != nullptr) && mesh_present) {
+          if (time_filter && mesh_present) {
             model::active_timed_tracklength_tallies.push_back(i);
-            model::add_to_time_grid(time_filter->bins());
+            add_to_time_grid(time_filter->bins());
           } else {
             model::active_tracklength_tallies.push_back(i);
           }
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index feb0235f2e3..04b047b6497 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2496,9 +2496,8 @@ void score_timed_tracklength_tally(Particle& p, double total_distance)
   double distance_traveled = 0.0;
   while (distance_traveled < total_distance) {
 
-    double distance =
-      std::min(model::distance_to_time_boundary(p.time(), speed),
-        total_distance - distance_traveled);
+    double distance = std::min(distance_to_time_boundary(p.time(), speed),
+      total_distance - distance_traveled);
     double dt = distance / speed;
 
     // Save particle last state for tracklength tallies
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index f8748280765..8b6762d9cd6 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -653,116 +653,67 @@ def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Run the model; this should not cause an infinite loop
     model.run()
 
-    
-def test_filter_time_mesh(run_in_tmpdir):
-    """
-    Test that TimeFilter+MeshFilter with collision estimator agree with tracklength estimator.
-    """
 
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-    #Set Material
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+def test_filter_time_mesh(run_in_tmpdir):
+    """Test combination of TimeFilter and MeshFilter"""
 
+    # Define material
     mat = openmc.Material()
     mat.add_nuclide('Fe56', 1.0)
     mat.set_density('g/cm3', 7.8)
 
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-    #Set geometry
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-
-    #Instantiate ZCylinder surfaces
-    surf_Z1 = openmc.XPlane(surface_id=1, x0=-1e10, boundary_type="reflective")
-    surf_Z2 = openmc.XPlane(surface_id=2, x0=1e10, boundary_type="reflective")
-
-    #Instantiate Cells
-    cell_F = openmc.Cell(cell_id=1, name="F")
-
-    #Use surface half - spaces to define regions
-    cell_F.region = +surf_Z1 & -surf_Z2
-
-    #Register Materials with Cells
-    cell_F.fill = mat
-
-    #Instantiate Universes
-    root = openmc.Universe(universe_id=0, name="root universe", cells=[cell_F])
-
-    #Instantiate a Geometry, register the root Universe, and export to XML
-    geometry = openmc.Geometry(root)
-
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-    #Settings
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-
-    #Instantiate a Settings object, set all runtime parameters, and export to XML
-    settings_file = openmc.Settings()
-    settings_file.run_mode = "fixed source"
-    settings_file.particles = 10000
-    settings_file.batches = 20
-    settings_file.output = {"tallies": False}
-    settings_file.cutoff = {"time_neutron": 1E-7}
+    # Define geometry
+    surf_Z1 = openmc.XPlane(x0=-1e10, boundary_type="reflective")
+    surf_Z2 = openmc.XPlane(x0=1e10, boundary_type="reflective")
+    cell_F = openmc.Cell(fill=mat, region=+surf_Z1 & -surf_Z2)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell_F])
 
-    #Create an initial uniform spatial source distribution over fissionable zones
-    delta_dist = openmc.stats.Point()
-    isotropic = openmc.stats.Isotropic()
-    settings_file.source = openmc.IndependentSource(space=delta_dist, angle=isotropic)
+    # Define settings
+    model.settings.run_mode = "fixed source"
+    model.settings.particles = 1000
+    model.settings.batches = 20
+    model.settings.output = {"tallies": False}
+    model.settings.cutoff = {"time_neutron": 1e-7}
 
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-    #Set tallies
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
+    # Define tallies
 
-    #Create a mesh filter that can be used in a tally
+    # Create a mesh filter that can be used in a tally
     mesh = openmc.RegularMesh()
     mesh.dimension = (21, 1, 1)
     mesh.lower_left = (-20.5, -1e10, -1e10)
     mesh.upper_right = (20.5, 1e10, 1e10)
-    time_grid = np.linspace(0.0, 1E-7, 21)
+    time_grid = np.linspace(0.0, 1e-7, 21)
 
     mesh_filter = openmc.MeshFilter(mesh)
     time_filter = openmc.TimeFilter(time_grid)
 
-    #Now use the mesh filter in a tally and indicate what scores are desired
+    # Now use the mesh filter in a tally and indicate what scores are desired
     tally1 = openmc.Tally(name="collision")
     tally1.estimator = "collision"
     tally1.filters = [time_filter, mesh_filter]
     tally1.scores = ["flux"]
-
     tally2 = openmc.Tally(name="tracklength")
     tally2.estimator = "tracklength"
     tally2.filters = [time_filter, mesh_filter]
     tally2.scores = ["flux"]
+    model.tallies = openmc.Tallies([tally1, tally2])
 
-    #Instantiate a tallies collection
-    tallies = openmc.Tallies([tally1, tally2])
-
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-    #Set the model
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-
-    model = openmc.Model()
-    model.geometry = geometry
-    model.settings = settings_file
-    model.tallies = tallies
-
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-    #Run and post - process
-    #== == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == =
-
-    sp_filename = model.run()
+    # Run and post-process
+    model.run(apply_tally_results=True)
 
-    #Get radial flux distribution
-    with openmc.StatePoint(sp_filename) as sp:
-        flux_collision = sp.tallies[tally1.id].mean.ravel()
-        flux_collision_unc = sp.tallies[tally1.id].std_dev.ravel()
-        flux_tracklength = sp.tallies[tally2.id].mean.ravel()
-        flux_tracklength_unc = sp.tallies[tally2.id].std_dev.ravel()
+    # Get radial flux distribution
+    flux_collision = tally1.mean.ravel()
+    flux_collision_unc = tally1.std_dev.ravel()
+    flux_tracklength = tally2.mean.ravel()
+    flux_tracklength_unc = tally2.std_dev.ravel()
 
-    #Construct arrays with uncertainties
+    # Construct arrays with uncertainties
     collision = uarray(flux_collision, flux_collision_unc)
     tracklength = uarray(flux_tracklength, flux_tracklength_unc)
     delta = collision - tracklength
 
-    #Check that difference is within uncertainty
+    # Check that difference is within uncertainty
     diff = nominal_values(delta)
     std_dev = std_devs(delta)
     assert np.all(diff < 3*std_dev)

From d614d3d39af973ddd233779a31f5864abb889653 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 12 Sep 2025 13:16:05 -0500
Subject: [PATCH 10/10] Use chi-squared test instead of comparing std devs

---
 tests/unit_tests/test_mesh.py | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 8b6762d9cd6..f4e7bb8d797 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -3,6 +3,7 @@
 from pathlib import Path
 
 import numpy as np
+from scipy.stats import chi2
 import pytest
 import openmc
 import openmc.lib
@@ -482,7 +483,7 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
     np.testing.assert_almost_equal(mean, ref_data)
 
     # attempt to apply a dataset with an improper size to a VTK write
-    with pytest.raises(ValueError, match='Cannot apply dataset "mean"') as e:
+    with pytest.raises(ValueError, match='Cannot apply dataset "mean"'):
         simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename)
 
 def test_mesh_get_homogenized_materials():
@@ -713,7 +714,21 @@ def test_filter_time_mesh(run_in_tmpdir):
     tracklength = uarray(flux_tracklength, flux_tracklength_unc)
     delta = collision - tracklength
 
-    # Check that difference is within uncertainty
+    # Compute differences and standard deviations
     diff = nominal_values(delta)
     std_dev = std_devs(delta)
-    assert np.all(diff < 3*std_dev)
+
+    # Exclude zero-uncertainty bins
+    mask = std_dev > 0.0
+    dof = int(np.sum(mask))
+
+    # Global chi-square consistency test between collision and tracklength
+    # estimators. Target false positive rate ~1e-4 (1 in 10,000)
+    z = diff[mask] / std_dev[mask]
+    chi2_stat = np.sum(z * z)
+    alpha = 1.0e-4
+    crit = chi2.ppf(1 - alpha, dof)
+    assert chi2_stat < crit, (
+        f"Collision vs tracklength tallies disagree: chi2={chi2_stat:.2f} "
+        f">= {crit=:.2f} ({dof=}, {alpha=})"
+    )