From afd9d0607445f00c6de97b1ff866521c319722b0 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Fri, 12 Sep 2025 22:23:17 +0200 Subject: [PATCH] Fixed a bug when combining TimeFilter, MeshFilter, and tracklength estimator (#3525) Co-authored-by: Paul Romano --- include/openmc/tallies/tally.h | 10 ++++ include/openmc/tallies/tally_scoring.h | 10 ++++ src/particle.cpp | 5 ++ src/tallies/tally.cpp | 80 ++++++++++++++++++++----- src/tallies/tally_scoring.cpp | 59 ++++++++++++++++-- tests/unit_tests/test_mesh.py | 83 +++++++++++++++++++++++++- 6 files changed, 227 insertions(+), 20 deletions(-) diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 8c088c4608..3beeb9d5ac 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -203,11 +203,14 @@ extern vector> tallies; extern vector active_tallies; extern vector active_analog_tallies; extern vector active_tracklength_tallies; +extern vector active_timed_tracklength_tallies; extern vector active_collision_tallies; extern vector active_meshsurf_tallies; extern vector active_surface_tallies; extern vector active_pulse_height_tallies; extern vector pulse_height_cells; +extern vector 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(); diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 28f1f16223..c3ab779e6a 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/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 diff --git a/src/particle.cpp b/src/particle.cpp index fd4df14c3f..f25ada02c8 100644 --- a/src/particle.cpp +++ b/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); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ba18386026..ae0bffe6eb 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" @@ -38,9 +40,10 @@ #include "xtensor/xview.hpp" #include -#include // for max +#include // for max, set_union #include -#include // for size_t +#include // for size_t +#include // for back_inserter #include namespace openmc { @@ -56,11 +59,13 @@ vector> tallies; vector active_tallies; vector active_analog_tallies; vector active_tracklength_tallies; +vector active_timed_tracklength_tallies; vector active_collision_tallies; vector active_meshsurf_tallies; vector active_surface_tallies; vector active_pulse_height_tallies; vector pulse_height_cells; +vector 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& 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 grid) +{ + if (grid.empty()) + return; + + // Create new vector with enough space to hold old and new grid points + vector 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() || + tally.get_filter()); + auto time_filter = tally.get_filter(); 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) { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e73fb90f31..04b047b649 100644 --- a/src/tallies/tally_scoring.cpp +++ b/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& 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 diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 1930e3264e..f0f289408a 100644 --- a/tests/unit_tests/test_mesh.py +++ b/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=})" + )