From a4ae6717518bc8ce562f33abf80a60f13f985f67 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 13:23:10 +0000 Subject: [PATCH 1/6] Assert that cross section routines are called for the right particle Material::calculate_xs dispatches on the particle type, so the neutron and photon routines below it each have a particle type they assume. That assumption is not written down anywhere and is not checked, and one caller does not hold it: the micro cross section branch in score_tracklength_tally_general and its collision counterpart call Particle::update_neutron_xs for whatever particle is being scored, whenever the tallied nuclide is absent from the scoring region. For a photon that evaluates a neutron cross section at a photon's energy and stores it in a cache only neutron scoring ever reads, so the value is discarded and the only trace is the work done to produce it. Guard those two call sites on the particle type, and assert the assumption in the four routines that hold it: Material::calculate_neutron_xs, Material::calculate_photon_xs, Particle::update_neutron_xs and PhotonInteraction::calculate_xs. Note that under temperature interpolation the discarded evaluation consumed a pseudorandom number, since Nuclide::calculate_xs samples between the bracketing temperatures. Removing it shifts the random number stream, so a model combining photon transport, temperature interpolation and a nuclide bin tally with multiply_density off will no longer reproduce earlier results bit for bit. No test in the suite combines those. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01JxR1sQWWnALXj39rX12v7c --- src/material.cpp | 4 ++++ src/particle.cpp | 4 ++++ src/photon.cpp | 4 ++++ src/tallies/tally_scoring.cpp | 38 ++++++++++++++++++++++------------- 4 files changed, 36 insertions(+), 14 deletions(-) diff --git a/src/material.cpp b/src/material.cpp index 52cc0ec87bd..4e11c5e6843 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -827,6 +827,8 @@ void Material::calculate_xs(Particle& p) const void Material::calculate_neutron_xs(Particle& p) const { + assert(p.type().is_neutron()); + // Find energy index on energy grid int neutron = ParticleType::neutron().transport_index(); int i_grid = @@ -901,6 +903,8 @@ void Material::calculate_neutron_xs(Particle& p) const void Material::calculate_photon_xs(Particle& p) const { + assert(p.type().is_photon()); + p.macro_xs().coherent = 0.0; p.macro_xs().incoherent = 0.0; p.macro_xs().photoelectric = 0.0; diff --git a/src/particle.cpp b/src/particle.cpp index 198da003a79..ac15d394695 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,5 +1,7 @@ #include "openmc/particle.h" +#include + #include // copy, min #include // log, abs @@ -964,6 +966,8 @@ void Particle::write_restart() const void Particle::update_neutron_xs( int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs) { + assert(type().is_neutron()); + // Get microscopic cross section cache auto& micro = this->neutron_xs(i_nuclide); diff --git a/src/photon.cpp b/src/photon.cpp index 72c970bc15f..b51665695e4 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -1,5 +1,7 @@ #include "openmc/photon.h" +#include + #include "openmc/array.h" #include "openmc/bremsstrahlung.h" #include "openmc/constants.h" @@ -710,6 +712,8 @@ void PhotonInteraction::compton_doppler( void PhotonInteraction::calculate_xs(Particle& p) const { + assert(p.type().is_photon()); + // Perform binary search on the element energy grid in order to determine // which points to interpolate between int n_grid = energy_.size(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 27e6671833a..93c63592eb0 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2447,14 +2447,19 @@ void score_tracklength_tally_general( if (!tally.multiply_density()) atom_density = 1.0; } else if (!tally.multiply_density()) { - // Determine log union grid index - if (i_log_union == C_NONE) { - int neutron = ParticleType::neutron().transport_index(); - i_log_union = std::log(p.E() / data::energy_min[neutron]) / - simulation::log_spacing; + // The cache being filled here holds neutron cross sections, and is + // only ever read back by neutron scoring, so there is nothing to + // update for any other kind of particle + if (p.type().is_neutron()) { + // Determine log union grid index + if (i_log_union == C_NONE) { + int neutron = ParticleType::neutron().transport_index(); + i_log_union = std::log(p.E() / data::energy_min[neutron]) / + simulation::log_spacing; + } + // Update micro xs cache + p.update_neutron_xs(i_nuclide, i_log_union); } - // Update micro xs cache - p.update_neutron_xs(i_nuclide, i_log_union); atom_density = 1.0; } } @@ -2577,14 +2582,19 @@ void score_collision_tally(Particle& p) if (!tally.multiply_density()) atom_density = 1.0; } else if (!tally.multiply_density()) { - // Determine log union grid index - if (i_log_union == C_NONE) { - int neutron = ParticleType::neutron().transport_index(); - i_log_union = std::log(p.E() / data::energy_min[neutron]) / - simulation::log_spacing; + // The cache being filled here holds neutron cross sections, and is + // only ever read back by neutron scoring, so there is nothing to + // update for any other kind of particle + if (p.type().is_neutron()) { + // Determine log union grid index + if (i_log_union == C_NONE) { + int neutron = ParticleType::neutron().transport_index(); + i_log_union = std::log(p.E() / data::energy_min[neutron]) / + simulation::log_spacing; + } + // Update micro xs cache + p.update_neutron_xs(i_nuclide, i_log_union); } - // Update micro xs cache - p.update_neutron_xs(i_nuclide, i_log_union); atom_density = 1.0; } } From e94fcbdade107903a3e8887e02546a752fb9a4e5 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 13:52:27 +0000 Subject: [PATCH 2/6] Keep the micro cross section caches addressed correctly Each particle carries two micro cross section caches: neutron_xs_, sized by data::nuclides, and photon_xs_, sized by data::elements because photon data does not distinguish isotopes. Which particle type each applies to, and which index space addresses it, are assumptions nothing states or checks, and three places get them wrong. Tally scoring indexes the photon cache with a nuclide index. The two index spaces coincide only when every element contributes exactly one nuclide, so a material of U235, U238 and O16 scores oxygen's cross section for U238 and reads past the end of the vector for O16. A tally may also name a nuclide no material contains, which adds a nuclide without adding an element, so even a single nuclide material reaches out of bounds. Add data::nuclide_to_element, recorded where the element is already being resolved, and resolve the nuclide once at the top of score_general_ce_nonanalog, where the particle and the nuclide are fixed for the whole call. The four lookups then address the cache by element index the way every other caller in the tree already does, so there is no second way to reach it and nothing new to get wrong. The two accessors name their parameters i_element and i_nuclide rather than i. The branch that refreshes a stale cache when the tallied nuclide is absent from the scoring region only handles neutrons, so a photon scored whatever energy it last collided at. Give it Particle::update_photon_xs, which refreshes an element's entry unless it already holds the particle's energy. It takes an element index rather than a nuclide one, unlike update_neutron_xs next to it, because that is the index space photon data is tabulated in; the caller that has a nuclide resolves it the same way scoring does. That same branch called Particle::update_neutron_xs for whatever particle was being scored. For a photon that evaluated a neutron cross section at a photon's energy and stored it where only neutron scoring reads, so the value was discarded. Guard it on the particle type, and assert the particle type in the four routines that assume one: Material::calculate_neutron_xs, Material::calculate_photon_xs, Particle::update_neutron_xs and PhotonInteraction::calculate_xs. The two estimators carried identical copies of that branch, so the guard would have doubled it. Both now call one update_absent_nuclide_xs, which is shorter than what either of them started with. Note that under temperature interpolation the discarded neutron evaluation consumed a pseudorandom number, since Nuclide::calculate_xs samples between the bracketing temperatures. Removing it shifts the random number stream, so a model combining photon transport, temperature interpolation and a nuclide bin tally with multiply_density off will not reproduce earlier results bit for bit. The absent nuclide test needs multiply_density off. An absent nuclide has an atom density of zero, so with it on the score is zero whatever the cache held, which is both a test that cannot fail for the right reason and the branch that refreshes the cache in the first place. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01JxR1sQWWnALXj39rX12v7c --- include/openmc/particle.h | 9 +++ include/openmc/particle_data.h | 9 ++- include/openmc/photon.h | 8 +++ src/nuclide.cpp | 8 +++ src/particle.cpp | 10 +++ src/photon.cpp | 2 + src/tallies/tally_scoring.cpp | 72 +++++++++++-------- tests/unit_tests/test_photon_micro_xs.py | 88 ++++++++++++++++++++++++ 8 files changed, 173 insertions(+), 33 deletions(-) create mode 100644 tests/unit_tests/test_photon_micro_xs.py diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 34ffc18cf74..ba52070fe61 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -125,6 +125,15 @@ class Particle : public ParticleData { //! \param[in] ncrystal_xs Thermal scattering xs from NCrystal void update_neutron_xs(int i_nuclide, int i_grid = C_NONE, int i_sab = C_NONE, double sab_frac = 0.0, double ncrystal_xs = -1.0); + + //! Update the microscopic cross section cache of an element + // + //! Evaluates the element's cross sections at the particle's energy unless + //! the cache already holds them, which is what transport leaves behind for + //! the elements present where the particle last collided. + //! + //! \param[in] i_element Index in data::elements + void update_photon_xs(int i_element); }; //============================================================================ diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 24b7eb53c0a..3dc38a39659 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -584,11 +584,14 @@ class ParticleData : public GeometryState { // Cross section caches // Microscopic neutron cross sections - NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; } - const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; } + NuclideMicroXS& neutron_xs(int i_nuclide) { return neutron_xs_[i_nuclide]; } + const NuclideMicroXS& neutron_xs(int i_nuclide) const + { + return neutron_xs_[i_nuclide]; + } // Microscopic photon cross sections - ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; } + ElementMicroXS& photon_xs(int i_element) { return photon_xs_[i_element]; } // Macroscopic cross sections MacroXS& macro_xs() { return macro_xs_; } diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 9c65b482ba0..cbba70d90d2 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -185,6 +185,14 @@ extern tensor::Tensor extern std::unordered_map element_map; extern vector> elements; +//! Index in \ref elements of the element each nuclide belongs to +// +//! Photon data is tabulated per element rather than per nuclide, so the photon +//! micro cross section cache is indexed by element while the neutron one is +//! indexed by nuclide. Anything holding a nuclide index needs this to reach the +//! photon data. Entries are C_NONE when photon transport is off. +extern vector nuclide_to_element; + } // namespace data } // namespace openmc diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 17d6e952c39..e3d3084f5fb 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1168,6 +1168,14 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n) close_group(group); file_close(file_id); } + + // Record which element this nuclide belongs to, so that a nuclide index + // can reach the photon data it shares with the other isotopes + int i_nuclide = data::nuclide_map.at(name); + if (data::nuclide_to_element.size() <= i_nuclide) { + data::nuclide_to_element.resize(i_nuclide + 1, C_NONE); + } + data::nuclide_to_element[i_nuclide] = data::element_map.at(element); } } return 0; diff --git a/src/particle.cpp b/src/particle.cpp index ac15d394695..2f627bdeb62 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -986,6 +986,16 @@ void Particle::update_neutron_xs( } } +void Particle::update_photon_xs(int i_element) +{ + assert(type().is_photon()); + + // If the cache doesn't match, recalculate micro xs + if (this->E() != this->photon_xs(i_element).last_E) { + data::elements[i_element]->calculate_xs(*this); + } +} + //============================================================================== // Non-method functions //============================================================================== diff --git a/src/photon.cpp b/src/photon.cpp index b51665695e4..3674d69f787 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -38,6 +38,7 @@ tensor::Tensor compton_profile_pz; std::unordered_map element_map; vector> elements; +vector nuclide_to_element; } // namespace data @@ -1120,6 +1121,7 @@ std::pair klein_nishina(double alpha, uint64_t* seed) void free_memory_photon() { data::elements.clear(); + data::nuclide_to_element.clear(); data::compton_profile_pz.resize({0}); data::ttb_e_grid.resize({0}); data::ttb_k_grid.resize({0}); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 93c63592eb0..1972c189593 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -177,6 +177,33 @@ void score_fission_delayed_dg(int i_tally, int d_bin, double score, dg_match.bins_[i_bin] = original_bin; } +//! Helper function to refresh the micro xs cache of an absent nuclide +// +//! A tally binned by nuclide with multiply_density off scores a nuclide the +//! scoring region need not contain, and transport only evaluates the cross +//! sections of what is actually there. Each cache is read back solely by +//! scoring of the particle type it holds cross sections for, so there is +//! nothing to refresh for any other kind of particle. +// +//! \param[in,out] p Particle being scored +//! \param[in] i_nuclide Index in data::nuclides +//! \param[in,out] i_log_union Index on the log union grid, or C_NONE if it has +//! not been determined yet +void update_absent_nuclide_xs(Particle& p, int i_nuclide, int& i_log_union) +{ + if (p.type().is_neutron()) { + // Determine log union grid index + if (i_log_union == C_NONE) { + int neutron = ParticleType::neutron().transport_index(); + i_log_union = + std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing; + } + p.update_neutron_xs(i_nuclide, i_log_union); + } else if (p.type().is_photon()) { + p.update_photon_xs(data::nuclide_to_element[i_nuclide]); + } +} + //! Helper function to retrieve fission q value from a nuclide double get_nuc_fission_q(const Nuclide& nuc, const Particle& p, int score_bin) @@ -586,6 +613,13 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, // Get the pre-collision energy of the particle. auto E = p.E_last(); + // Photon data is tabulated per element, so the cache holding photon micro + // cross sections is indexed by element while i_nuclide is not. Resolve it + // once here; C_NONE whenever no score below can ask for one. + int i_element = (i_nuclide >= 0 && p.type().is_photon()) + ? data::nuclide_to_element[i_nuclide] + : C_NONE; + for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; auto score_index = start_index + i; @@ -601,7 +635,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (p.type().is_neutron()) { score = p.neutron_xs(i_nuclide).total * atom_density * flux; } else if (p.type().is_photon()) { - score = p.photon_xs(i_nuclide).total * atom_density * flux; + score = p.photon_xs(i_element).total * atom_density * flux; } } else { score = p.macro_xs().total * flux; @@ -625,7 +659,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const auto& micro = p.neutron_xs(i_nuclide); score = (micro.total - micro.absorption) * atom_density * flux; } else { - const auto& micro = p.photon_xs(i_nuclide); + const auto& micro = p.photon_xs(i_element); score = (micro.coherent + micro.incoherent) * atom_density * flux; } } else { @@ -645,7 +679,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (p.type().is_neutron()) { score = p.neutron_xs(i_nuclide).absorption * atom_density * flux; } else { - const auto& xs = p.photon_xs(i_nuclide); + const auto& xs = p.photon_xs(i_element); score = (xs.total - xs.coherent - xs.incoherent) * atom_density * flux; } @@ -1052,7 +1086,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, continue; if (i_nuclide >= 0) { - const auto& micro = p.photon_xs(i_nuclide); + const auto& micro = p.photon_xs(i_element); double xs = (score_bin == COHERENT) ? micro.coherent : (score_bin == INCOHERENT) ? micro.incoherent : (score_bin == PHOTOELECTRIC) ? micro.photoelectric @@ -2447,19 +2481,8 @@ void score_tracklength_tally_general( if (!tally.multiply_density()) atom_density = 1.0; } else if (!tally.multiply_density()) { - // The cache being filled here holds neutron cross sections, and is - // only ever read back by neutron scoring, so there is nothing to - // update for any other kind of particle - if (p.type().is_neutron()) { - // Determine log union grid index - if (i_log_union == C_NONE) { - int neutron = ParticleType::neutron().transport_index(); - i_log_union = std::log(p.E() / data::energy_min[neutron]) / - simulation::log_spacing; - } - // Update micro xs cache - p.update_neutron_xs(i_nuclide, i_log_union); - } + // Update micro xs cache + update_absent_nuclide_xs(p, i_nuclide, i_log_union); atom_density = 1.0; } } @@ -2582,19 +2605,8 @@ void score_collision_tally(Particle& p) if (!tally.multiply_density()) atom_density = 1.0; } else if (!tally.multiply_density()) { - // The cache being filled here holds neutron cross sections, and is - // only ever read back by neutron scoring, so there is nothing to - // update for any other kind of particle - if (p.type().is_neutron()) { - // Determine log union grid index - if (i_log_union == C_NONE) { - int neutron = ParticleType::neutron().transport_index(); - i_log_union = std::log(p.E() / data::energy_min[neutron]) / - simulation::log_spacing; - } - // Update micro xs cache - p.update_neutron_xs(i_nuclide, i_log_union); - } + // Update micro xs cache + update_absent_nuclide_xs(p, i_nuclide, i_log_union); atom_density = 1.0; } } diff --git a/tests/unit_tests/test_photon_micro_xs.py b/tests/unit_tests/test_photon_micro_xs.py new file mode 100644 index 00000000000..e46c1ed07e0 --- /dev/null +++ b/tests/unit_tests/test_photon_micro_xs.py @@ -0,0 +1,88 @@ +"""Photon micro cross section tallies are looked up per element. + +Photon data does not distinguish isotopes, so the cache holding photon micro +cross sections is indexed by element while a tally's nuclide bins are indexed by +nuclide. Two isotopes of the same element present at the same atom density must +therefore score the same microscopic cross section, whatever that cross section +happens to be -- which makes this check independent of the data library. +""" + +import openmc +import pytest + + +def test_isotopes_of_one_element_score_alike(run_in_tmpdir): + openmc.reset_auto_ids() + + mat = openmc.Material() + mat.add_nuclide('U235', 1.0) + mat.add_nuclide('U238', 1.0) + mat.set_density('g/cm3', 10.0) + + sphere = openmc.Sphere(r=5.0, boundary_type='vacuum') + model = openmc.Model() + model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-sphere)]) + model.settings.run_mode = 'fixed source' + model.settings.photon_transport = True + model.settings.particles = 200 + model.settings.batches = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1.0e6], [1.0]), + particle='photon') + + tally = openmc.Tally() + tally.nuclides = ['U235', 'U238'] + tally.scores = ['total'] + model.tallies = openmc.Tallies([tally]) + + model.run(apply_tally_results=True) + + u235, u238 = tally.mean.ravel() + assert u235 > 0.0 + assert u238 == pytest.approx(u235, rel=1e-12) + + +def test_nuclide_absent_from_the_material(run_in_tmpdir): + """A tally may name a nuclide no material contains. + + Doing so adds it to the global nuclide list without adding an element, so + its index reaches past the end of the element-indexed photon cache unless + the lookup goes through the nuclide-to-element mapping. + + multiply_density has to be off for the score to say anything: an absent + nuclide has an atom density of zero, which multiplies away whatever came + out of the cache. With it off the score is the microscopic cross section + itself, and it is also the branch that refreshes the cache for a nuclide + the scoring region does not contain. + """ + openmc.reset_auto_ids() + + mat = openmc.Material() + mat.add_nuclide('U235', 1.0) + mat.set_density('g/cm3', 10.0) + + sphere = openmc.Sphere(r=5.0, boundary_type='vacuum') + model = openmc.Model() + model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-sphere)]) + model.settings.run_mode = 'fixed source' + model.settings.photon_transport = True + model.settings.particles = 200 + model.settings.batches = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1.0e6], [1.0]), + particle='photon') + + tally = openmc.Tally() + tally.nuclides = ['U235', 'U238'] # U238 is not in the material + tally.scores = ['total'] + tally.multiply_density = False + model.tallies = openmc.Tallies([tally]) + + model.run(apply_tally_results=True) + + # Same element as U235, so the same microscopic cross section + u235, u238 = tally.mean.ravel() + assert u235 > 0.0 + assert u238 == pytest.approx(u235, rel=1e-12) From 80f02fd05decb3b170b631a2eda491947fe6763c Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 6 Oct 2026 17:35:25 +0300 Subject: [PATCH 3/6] Apply batched suggestions from code review Co-authored-by: Jonathan Shimwell Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com> --- include/openmc/photon.h | 2 +- src/particle.cpp | 5 +- src/photon.cpp | 2 +- tests/unit_tests/test_photon_micro_xs.py | 69 ++++++++++++++++++++++-- 4 files changed, 70 insertions(+), 8 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index cbba70d90d2..8f15b62209e 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -190,7 +190,7 @@ extern vector> elements; //! Photon data is tabulated per element rather than per nuclide, so the photon //! micro cross section cache is indexed by element while the neutron one is //! indexed by nuclide. Anything holding a nuclide index needs this to reach the -//! photon data. Entries are C_NONE when photon transport is off. +//! photon data. Empty when photon transport is off. extern vector nuclide_to_element; } // namespace data diff --git a/src/particle.cpp b/src/particle.cpp index 2f627bdeb62..5d9c7a02510 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,9 +1,8 @@ #include "openmc/particle.h" -#include - #include // copy, min -#include // log, abs +#include +#include // log, abs #include diff --git a/src/photon.cpp b/src/photon.cpp index 3674d69f787..ab16dcafcad 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -1,6 +1,5 @@ #include "openmc/photon.h" -#include #include "openmc/array.h" #include "openmc/bremsstrahlung.h" @@ -18,6 +17,7 @@ #include "openmc/tensor.h" +#include #include #include #include diff --git a/tests/unit_tests/test_photon_micro_xs.py b/tests/unit_tests/test_photon_micro_xs.py index e46c1ed07e0..1cef989c30a 100644 --- a/tests/unit_tests/test_photon_micro_xs.py +++ b/tests/unit_tests/test_photon_micro_xs.py @@ -4,7 +4,7 @@ cross sections is indexed by element while a tally's nuclide bins are indexed by nuclide. Two isotopes of the same element present at the same atom density must therefore score the same microscopic cross section, whatever that cross section -happens to be -- which makes this check independent of the data library. +happens to be, which makes this check independent of the data library. """ import openmc @@ -53,8 +53,9 @@ def test_nuclide_absent_from_the_material(run_in_tmpdir): multiply_density has to be off for the score to say anything: an absent nuclide has an atom density of zero, which multiplies away whatever came out of the cache. With it off the score is the microscopic cross section - itself, and it is also the branch that refreshes the cache for a nuclide - the scoring region does not contain. + itself. The element is still present in the material, so transport keeps + its cache current and this does not exercise the refresh for an element + the scoring region lacks. test_element_absent_from_the_scoring_region does. """ openmc.reset_auto_ids() @@ -86,3 +87,65 @@ def test_nuclide_absent_from_the_material(run_in_tmpdir): u235, u238 = tally.mean.ravel() assert u235 > 0.0 assert u238 == pytest.approx(u235, rel=1e-12) + + +def test_element_absent_from_the_scoring_region(run_in_tmpdir): + """The cache of an element absent from the scoring region is refreshed. + + In the test above the tallied nuclide shares its element with the material, + so transport keeps that element's cache current and the refresh is never + needed. Here the inner sphere holds no silicon at all. Uncollided photons + crossing it are all at the source energy and have not yet reached the + silicon shell, so nothing but the tally evaluates silicon for them. The + microscopic total over the flux in the uncollided energy bin is then the + silicon cross section at the source energy, and must match the same ratio + in the shell, which does contain silicon. + """ + openmc.reset_auto_ids() + + iron = openmc.Material() + iron.add_nuclide('Fe56', 1.0) + iron.set_density('g/cm3', 0.1) + + silicon = openmc.Material() + silicon.add_nuclide('Si28', 1.0) + silicon.set_density('g/cm3', 0.1) + + inner = openmc.Sphere(r=5.0) + outer = openmc.Sphere(r=10.0, boundary_type='vacuum') + iron_cell = openmc.Cell(fill=iron, region=-inner) + silicon_cell = openmc.Cell(fill=silicon, region=+inner & -outer) + + model = openmc.Model() + model.geometry = openmc.Geometry([iron_cell, silicon_cell]) + model.settings.run_mode = 'fixed source' + model.settings.photon_transport = True + model.settings.particles = 200 + model.settings.batches = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1.0e6], [1.0]), + particle='photon') + + filters = [ + openmc.CellFilter([iron_cell, silicon_cell]), + openmc.EnergyFilter([0.999e6, 1.001e6]), + ] + flux = openmc.Tally() + flux.filters = filters + flux.scores = ['flux'] + + total = openmc.Tally() + total.filters = filters + total.nuclides = ['Si28'] + total.scores = ['total'] + total.multiply_density = False + model.tallies = openmc.Tallies([flux, total]) + + model.run(apply_tally_results=True) + + iron_flux, silicon_flux = flux.mean.ravel() + iron_total, silicon_total = total.mean.ravel() + assert silicon_total > 0.0 + assert iron_total / iron_flux == pytest.approx( + silicon_total / silicon_flux, rel=1e-10) From c69b5d96a4f361d90150f246f9233f3431c8d480 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 6 Oct 2026 17:47:48 +0300 Subject: [PATCH 4/6] clang format --- src/photon.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/photon.cpp b/src/photon.cpp index ab16dcafcad..11e020206d7 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -1,6 +1,5 @@ #include "openmc/photon.h" - #include "openmc/array.h" #include "openmc/bremsstrahlung.h" #include "openmc/constants.h" From 553489d9b714eaa67780255209feb7803e4b5b2c Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 6 Oct 2026 18:16:17 +0200 Subject: [PATCH 5/6] Test the absent element refresh with both estimators Run test_element_absent_from_the_scoring_region with the collision estimator as well as tracklength. The collision estimator scores after the collision, when photoelectric absorption has set the photon energy to zero, so refreshing the cache at p.E() gives NaN. The source is at 100 keV so photoelectric absorption in the iron is common enough to catch it. Also drop a redeclaration of i_nuclide in openmc_load_nuclide that shadowed the identical one a few lines above. Co-authored-by: jon-proximafusion <174589458+jon-proximafusion@users.noreply.github.com> --- src/nuclide.cpp | 1 - tests/unit_tests/test_photon_micro_xs.py | 14 +++++++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index e3d3084f5fb..f3e0420ab4c 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1171,7 +1171,6 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n) // Record which element this nuclide belongs to, so that a nuclide index // can reach the photon data it shares with the other isotopes - int i_nuclide = data::nuclide_map.at(name); if (data::nuclide_to_element.size() <= i_nuclide) { data::nuclide_to_element.resize(i_nuclide + 1, C_NONE); } diff --git a/tests/unit_tests/test_photon_micro_xs.py b/tests/unit_tests/test_photon_micro_xs.py index 1cef989c30a..bab69043197 100644 --- a/tests/unit_tests/test_photon_micro_xs.py +++ b/tests/unit_tests/test_photon_micro_xs.py @@ -89,7 +89,8 @@ def test_nuclide_absent_from_the_material(run_in_tmpdir): assert u238 == pytest.approx(u235, rel=1e-12) -def test_element_absent_from_the_scoring_region(run_in_tmpdir): +@pytest.mark.parametrize('estimator', ['tracklength', 'collision']) +def test_element_absent_from_the_scoring_region(run_in_tmpdir, estimator): """The cache of an element absent from the scoring region is refreshed. In the test above the tallied nuclide shares its element with the material, @@ -100,6 +101,11 @@ def test_element_absent_from_the_scoring_region(run_in_tmpdir): microscopic total over the flux in the uncollided energy bin is then the silicon cross section at the source energy, and must match the same ratio in the shell, which does contain silicon. + + The collision estimator scores after the collision, when photoelectric + absorption has already set the photon energy to zero, so the refresh has to + use the pre-collision energy. The source is at 100 keV so that photoelectric + absorption in the iron is common enough to catch that. """ openmc.reset_auto_ids() @@ -124,22 +130,24 @@ def test_element_absent_from_the_scoring_region(run_in_tmpdir): model.settings.batches = 2 model.settings.source = openmc.IndependentSource( space=openmc.stats.Point(), - energy=openmc.stats.Discrete([1.0e6], [1.0]), + energy=openmc.stats.Discrete([1.0e5], [1.0]), particle='photon') filters = [ openmc.CellFilter([iron_cell, silicon_cell]), - openmc.EnergyFilter([0.999e6, 1.001e6]), + openmc.EnergyFilter([0.999e5, 1.001e5]), ] flux = openmc.Tally() flux.filters = filters flux.scores = ['flux'] + flux.estimator = estimator total = openmc.Tally() total.filters = filters total.nuclides = ['Si28'] total.scores = ['total'] total.multiply_density = False + total.estimator = estimator model.tallies = openmc.Tallies([flux, total]) model.run(apply_tally_results=True) From ec63c25be229b1e3ae99a86b478ea2b4170a47fe Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 6 Oct 2026 23:29:37 +0300 Subject: [PATCH 6/6] fix collision estimator --- src/tallies/tally_scoring.cpp | 16 ++++- tests/unit_tests/test_neutron_micro_xs.py | 76 +++++++++++++++++++++++ tests/unit_tests/test_photon_micro_xs.py | 6 +- 3 files changed, 93 insertions(+), 5 deletions(-) create mode 100644 tests/unit_tests/test_neutron_micro_xs.py diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 1972c189593..d808a315b8e 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -185,12 +185,22 @@ void score_fission_delayed_dg(int i_tally, int d_bin, double score, //! scoring of the particle type it holds cross sections for, so there is //! nothing to refresh for any other kind of particle. // +//! Scores are evaluated at the pre-collision energy, so the cross sections are +//! too. The collision estimator runs after the collision, when p.E() is +//! already the outgoing energy (zero after photoelectric absorption or pair +//! production), and the cross section routines read p.E() throughout, so the +//! particle is evaluated at E_last() and its energy restored afterwards. On the +//! tracklength path the two are equal. +// //! \param[in,out] p Particle being scored //! \param[in] i_nuclide Index in data::nuclides -//! \param[in,out] i_log_union Index on the log union grid, or C_NONE if it has -//! not been determined yet +//! \param[in,out] i_log_union Index on the log union grid at E_last(), or +//! C_NONE if it has not been determined yet void update_absent_nuclide_xs(Particle& p, int i_nuclide, int& i_log_union) { + double E = p.E(); + p.E() = p.E_last(); + if (p.type().is_neutron()) { // Determine log union grid index if (i_log_union == C_NONE) { @@ -202,6 +212,8 @@ void update_absent_nuclide_xs(Particle& p, int i_nuclide, int& i_log_union) } else if (p.type().is_photon()) { p.update_photon_xs(data::nuclide_to_element[i_nuclide]); } + + p.E() = E; } //! Helper function to retrieve fission q value from a nuclide diff --git a/tests/unit_tests/test_neutron_micro_xs.py b/tests/unit_tests/test_neutron_micro_xs.py new file mode 100644 index 00000000000..2496b130292 --- /dev/null +++ b/tests/unit_tests/test_neutron_micro_xs.py @@ -0,0 +1,76 @@ +"""Neutron micro cross sections of a nuclide absent from the scoring region. + +A tally binned by nuclide with multiply_density off scores nuclides the +scoring region need not contain, and the cross sections of those are evaluated +at scoring time rather than by transport. This checks that they are evaluated +at the energy the score belongs to, by comparing against a region that does +contain the nuclide, which makes the check independent of the data library. +""" + +import openmc +import pytest + + +@pytest.mark.parametrize('estimator', ['tracklength', 'collision']) +def test_nuclide_absent_from_the_scoring_region(run_in_tmpdir, estimator): + """The cache of a nuclide absent from the scoring region is refreshed. + + The inner sphere holds no silicon at all, so nothing but the tally + evaluates silicon there. The energy bin is narrow enough that only + neutrons at the source energy fall in it, so the microscopic total over + the flux in that bin is the silicon cross section at the source energy, + and must match the same ratio in the shell, which does contain silicon. + + The collision estimator scores after the collision, when elastic + scattering off iron has already changed the neutron energy, so the + refresh has to use the pre-collision energy. + """ + openmc.reset_auto_ids() + + iron = openmc.Material() + iron.add_nuclide('Fe56', 1.0) + iron.set_density('g/cm3', 7.87) + + silicon = openmc.Material() + silicon.add_nuclide('Si28', 1.0) + silicon.set_density('g/cm3', 2.33) + + inner = openmc.Sphere(r=5.0) + outer = openmc.Sphere(r=10.0, boundary_type='vacuum') + iron_cell = openmc.Cell(fill=iron, region=-inner) + silicon_cell = openmc.Cell(fill=silicon, region=+inner & -outer) + + model = openmc.Model() + model.geometry = openmc.Geometry([iron_cell, silicon_cell]) + model.settings.run_mode = 'fixed source' + model.settings.particles = 200 + model.settings.batches = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1.0e6], [1.0])) + + filters = [ + openmc.CellFilter([iron_cell, silicon_cell]), + openmc.EnergyFilter([1.0e6*(1 - 1e-9), 1.0e6*(1 + 1e-9)]), + ] + flux = openmc.Tally() + flux.filters = filters + flux.scores = ['flux'] + flux.estimator = estimator + + total = openmc.Tally() + total.filters = filters + total.nuclides = ['Si28'] + total.scores = ['total'] + total.multiply_density = False + total.estimator = estimator + model.tallies = openmc.Tallies([flux, total]) + + model.run(apply_tally_results=True) + + iron_flux, silicon_flux = flux.mean.ravel() + iron_total, silicon_total = total.mean.ravel() + assert iron_flux > 0.0 + assert silicon_total > 0.0 + assert iron_total / iron_flux == pytest.approx( + silicon_total / silicon_flux, rel=1e-10) diff --git a/tests/unit_tests/test_photon_micro_xs.py b/tests/unit_tests/test_photon_micro_xs.py index bab69043197..19dcf29bcf0 100644 --- a/tests/unit_tests/test_photon_micro_xs.py +++ b/tests/unit_tests/test_photon_micro_xs.py @@ -53,9 +53,9 @@ def test_nuclide_absent_from_the_material(run_in_tmpdir): multiply_density has to be off for the score to say anything: an absent nuclide has an atom density of zero, which multiplies away whatever came out of the cache. With it off the score is the microscopic cross section - itself. The element is still present in the material, so transport keeps - its cache current and this does not exercise the refresh for an element - the scoring region lacks. test_element_absent_from_the_scoring_region does. + itself. The element is still present in the material, so transport keeps + its cache current and this does not exercise the refresh for an element + the scoring region lacks. test_element_absent_from_the_scoring_region does. """ openmc.reset_auto_ids()