Fixed a bug when combining TimeFilter, MeshFilter, and tracklength estimator (#3525)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

include/openmc/tallies/tally.h

@@ -203,11 +203,14 @@ 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;
 extern vector<int> active_pulse_height_tallies;
 extern vector<int> pulse_height_cells;
+extern vector<double> time_grid;
+
 } // namespace model
 
 namespace simulation {
@@ -239,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();
 

include/openmc/tallies/tally_scoring.h

@@ -91,6 +91,16 @@ 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 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 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.
 //
 //! \param p The particle being tracked

src/particle.cpp

@@ -252,6 +252,11 @@ 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()) {
     score_tracklength_tally(*this, distance);

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"
@@ -38,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 {
@@ -56,11 +59,13 @@ 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;
 vector<int> active_pulse_height_tallies;
 vector<int> pulse_height_cells;
+vector<double> time_grid;
 } // namespace model
 
 namespace simulation {
@@ -243,8 +248,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 +323,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 +499,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 +557,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 +991,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;
@@ -1064,21 +1071,59 @@ 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();
   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();
   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]};
 
     if (tally.active_) {
       model::active_tallies.push_back(i);
+      bool mesh_present = (tally.get_filter<MeshFilter>() ||
+                           tally.get_filter<MeshMaterialFilter>());
+      auto time_filter = tally.get_filter<TimeFilter>();
       switch (tally.type_) {
 
       case TallyType::VOLUME:
@@ -1087,7 +1132,12 @@ void setup_active_tallies()
           model::active_analog_tallies.push_back(i);
           break;
         case TallyEstimator::TRACKLENGTH:
-          model::active_tracklength_tallies.push_back(i);
+          if (time_filter && mesh_present) {
+            model::active_timed_tracklength_tallies.push_back(i);
+            add_to_time_grid(time_filter->bins());
+          } else {
+            model::active_tracklength_tallies.push_back(i);
+          }
           break;
         case TallyEstimator::COLLISION:
           model::active_collision_tallies.push_back(i);
@@ -1123,10 +1173,12 @@ 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();
   model::active_pulse_height_tallies.clear();
+  model::time_grid.clear();
 
   model::tally_map.clear();
 }
@@ -1465,8 +1517,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)
 {

src/tallies/tally_scoring.cpp

@@ -2404,15 +2404,13 @@ void score_analog_tally_mg(Particle& p)
     match.bins_present_ = false;
 }
 
-void score_tracklength_tally(Particle& p, double distance)
+void score_tracklength_tally_general(
+  Particle& p, double flux, const vector<int>& tallies)
 {
-  // Determine the tracklength estimate of the flux
-  double flux = p.wgt() * distance;
-
   // 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
@@ -2481,6 +2479,57 @@ void score_tracklength_tally(Particle& p, double distance)
     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(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

tests/unit_tests/test_mesh.py

@@ -3,10 +3,12 @@ from tempfile import TemporaryDirectory
 from pathlib import Path
 
 import numpy as np
+from scipy.stats import chi2
 import pytest
 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)])
@@ -481,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():
@@ -678,3 +680,82 @@ 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 combination of TimeFilter and MeshFilter"""
+
+    # Define material
+    mat = openmc.Material()
+    mat.add_nuclide('Fe56', 1.0)
+    mat.set_density('g/cm3', 7.8)
+
+    # 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])
+
+    # 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}
+
+    # Define 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"]
+    model.tallies = openmc.Tallies([tally1, tally2])
+
+    # Run and post-process
+    model.run(apply_tally_results=True)
+
+    # 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
+    collision = uarray(flux_collision, flux_collision_unc)
+    tracklength = uarray(flux_tracklength, flux_tracklength_unc)
+    delta = collision - tracklength
+
+    # Compute differences and standard deviations
+    diff = nominal_values(delta)
+    std_dev = std_devs(delta)
+
+    # 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=})"
+    )