Fix collision track feature for photon transport (#3946)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Patrick Shriwise 2026-06-08 18:37:34 -05:00 committed by GitHub
parent 4f6a25e00a
commit ea6ba328c9
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
32 changed files with 167 additions and 196 deletions

View file

@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
(``collision_track.N.h5``) or when the run is performed in parallel. The file
contains the information needed to reconstruct each recorded collision.
The current revision of the collision track file format is 1.1.
The current revision of the collision track file format is 1.2.
**/**
@ -33,7 +33,7 @@ The current revision of the collision track file format is 1.1.
- ``event_mt`` (*int*) -- ENDF MT number identifying the reaction.
- ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events).
- ``cell_id`` (*int*) -- ID of the cell in which the collision occurred.
- ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format).
- ``nuclide_id`` (*int*) -- PDG number of the nuclide (100ZZZAAAM).
- ``material_id`` (*int*) -- ID of the material containing the collision site.
- ``universe_id`` (*int*) -- ID of the universe containing the collision site.
- ``n_collision`` (*int*) -- Collision counter for the particle history.

View file

@ -98,6 +98,11 @@ sub-elements:
A list of strings representing the nuclide, to define specific
define specific target nuclide collisions to be banked.
.. note::
Electron and positron collision-track events are not associated with
a specific nuclide. If a ``nuclides`` entry is specified, these events
are omitted.
*Default*: None
:reactions:

View file

@ -792,6 +792,11 @@ collision_track.h5 file at the end of the simulation. The file contains
300 recorded collisions that occurred in materials with IDs 1 or 2, involving
fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells
with IDs 5 and 12.
.. note::
Electron and positron collision-track events are not associated with a
specific nuclide. If a ``nuclides`` entry is specified, these events are omitted.
The file can be read using :func:`openmc.read_collision_track_file`.
The example below shows how to extract the data from the collision_track
feature and displays the fields stored in the file:

View file

@ -35,7 +35,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
// ============================================================================
// ADJUSTABLE PARAMETERS

View file

@ -84,6 +84,9 @@ public:
double collapse_rate(int MT, double temperature, span<const double> energy,
span<const double> flux) const;
//! Return a ParticleType object representing this nuclide
ParticleType particle_type() const { return {Z_, A_, metastable_}; }
//============================================================================
// Data members
std::string name_; //!< Name of nuclide, e.g. "U235"

View file

@ -200,8 +200,13 @@ void collision_track_record(Particle& particle)
return;
int cell_id = model::cells[cell_index]->id_;
const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
std::string nuclide = nuclide_ptr->name_;
std::string nuclide {};
int nuclide_id = 0;
if (particle.event_nuclide() != NUCLIDE_NONE) {
const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
nuclide = nuclide_ptr->name_;
nuclide_id = nuclide_ptr->particle_type().pdg_number();
}
int universe_id = model::universes[particle.lowest_coord().universe()]->id_;
double delta_E = particle.E_last() - particle.E();
int material_index = particle.material();
@ -224,8 +229,7 @@ void collision_track_record(Particle& particle)
site.event_mt = particle.event_mt();
site.delayed_group = particle.delayed_group();
site.cell_id = cell_id;
site.nuclide_id =
10000 * nuclide_ptr->Z_ + 10 * nuclide_ptr->A_ + nuclide_ptr->metastable_;
site.nuclide_id = nuclide_id;
site.material_id = material_id;
site.universe_id = universe_id;
site.n_collision = particle.n_collision();

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<reactions>(n,fission) 101</reactions>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<cell_ids>22</cell_ids>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<material_ids>1</material_ids>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<nuclides>O16 U235</nuclides>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -52,7 +52,7 @@
<collision_track>
<cell_ids>22</cell_ids>
<universe_ids>77</universe_ids>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -51,7 +51,7 @@
</source>
<collision_track>
<deposited_E_threshold>550000.0</deposited_E_threshold>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -38,9 +38,9 @@
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>80</particles>
<batches>5</batches>
<inactive>1</inactive>
<inactive>4</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
@ -56,7 +56,7 @@
<material_ids>1 11</material_ids>
<nuclides>U238 U235 H1 U234</nuclides>
<deposited_E_threshold>100000.0</deposited_E_threshold>
<max_collisions>300</max_collisions>
<max_collisions>100</max_collisions>
</collision_track>
<seed>1</seed>
</settings>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -1,57 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="11.0" units="g/cm3"/>
<nuclide name="U234" ao="0.0004524"/>
<nuclide name="U235" ao="0.0506068"/>
<nuclide name="U238" ao="0.948709"/>
<nuclide name="U236" ao="0.0002318"/>
<nuclide name="O16" ao="2.0"/>
</material>
<material id="11">
<density value="1.0" units="g/cm3"/>
<nuclide name="H1" ao="2.0"/>
<nuclide name="O16" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
<cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
<cell id="22" material="1" region="-1 2 -3" universe="77"/>
<cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
<surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
<surface id="2" type="z-plane" coeffs="-2.0"/>
<surface id="3" type="z-plane" coeffs="2.0"/>
<surface id="4" type="x-plane" coeffs="-3.0"/>
<surface id="5" type="x-plane" coeffs="3.0"/>
<surface id="6" type="y-plane" coeffs="-3.0"/>
<surface id="7" type="y-plane" coeffs="3.0"/>
<surface id="8" type="z-plane" coeffs="-3.0"/>
<surface id="9" type="z-plane" coeffs="3.0"/>
<surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
<surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
<surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
<surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
<surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
<surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>5</batches>
<inactive>1</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
</space>
<constraints>
<fissionable>true</fissionable>
</constraints>
</source>
<collision_track>
<max_collisions>200</max_collisions>
</collision_track>
<seed>1</seed>
</settings>
</model>

View file

@ -1,2 +0,0 @@
k-combined:
5.642735E-02 1.494035E-02

View file

@ -59,28 +59,13 @@ TODO:
"""
import os
import openmc
import openmc.lib
import pytest
from tests.testing_harness import CollisionTrackTestHarness
from tests.regression_tests import config
@pytest.fixture(scope="function")
def two_threads(monkeypatch):
"""Set the number of OMP threads to 2 for the test."""
monkeypatch.setenv("OMP_NUM_THREADS", "2")
@pytest.fixture(scope="function")
def single_process(monkeypatch):
"""Set the number of MPI process to 1 for the test."""
monkeypatch.setitem(config, "mpi_np", "1")
@pytest.fixture(scope="module")
def model_1():
"""Cylindrical core contained in a first box which is contained in a larger box.
@ -181,9 +166,9 @@ def model_1():
# =============================================================================
model.settings = openmc.Settings()
model.settings.particles = 100
model.settings.particles = 80
model.settings.batches = 5
model.settings.inactive = 1
model.settings.inactive = 4
model.settings.seed = 1
bounds = [
@ -203,19 +188,19 @@ def model_1():
@pytest.mark.parametrize(
"folder, model_name, parameter",
[("case_1_Reactions", "model_1", {"max_collisions": 300, "reactions": ["(n,fission)", 101]}),
[("case_1_Reactions", "model_1", {"max_collisions": 100, "reactions": ["(n,fission)", 101]}),
("case_2_Cell_ID", "model_1", {
"max_collisions": 300, "cell_ids": [22]}),
"max_collisions": 100, "cell_ids": [22]}),
("case_3_Material_ID", "model_1", {
"max_collisions": 300, "material_ids": [1]}),
"max_collisions": 100, "material_ids": [1]}),
("case_4_Nuclide_ID", "model_1", {
"max_collisions": 300, "nuclides": ["O16", "U235"]}),
"max_collisions": 100, "nuclides": ["O16", "U235"]}),
("case_5_Universe_ID", "model_1", {
"max_collisions": 300, "cell_ids": [22], "universe_ids": [77]}),
"max_collisions": 100, "cell_ids": [22], "universe_ids": [77]}),
("case_6_deposited_energy_threshold", "model_1", {
"max_collisions": 300, "deposited_E_threshold": 5.5e5}),
"max_collisions": 100, "deposited_E_threshold": 5.5e5}),
("case_7_all_parameters_used_together", "model_1", {
"max_collisions": 300,
"max_collisions": 100,
"reactions": ["elastic", 18, "(n,disappear)"],
"material_ids": [1, 11],
"universe_ids": [77],
@ -235,21 +220,3 @@ def test_collision_track_several_cases(
"statepoint.5.h5", model=model, workdir=folder
)
harness.main()
@pytest.mark.skipif(config["event"], reason="Results from history-based mode.")
def test_collision_track_2threads(model_1, two_threads, single_process):
# This test checks that the `max_collisions` setting is honored:
# no collisions beyond the specified limit should be recorded.
#
# For the result to be reproducible, the number of threads and
# the transport mode (history vs. event) must remain fixed.
assert os.environ["OMP_NUM_THREADS"] == "2"
assert config["mpi_np"] == "1"
model_1.settings.collision_track = {
"max_collisions": 200
}
harness = CollisionTrackTestHarness(
"statepoint.5.h5", model=model_1, workdir="case_8_2threads"
)
harness.main()

View file

@ -546,19 +546,19 @@ class CollisionTrackTestHarness(PyAPITestHarness):
def _test_output_created(self):
"""Make sure collision_track.h5 has also been created."""
super()._test_output_created()
if self._model.settings.collision_track:
assert os.path.exists(
"collision_track.h5"
), "collision_track file has not been created."
def _compare_output(self):
def _compare_results(self):
"""Compare collision_track.h5 files."""
if self._model.settings.collision_track:
collision_track_true = self._return_collision_track_data(
"collision_track_true.h5")
collision_track_test = self._return_collision_track_data(
"collision_track.h5")
assert collision_track_true.shape == collision_track_test.shape
np.testing.assert_allclose(
collision_track_true, collision_track_test, rtol=1e-07)
@ -582,15 +582,18 @@ class CollisionTrackTestHarness(PyAPITestHarness):
def _overwrite_results(self):
"""Also add the 'collision_track.h5' file during overwriting."""
super()._overwrite_results()
if os.path.exists("collision_track.h5"):
shutil.copyfile("collision_track.h5", "collision_track_true.h5")
def _write_results(self, results_string):
# The result file for this test are written by the OpenMC executable itself
pass
@staticmethod
def _return_collision_track_data(filepath):
"""
Read a collision_track file and return a sorted array composed
of flatten arrays of collision information.
Read a collision_track file and return a sorted array composed of
flattened collision records.
Parameters
----------
@ -600,42 +603,58 @@ class CollisionTrackTestHarness(PyAPITestHarness):
Returns
-------
data : np.array
Sorted array composed of flatten arrays of collision_track data for
each collision information
Sorted array composed of flattened collision-track records.
"""
data = []
keys = []
# Read source file
source = openmc.read_collision_track_file(filepath)
for src in source:
r = src['r']
u = src['u']
e = src['E']
de = src['dE']
time = src['time']
wgt = src['wgt']
delayed_group = src['delayed_group']
cell_id = src['cell_id']
nuclide_id = src['nuclide_id']
material_id = src['material_id']
universe_id = src['universe_id']
n_collision = src['n_collision']
event_mt = src['event_mt']
key = (
f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}"
f"{e:.10e} {de:.10e} {time:.10e} {wgt:.10e} {event_mt} {delayed_group} {cell_id}"
f"{nuclide_id} {material_id} {universe_id} {n_collision} "
)
keys.append(key)
values = [*r, *u, e, de, time, wgt, event_mt,
delayed_group, cell_id, nuclide_id, material_id,
universe_id, n_collision]
assert len(values) == 17
data.append(values)
columns = [
source['r']['x'],
source['r']['y'],
source['r']['z'],
source['u']['x'],
source['u']['y'],
source['u']['z'],
source['E'],
source['dE'],
source['time'],
source['wgt'],
source['event_mt'],
source['delayed_group'],
source['cell_id'],
source['nuclide_id'],
source['material_id'],
source['universe_id'],
source['n_collision'],
source['particle'],
source['parent_id'],
source['progeny_id'],
]
data = np.column_stack(columns)
data = np.array(data)
keys = np.array(keys)
sorted_idx = np.argsort(keys, kind='stable')
# Sort by the complete record, prioritizing stable integer identifiers
# before floating-point fields. This removes dependence on the order in
# which threads append otherwise reproducible collision records.
sort_columns = [
source['parent_id'],
source['progeny_id'],
source['n_collision'],
source['particle'],
source['cell_id'],
source['material_id'],
source['universe_id'],
source['nuclide_id'],
source['event_mt'],
source['delayed_group'],
source['r']['x'],
source['r']['y'],
source['r']['z'],
source['u']['x'],
source['u']['y'],
source['u']['z'],
source['E'],
source['dE'],
source['time'],
source['wgt'],
]
sorted_idx = np.lexsort(tuple(reversed(sort_columns)))
return data[sorted_idx]

View file

@ -34,6 +34,7 @@ def geometry():
{"max_collisions": 200, "mcpl": True}
],
ids=str
)
def test_xml_serialization(parameter, run_in_tmpdir):
"""Check that the different use cases can be written and read in XML."""
@ -45,7 +46,7 @@ def test_xml_serialization(parameter, run_in_tmpdir):
assert read_settings.collision_track == parameter
@pytest.fixture(scope="module")
@pytest.fixture
def model():
"""Simple hydrogen sphere divided in two hemispheres
by a z-plane to form 2 cells."""
@ -127,3 +128,43 @@ def test_format_similarity(run_in_tmpdir, model):
np.testing.assert_allclose(data_h5, data_mcpl, rtol=1e-05)
# tolerance not that low due to the strings that is saved in MCPL,
# not enough precision!
def test_photon_particles(run_in_tmpdir, model):
"""Test that the collision track can be used to track photon particles."""
model.settings.collision_track = {"max_collisions": 200, "cell_ids": [1, 2]}
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Box(*model.geometry.bounding_box),
energy=openmc.stats.delta_function(1e5),
particle='photon'
)
model.run()
with h5py.File("collision_track.h5", "r") as f:
source = f["collision_track_bank"]
assert len(source) < 200
allowed_particles = (openmc.ParticleType.PHOTON, openmc.ParticleType.ELECTRON)
for point in source:
particle_type = openmc.ParticleType(point['particle'])
assert particle_type in allowed_particles
if particle_type == openmc.ParticleType.ELECTRON:
assert point['nuclide_id'] == 0
def test_collision_track_two_threads(model, run_in_tmpdir):
# This test checks that the `max_collisions` setting is honored:
# no collisions beyond the specified limit should be recorded.
#
# The exact set of events in the capped bank is not reproducible with
# multiple threads because the bank stores whichever thread appends first
# until capacity is reached.
model.settings.collision_track = {"max_collisions": 200}
model.run(threads=2, particles=500)
collision_track = openmc.read_collision_track_hdf5("collision_track.h5")
assert len(collision_track) == 200