From 55fd15a8295109b7d451c36181395c5e1dee27ce Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 16:50:32 -0400 Subject: [PATCH 1/4] Mark multigroup source files and keep their groups in FileSource In a multigroup simulation the source bank stores the energy group index of each site in E, and so do the source point files written from it (source..h5 and source.h5). FileSource read E as an energy in eV and converted it to a group: with two groups, a site of group 1 was started in group 0, and a site of group 0 was rejected by the default energy bounds, which require E > 0, or, with wider bounds, given the invalid group index G when the lowest group edge is above zero. These source point files now have an energy_mode attribute equal to "multi-group". FileSource reads it, or the energy_mode dataset of a state point, and then keeps the stored group of each site: it checks that E is an integer from 0 to G - 1, skips the conversion from eV and does not apply energy bounds to these sites, with a warning when bounds are given. A continuous-energy simulation refuses such a file, and an unrecognized energy_mode value is an error. Files without the marker, such as surface source files, which store energies, older source files and files written by openmc.write_source_file, are read in eV as before. The decision is made per source, so the other sources of a simulation are converted as before. The file format version is unchanged: the attribute is optional, and older versions ignore it. --- docs/source/io_formats/source.rst | 10 ++++- docs/source/io_formats/statepoint.rst | 6 ++- docs/source/usersguide/settings.rst | 10 +++++ include/openmc/source.h | 7 ++++ include/openmc/state_point.h | 9 +++- openmc/source.py | 33 +++++++++++---- src/simulation.cpp | 4 +- src/source.cpp | 60 +++++++++++++++++++++++++-- src/state_point.cpp | 9 ++-- 9 files changed, 127 insertions(+), 21 deletions(-) diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 2e81e0a3f2f..7d3aca15ebd 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -17,6 +17,12 @@ following the same format. :Attributes: - **filetype** (*char[]*) -- String indicating the type of file. - **version** (*int[2]*) -- Major and minor version of the source file format. + - **energy_mode** (*char[]*) -- 'multi-group', indicating that + ``E`` holds energy group indices. Only present in source point + files written by a multigroup simulation. Surface source files + and files written by a continuous-energy simulation or by + :func:`openmc.write_source_file` do not have this attribute, and + their ``E`` is an energy in [eV]. :Datasets: @@ -25,4 +31,6 @@ following the same format. ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, which represent the position, direction, energy, time, weight, delayed group, surface ID, and particle type (PDG number), - respectively. + respectively. When the **energy_mode** attribute is + 'multi-group', ``E`` is the energy group index of the particle + (0 for the highest-energy group) rather than its energy in [eV]. diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 6e0561e2481..49cb9f4e910 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -56,8 +56,10 @@ The current version of the statepoint file format is 18.2. ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and ``particle``, which represent the position, direction, energy, time, weight, delayed group, surface ID, and particle type - (PDG number), respectively. Only present when `run_mode` is - 'eigenvalue'. + (PDG number), respectively. When `energy_mode` is 'multi-group', + ``E`` is the energy group index of the particle (0 for the + highest-energy group) rather than its energy in [eV]. Only present + when `run_mode` is 'eigenvalue'. **/tallies/** diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5fad64f3c10..9c797c6e644 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -352,6 +352,16 @@ source file can be manually generated with the :func:`openmc.write_source_file` function. This is particularly useful for coupling OpenMC with another program that generates a source to be used in OpenMC. +In a multigroup simulation, the source point files and state points store the +energy group index of each site rather than its energy, and source point files +record this with an ``energy_mode`` attribute (see :ref:`io_source`). When a +multigroup simulation reads such a file, each particle starts in its stored +group; energy bounds given as source constraints are not applied to these sites, +and a warning says so. A continuous-energy simulation cannot use such a file. +Files without this record, such as surface source files and files written by +:func:`openmc.write_source_file`, are read as energies in [eV], which a +multigroup simulation converts to groups. + Surface Sources +++++++++++++++ diff --git a/include/openmc/source.h b/include/openmc/source.h index 63e2e244f93..d6710bd1b01 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -94,6 +94,10 @@ class Source { bool satisfies_spatial_constraints(Position r) const; + //! Whether the energy of a sampled site is a multigroup energy group index + //! rather than an energy in [eV] + virtual bool energy_is_group_index() const { return false; } + protected: // Strategy used for rejecting sites when constraints are applied. KILL means // that sites are always accepted but if they don't satisfy constraints, they @@ -183,11 +187,14 @@ class FileSource : public Source { void load_sites_from_file( const std::string& path); //!< Load source sites from file + bool energy_is_group_index() const override { return group_index_; } + protected: SourceSite sample(uint64_t* seed) const override; private: vector sites_; //!< Source sites + bool group_index_ {false}; //!< Whether site energies are group indices }; //============================================================================== diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index fb1aaf7b985..4646a701599 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -36,11 +36,16 @@ void load_state_point(); // // The source_bank variable is used as work space if MPI is used, // so it cannot be given as a const span. +// +// If group_indices is true, the energies of the sites are multigroup energy +// group indices, which the file records with an "energy_mode" attribute equal +// to "multi-group" so that a FileSource reads them as group indices. MCPL +// files do not record it. void write_h5_source_point(const char* filename, span source_bank, - const vector& bank_index); + const vector& bank_index, bool group_indices = false); void write_source_point(std::string, span source_bank, - const vector& bank_index, bool use_mcpl); + const vector& bank_index, bool use_mcpl, bool group_indices = false); // This appends a source bank specification to an HDF5 file // that's already open. It is used internally by write_source_point. diff --git a/openmc/source.py b/openmc/source.py index e68bdd5cf25..61ba2e1c3a2 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -816,11 +816,13 @@ class FileSource(SourceBase): :class:`openmc.Universe` for which sampled sites must be within. For 'time_bounds' and 'energy_bounds', the corresponding value is a sequence of floats giving the lower and upper bounds on time in [s] or energy in - [eV] that the sampled particle must be within. For 'fissionable', the - value is a bool indicating that only sites in fissionable material - should be accepted. The 'rejection_strategy' indicates what should - happen when a source particle is rejected: either 'resample' (pick a new - particle) or 'kill' (accept and terminate). + [eV] that the sampled particle must be within. Energy bounds are not + applied to the sites of a file that holds multigroup energy group + indices (see below). For 'fissionable', the value is a bool indicating + that only sites in fissionable material should be accepted. The + 'rejection_strategy' indicates what should happen when a source + particle is rejected: either 'resample' (pick a new particle) or 'kill' + (accept and terminate). Attributes ---------- @@ -835,6 +837,15 @@ class FileSource(SourceBase): 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + Notes + ----- + Source point files and state points written by a multigroup simulation + store the energy group index of each site rather than its energy (see + :ref:`io_source`). A multigroup simulation starts each particle from such a + file in its stored group, and a continuous-energy simulation cannot use + such a file. Other files, such as surface source files and files written by + :func:`openmc.write_source_file`, are read as energies in [eV]. + """ def __init__( @@ -1337,7 +1348,8 @@ class SourceParticle: u : iterable of float Directional cosines E : float - Energy of particle in [eV] + Energy of particle in [eV]. For a particle read from a source file + written by a multigroup simulation, the energy group index. time : float Time of particle in [s] wgt : float @@ -1435,6 +1447,10 @@ class ParticleList(list): def from_hdf5(cls, filename: PathLike) -> ParticleList: """Create particle list from an HDF5 file. + Energies are returned as stored. For a source file written by a + multigroup simulation (``energy_mode`` attribute 'multi-group'), they + are energy group indices. + Parameters ---------- filename : path-like @@ -1548,7 +1564,10 @@ def export_to_hdf5(self, filename: PathLike, **kwargs): """Export particle list to an HDF5 file. This method write out an .h5 file that can be used as a source file in - conjunction with the :class:`openmc.FileSource` class. + conjunction with the :class:`openmc.FileSource` class. The file has no + ``energy_mode`` attribute, so OpenMC reads its energies in [eV], also + when the particles were read from a source file written by a multigroup + simulation. Parameters ---------- diff --git a/src/simulation.cpp b/src/simulation.cpp index 4cef7a8e628..ecf8e4b349e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -573,7 +573,7 @@ void finalize_batch() settings::path_output, simulation::current_batch, w); span bankspan(simulation::source_bank); write_source_point(source_point_filename, bankspan, - simulation::work_index, settings::source_mcpl_write); + simulation::work_index, settings::source_mcpl_write, !settings::run_CE); } // Write a continously-overwritten source point if requested. @@ -581,7 +581,7 @@ void finalize_batch() auto filename = settings::path_output + "source"; span bankspan(simulation::source_bank); write_source_point(filename, bankspan, simulation::work_index, - settings::source_mcpl_write); + settings::source_mcpl_write, !settings::run_CE); } } diff --git a/src/source.cpp b/src/source.cpp index f951722cc18..ff26e3457b2 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -214,10 +214,12 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const if (constraints_applied()) { accepted = true; } else { - // Check whether sampled site satisfies constraints - accepted = satisfies_spatial_constraints(site.r) && - satisfies_energy_constraints(site.E) && - satisfies_time_constraints(site.time); + // Check whether sampled site satisfies constraints. Energy bounds are + // not applied to multigroup energy group indices. + accepted = + satisfies_spatial_constraints(site.r) && + (energy_is_group_index() || satisfies_energy_constraints(site.E)) && + satisfies_time_constraints(site.time); if (!accepted) { ++n_local_reject; @@ -524,6 +526,44 @@ void FileSource::load_sites_from_file(const std::string& path) fatal_error("Specified starting source file not a source file type."); } + // Source files written by a multigroup simulation store an energy group + // index rather than an energy for each site, which an "energy_mode" + // attribute equal to "multi-group" indicates. A state point records the + // energy mode of its simulation as a dataset. Without either, the sites + // hold energies in [eV]. + std::string energy_mode {"continuous-energy"}; + if (filetype == "source") { + if (attribute_exists(file_id, "energy_mode")) { + read_attribute(file_id, "energy_mode", energy_mode); + } + } else if (object_exists(file_id, "energy_mode")) { + read_dataset(file_id, "energy_mode", energy_mode); + } + if (energy_mode == "multi-group") { + group_index_ = true; + } else if (energy_mode != "continuous-energy") { + fatal_error( + fmt::format("Source file '{}' has an unrecognized energy mode '{}'.", + path, energy_mode)); + } + + if (group_index_) { + if (settings::run_CE) { + fatal_error(fmt::format( + "Source file '{}' holds multigroup energy group indices and cannot " + "be used in a continuous-energy simulation.", + path)); + } + if (mpi::master && + (energy_bounds_.first != 0.0 || + energy_bounds_.second != std::numeric_limits::max())) { + warning(fmt::format( + "Energy bounds are not applied to the sites of source file '{}', " + "which holds multigroup energy group indices.", + path)); + } + } + // Read in the source particles read_source_bank(file_id, sites_, false); @@ -1222,6 +1262,18 @@ SourceSite sample_external_source(uint64_t* seed) model::external_sources[i]->strength() * n_sources / total_strength; } + // Sites from a source file written by a multigroup simulation already hold + // an energy group index, which only needs to be checked + if (model::external_sources[i]->energy_is_group_index()) { + int n_groups = data::mg.num_energy_groups_; + if (!(site.E >= 0.0 && site.E < n_groups && site.E == std::floor(site.E))) { + fatal_error(fmt::format("Source site has energy group index {}, which is " + "not valid for {} energy groups.", + site.E, n_groups)); + } + return site; + } + // If running in MG, convert site.E to group if (!settings::run_CE) { site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(), diff --git a/src/state_point.cpp b/src/state_point.cpp index cbfb13f7670..5131faaeb81 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -596,7 +596,7 @@ hid_t h5banktype(bool memory) } void write_source_point(std::string filename, span source_bank, - const vector& bank_index, bool use_mcpl) + const vector& bank_index, bool use_mcpl, bool group_indices) { std::string ext = use_mcpl ? "mcpl" : "h5"; @@ -616,12 +616,13 @@ void write_source_point(std::string filename, span source_bank, write_mcpl_source_point(filename.c_str(), source_bank, bank_index); } else { filename.append(".h5"); - write_h5_source_point(filename.c_str(), source_bank, bank_index); + write_h5_source_point( + filename.c_str(), source_bank, bank_index, group_indices); } } void write_h5_source_point(const char* filename, span source_bank, - const vector& bank_index) + const vector& bank_index, bool group_indices) { // When using parallel HDF5, the file is written to collectively by all // processes. With MPI-only, the file is opened and written by the master @@ -648,6 +649,8 @@ void write_h5_source_point(const char* filename, span source_bank, file_id = file_open(filename_.c_str(), 'w', true); write_attribute(file_id, "filetype", "source"); write_attribute(file_id, "version", VERSION_STATEPOINT); + if (group_indices) + write_attribute(file_id, "energy_mode", "multi-group"); } // Get pointer to source bank and write to file From 0590cf0b8d960f8a70fa36eb084524a3d283dfb0 Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 19:08:39 -0400 Subject: [PATCH 2/4] Test multigroup source files read by FileSource Add tests of the energy_mode marker that multigroup simulations write to their source point files and of how a FileSource reads them. The marker is on source..h5 and source.h5 only. Fixed source and eigenvalue simulations from such a file start every particle in its stored group, and sampling keeps that group with default and with widened energy bounds; one warning is printed when bounds are given, while domain and time constraints still apply. A multigroup state point is read through its energy_mode dataset. Group indices that are out of range, fractional or NaN, an unknown energy_mode value, and a marked file or multigroup state point in a continuous-energy simulation are refused. Files that hold energies are unmarked and read in eV as before: surface source files written by a multigroup simulation, source points of a continuous-energy simulation and files written by openmc.write_source_file, also when another source of the simulation holds group indices. With --mpi, the writer and two readers run on several processes. --- tests/unit_tests/test_source_file_mg.py | 563 ++++++++++++++++++++++++ 1 file changed, 563 insertions(+) create mode 100644 tests/unit_tests/test_source_file_mg.py diff --git a/tests/unit_tests/test_source_file_mg.py b/tests/unit_tests/test_source_file_mg.py new file mode 100644 index 00000000000..530e27e50c1 --- /dev/null +++ b/tests/unit_tests/test_source_file_mg.py @@ -0,0 +1,563 @@ +"""Tests of source files written and read by multigroup simulations. + +In a multigroup simulation the source bank stores the energy group index of +each site in E. The source point files written from it are marked with an +energy_mode attribute equal to "multi-group", and a FileSource of such a file, +or of a multigroup state point, keeps the stored group of each site instead of +converting E from eV. Files without the marker, such as surface source files +and files written by openmc.write_source_file, hold energies in eV. + +The tests that run OpenMC in a subprocess come first; the tests that sample +the source through openmc.lib follow. +""" +import math +from pathlib import Path +import re +import shutil +import sys + +import h5py +import numpy as np +import openmc +import openmc.lib +import pytest + +from tests.regression_tests import config + +# Group edges in [eV] and fission spectrum of each library. Group 0 is the +# highest-energy group. The medium is fissile only in the group that the +# fission spectrum feeds, so every fission site is created in that group. +LIBRARIES = { + '1g': ([1.0e-5, 2.0e7], [1.0]), + '1g-zero-edge': ([0.0, 2.0e7], [1.0]), + '2g-thermal': ([1.0e-5, 0.625, 2.0e7], [0.0, 1.0]), + '2g-fast': ([1.0e-5, 0.625, 2.0e7], [1.0, 0.0]), + '2g-fast-zero-edge': ([0.0, 0.625, 2.0e7], [1.0, 0.0]), +} + +# Energy bounds wider than every group structure above +WIDE_BOUNDS = (-1.0, 3.0e7) + +# Average energies of the two groups of the two-group libraries +LOW_AVERAGE = 0.5*(1.0e-5 + 0.625) +HIGH_AVERAGE = 0.5*(0.625 + 2.0e7) + + +def birth_group(name): + """Energy group in which the fission neutrons of a library are born""" + return LIBRARIES[name][1].index(1.0) + + +def make_model(name, directory='.', plane=None): + """Reflective 10 cm cube of one material, with its library + + The cube is divided into an inner cell [-2, 2]^3 and an outer cell or, if + a plane is given, into the two half-spaces of that plane. Each group has a + total cross section of 1 /cm, half of it absorption and half scattering + within the group, so particles never change group. The library is written + to mgxs.h5 in the given directory. Returns the model and its two cells. + + """ + edges, chi = LIBRARIES[name] + groups = openmc.mgxs.EnergyGroups(edges) + n = groups.num_groups + fission = 0.2*np.array(chi) + xsdata = openmc.XSdata('medium', groups) + xsdata.order = 0 + xsdata.set_total(np.ones(n)) + xsdata.set_absorption(np.full(n, 0.5)) + xsdata.set_scatter_matrix(np.diag(np.full(n, 0.5))[..., np.newaxis]) + xsdata.set_fission(fission) + xsdata.set_nu_fission(2.0*fission) + xsdata.set_chi(chi) + library = openmc.MGXSLibrary(groups) + library.add_xsdata(xsdata) + library.export_to_hdf5(Path(directory) / 'mgxs.h5') + + mat = openmc.Material() + mat.set_density('macro', 1.0) + mat.add_macroscopic('medium') + cube = openmc.model.RectangularParallelepiped( + -5.0, 5.0, -5.0, 5.0, -5.0, 5.0, boundary_type='reflective') + if plane is None: + inner = openmc.model.RectangularParallelepiped( + -2.0, 2.0, -2.0, 2.0, -2.0, 2.0) + cells = [openmc.Cell(fill=mat, region=-inner), + openmc.Cell(fill=mat, region=+inner & -cube)] + else: + cells = [openmc.Cell(fill=mat, region=-plane & -cube), + openmc.Cell(fill=mat, region=+plane & -cube)] + + model = openmc.Model() + model.materials = openmc.Materials([mat]) + model.materials.cross_sections = 'mgxs.h5' + model.geometry = openmc.Geometry(cells) + model.settings.energy_mode = 'multi-group' + model.settings.particles = 1000 + model.settings.batches = 3 + # Histories here take tens of events; the limit stops one that would run + # forever, such as a particle in an invalid group if a test fails + model.settings.max_particle_events = 1000 + return model, cells + + +def box_source(energies, probabilities): + return openmc.IndependentSource( + space=openmc.stats.Box((-5.0, -5.0, -5.0), (5.0, 5.0, 5.0)), + energy=openmc.stats.Discrete(energies, probabilities)) + + +def add_flux_tally(model, name): + tally = openmc.Tally() + tally.filters = [openmc.EnergyFilter(LIBRARIES[name][0])] + tally.scores = ['flux'] + model.tallies = [tally] + return tally + + +def run(model, cwd='.', mpi=False, **kwargs): + """Run OpenMC in a subprocess, with MPI if requested and enabled""" + if mpi and config['mpi']: + kwargs['mpi_args'] = [config['mpiexec'], '-n', config['mpi_np']] + kwargs.setdefault('output', False) + return model.run(cwd=cwd, openmc_exec=config['exe'], + event_based=config['event'], **kwargs) + + +def run_flux(model, name, **kwargs): + """Run the model and return its flux in each energy bin, lowest first""" + tally = add_flux_tally(model, name) + sp_path = run(model, **kwargs) + with openmc.StatePoint(sp_path) as sp: + return sp.tallies[tally.id].mean.ravel() + + +def write_source_points(directory, name, statepoint=False, mpi=False): + """Run a multigroup eigenvalue simulation that writes source points + + By default the sites of the last batch are written to source.3.h5 and + source.h5; with statepoint=True they are kept in statepoint.3.h5 instead. + + """ + model, _ = make_model(name, directory) + model.settings.run_mode = 'eigenvalue' + model.settings.inactive = 1 + # Start the particles in the group in which the medium is fissile + model.settings.source = box_source( + [0.1 if birth_group(name) == 1 else 1.0e6], [1.0]) + if statepoint: + model.settings.sourcepoint = {'batches': [3]} + else: + model.settings.sourcepoint = { + 'batches': [3], 'separate': True, 'overwrite': True} + run(model, cwd=directory, mpi=mpi) + + +@pytest.fixture(scope='module') +def mg_files(tmp_path_factory): + """Directory with the source points of a library, written once""" + directories = {} + + def get(name, statepoint=False): + if (name, statepoint) not in directories: + directory = tmp_path_factory.mktemp('mg_source') + write_source_points(directory, name, statepoint) + directories[name, statepoint] = directory + return directories[name, statepoint] + + return get + + +@pytest.fixture(scope='module') +def surface_source_file(tmp_path_factory): + """Surface source file written by a two-group fixed source simulation""" + directory = tmp_path_factory.mktemp('mg_surface_source') + plane = openmc.XPlane(0.0) + model, _ = make_model('2g-thermal', directory, plane=plane) + model.settings.run_mode = 'fixed source' + model.settings.source = box_source([0.1, 1.0e6], [0.5, 0.5]) + model.settings.surf_source_write = { + 'surface_ids': [plane.id], 'max_particles': 100_000} + run(model, cwd=directory) + return directory / 'surface_source.h5' + + +def copy_source(src, dst, energy_mode=None, E=None, time=None): + """Copy a source file, optionally changing its marker or its sites""" + shutil.copy(src, dst) + with h5py.File(dst, 'r+') as fh: + if energy_mode is not None: + fh.attrs['energy_mode'] = np.bytes_(energy_mode) + if E is not None or time is not None: + bank = fh['source_bank'][()] + if E is not None: + bank['E'] = E + if time is not None: + bank['time'] = time + fh['source_bank'][...] = bank + + +def read_bank(path): + with h5py.File(path, 'r') as fh: + return fh['source_bank'][()] + + +def site_keys(sites): + """Position, direction and weight of each site, as tuples""" + if sites.dtype['r'].names: + r = np.column_stack([sites['r'][c] for c in 'xyz']) + u = np.column_stack([sites['u'][c] for c in 'xyz']) + else: + r, u = sites['r'], sites['u'] + return [tuple(x) + tuple(y) + (w,) for x, y, w in + zip(r.tolist(), u.tolist(), sites['wgt'].tolist())] + + +def file_energies(path): + """Map the position, direction and weight of each file site to its E""" + bank = read_bank(path) + return dict(zip(site_keys(bank), bank['E'].tolist())) + + +def write_ev_file(path): + """Write 500 sites at 0.1 eV and 500 at 1 MeV, uniform in the cube""" + rng = np.random.default_rng(1) + positions = rng.uniform(-5.0, 5.0, (1000, 3)) + particles = [ + openmc.SourceParticle(r=tuple(r), u=(0.0, 0.0, 1.0), + E=0.1 if i < 500 else 1.0e6) + for i, r in enumerate(positions) + ] + openmc.write_source_file(particles, path) + + +def count_bounds_warnings(text): + return len(re.findall(r'WARNING: Energy bounds', ' '.join(text.split()))) + + +# --------------------------------------------------------------------------- +# Tests that run OpenMC in a subprocess + + +@pytest.mark.parametrize('name', LIBRARIES) +def test_source_points_marked(name, mg_files): + directory = mg_files(name) + for filename in ('source.3.h5', 'source.h5'): + with h5py.File(directory / filename, 'r') as fh: + assert set(fh.attrs) == {'filetype', 'version', 'energy_mode'} + assert fh.attrs['filetype'] == b'source' + assert fh.attrs['energy_mode'] == b'multi-group' + bank = fh['source_bank'][()] + assert bank.size == 1000 + assert np.all(bank['E'] == birth_group(name)) + assert np.array_equal(read_bank(directory / 'source.3.h5'), + read_bank(directory / 'source.h5')) + + +@pytest.mark.parametrize('name', LIBRARIES) +def test_fixed_source_from_file(name, mg_files, run_in_tmpdir): + shutil.copy(mg_files(name) / 'source.3.h5', 'source.3.h5') + model, _ = make_model(name) + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource('source.3.h5') + flux = run_flux(model, name) + + # The flux lies only in the energy bin of the stored group (bins are in + # increasing energy, groups in decreasing energy) + stored_bin = flux.size - 1 - birth_group(name) + assert flux[stored_bin] > 0.0 + assert np.all(np.delete(flux, stored_bin) == 0.0) + + +def test_eigenvalue_from_file(mg_files, run_in_tmpdir): + # The file holds sites in the low group, the only fissile one + shutil.copy(mg_files('2g-thermal') / 'source.3.h5', 'source.3.h5') + model, _ = make_model('2g-thermal') + model.settings.run_mode = 'eigenvalue' + model.settings.inactive = 1 + model.settings.source = openmc.FileSource('source.3.h5') + flux = run_flux(model, '2g-thermal') + assert flux[0] > 0.0 + assert flux[1] == 0.0 + + +@pytest.mark.parametrize('energy_mode', [None, 'continuous-energy']) +def test_ev_file_transport(energy_mode, run_in_tmpdir): + write_ev_file('ev.h5') + with h5py.File('ev.h5', 'r') as fh: + assert set(fh.attrs) == {'filetype', 'version'} + copy_source('ev.h5', 'source.h5', energy_mode=energy_mode) + model, _ = make_model('2g-thermal') + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource('source.h5') + # Sites at 0.1 eV start in the low group and sites at 1 MeV in the high + # group + assert np.all(run_flux(model, '2g-thermal') > 0.0) + + +def test_unknown_energy_mode(mg_files, run_in_tmpdir): + copy_source(mg_files('2g-thermal') / 'source.3.h5', 'source.h5', + energy_mode='foo') + model, _ = make_model('2g-thermal') + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource('source.h5') + with pytest.raises(RuntimeError, match="'foo'"): + run(model) + + +@pytest.mark.parametrize('E, pattern', [ + (2.0, r'group index 2(\.0)?\b'), + (0.5, r'group index 0\.5\b'), + (-1.0, r'group index -1(\.0)?\b'), + (math.nan, r'group index nan\b'), +], ids=['G', 'fraction', 'negative', 'nan']) +def test_invalid_group_index(E, pattern, mg_files, run_in_tmpdir): + copy_source(mg_files('2g-thermal') / 'source.3.h5', 'source.h5', E=E) + model, _ = make_model('2g-thermal') + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource('source.h5') + with pytest.raises(RuntimeError) as excinfo: + run(model) + message = str(excinfo.value) + assert re.search(pattern, message, re.IGNORECASE) + assert re.search(r'\b2 energy groups', message) + + +@pytest.mark.parametrize('filename', ['source.3.h5', 'statepoint.3.h5']) +def test_continuous_energy_refuses(filename, mg_files, run_in_tmpdir): + statepoint = filename.startswith('statepoint') + shutil.copy(mg_files('2g-thermal', statepoint) / filename, filename) + + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + mat.set_density('g/cm3', 0.1) + cube = openmc.model.RectangularParallelepiped( + -5.0, 5.0, -5.0, 5.0, -5.0, 5.0, boundary_type='vacuum') + model = openmc.Model() + model.materials = openmc.Materials([mat]) + model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-cube)]) + model.settings.run_mode = 'fixed source' + model.settings.particles = 100 + model.settings.batches = 2 + model.settings.source = openmc.FileSource(filename) + with pytest.raises(RuntimeError, + match='multigroup energy group indices'): + run(model) + + +def test_continuous_energy_source_points_unmarked(run_in_tmpdir): + fuel = openmc.Material() + fuel.add_nuclide('U235', 1.0) + fuel.add_nuclide('H1', 100.0) + fuel.set_density('g/cm3', 1.0) + sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + model = openmc.Model() + model.materials = openmc.Materials([fuel]) + model.geometry = openmc.Geometry([openmc.Cell(fill=fuel, region=-sphere)]) + model.settings.particles = 500 + model.settings.batches = 3 + model.settings.inactive = 1 + model.settings.sourcepoint = { + 'batches': [3], 'separate': True, 'overwrite': True} + run(model) + for filename in ('source.3.h5', 'source.h5'): + with h5py.File(filename, 'r') as fh: + assert set(fh.attrs) == {'filetype', 'version'} + + +def test_surface_source_file_unmarked(surface_source_file): + with h5py.File(surface_source_file, 'r') as fh: + assert set(fh.attrs) == {'filetype', 'version'} + E = fh['source_bank']['E'] + # Surface source sites hold the average energy of their group in [eV] + assert set(np.unique(E).tolist()) == {LOW_AVERAGE, HIGH_AVERAGE} + + +def test_surface_source_transport(surface_source_file, run_in_tmpdir): + shutil.copy(surface_source_file, 'surface_source.h5') + model, _ = make_model('2g-thermal', plane=openmc.XPlane(0.0)) + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource('surface_source.h5') + assert np.all(run_flux(model, '2g-thermal') > 0.0) + + +def test_source_points_with_mpi(run_in_tmpdir, capsys): + # With --mpi, the writer and the readers run on several processes + Path('writer').mkdir() + write_source_points('writer', '2g-thermal', mpi=True) + with h5py.File('writer/source.3.h5', 'r') as fh: + assert fh.attrs['energy_mode'] == b'multi-group' + bank = fh['source_bank'][()] + assert bank.size == 1000 + assert np.all(bank['E'] == 1.0) + + for directory in ('reader', 'wide'): + Path(directory).mkdir() + shutil.copy('writer/source.3.h5', directory) + + model, _ = make_model('2g-thermal', 'reader') + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource('source.3.h5') + flux = run_flux(model, '2g-thermal', cwd='reader', mpi=True) + assert flux[0] > 0.0 + assert flux[1] == 0.0 + + # Energy bounds are not applied to group indices, which one warning says + model, _ = make_model('2g-thermal', 'wide') + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource( + 'source.3.h5', constraints={'energy_bounds': WIDE_BOUNDS}) + capsys.readouterr() + run(model, cwd='wide', mpi=True, output=True) + output = capsys.readouterr().out + assert count_bounds_warnings(output) == 1 + if config['mpi']: + assert re.search(rf'MPI Processes \| {config["mpi_np"]}\b', output) + + +# --------------------------------------------------------------------------- +# Tests that sample the source through openmc.lib + + +def sample(model, n_samples=1000): + model.export_to_model_xml() + with openmc.lib.run_in_memory(output=False): + return openmc.lib.sample_external_source( + n_samples, prn_seed=1, as_array=True) + + +@pytest.mark.parametrize('bounds', [None, WIDE_BOUNDS], + ids=['default-bounds', 'wide-bounds']) +@pytest.mark.parametrize('name', LIBRARIES) +def test_sample_file_source(name, bounds, mg_files, run_in_tmpdir): + shutil.copy(mg_files(name) / 'source.3.h5', 'source.3.h5') + model, _ = make_model(name) + constraints = {'energy_bounds': bounds} if bounds else None + model.settings.source = openmc.FileSource( + 'source.3.h5', constraints=constraints) + sites = sample(model) + + # Every sampled site is a file site and keeps its stored group + stored = file_energies('source.3.h5') + assert all(key in stored for key in site_keys(sites)) + assert np.all(sites['E'] == birth_group(name)) + + +@pytest.mark.parametrize('bounds, n_warnings', [ + (None, 0), + ((0.0, sys.float_info.max), 0), + (WIDE_BOUNDS, 1), + ((0.0, 1.0e7), 1), + ((-1.0, sys.float_info.max), 1), +], ids=['none', 'default', 'wide', 'upper', 'lower']) +def test_energy_bounds_warning(bounds, n_warnings, mg_files, run_in_tmpdir, + capfd): + shutil.copy(mg_files('2g-thermal') / 'source.3.h5', 'source.3.h5') + model, _ = make_model('2g-thermal') + constraints = {'energy_bounds': bounds} if bounds else None + model.settings.source = openmc.FileSource( + 'source.3.h5', constraints=constraints) + model.export_to_model_xml() + capfd.readouterr() + with openmc.lib.run_in_memory(output=False): + pass + assert count_bounds_warnings(capfd.readouterr().err) == n_warnings + + +def test_domain_constraint(mg_files, run_in_tmpdir): + shutil.copy(mg_files('2g-thermal') / 'source.3.h5', 'source.3.h5') + model, (inner, _) = make_model('2g-thermal') + model.settings.source_rejection_fraction = 0.001 + model.settings.source = openmc.FileSource( + 'source.3.h5', constraints={'domains': [inner]}) + sites = sample(model) + assert np.all(np.abs(sites['r']) < 2.0) + assert np.all(sites['E'] == 1.0) + + +def test_time_constraint(mg_files, run_in_tmpdir): + # Time bounds still apply to sites that hold group indices + time = np.where(np.arange(1000) % 2 == 0, 1.0, 3.0) + copy_source(mg_files('2g-thermal') / 'source.3.h5', 'source.h5', + time=time) + model, _ = make_model('2g-thermal') + model.settings.source = openmc.FileSource( + 'source.h5', constraints={'time_bounds': [0.0, 2.0]}) + sites = sample(model) + assert np.all(sites['time'] == 1.0) + assert np.all(sites['E'] == 1.0) + + +def test_statepoint_source(mg_files, run_in_tmpdir): + path = mg_files('2g-thermal', statepoint=True) / 'statepoint.3.h5' + shutil.copy(path, 'statepoint.3.h5') + with h5py.File('statepoint.3.h5', 'r') as fh: + assert fh['energy_mode'][()] == b'multi-group' + assert np.all(fh['source_bank']['E'] == 1.0) + model, _ = make_model('2g-thermal') + model.settings.source = openmc.FileSource('statepoint.3.h5') + sites = sample(model) + assert np.all(sites['E'] == 1.0) + + +@pytest.mark.parametrize('energy_mode', [None, 'continuous-energy']) +def test_ev_file_sampling(energy_mode, run_in_tmpdir): + write_ev_file('ev.h5') + copy_source('ev.h5', 'source.h5', energy_mode=energy_mode) + model, _ = make_model('2g-thermal') + model.settings.source = openmc.FileSource('source.h5') + sites = sample(model) + + # The energies in eV are converted to groups: 0.1 eV is in group 1 + stored = file_energies('source.h5') + expected = [1.0 if stored[key] == 0.1 else 0.0 for key in site_keys(sites)] + assert np.array_equal(sites['E'], expected) + assert 0 < np.count_nonzero(sites['E']) < sites.size + + +def test_ev_file_energy_bounds(run_in_tmpdir, capfd): + # Energy bounds still apply to files that hold energies, with no warning + write_ev_file('source.h5') + model, _ = make_model('2g-thermal') + model.settings.source = openmc.FileSource( + 'source.h5', constraints={'energy_bounds': [1.0, 3.0e7]}) + capfd.readouterr() + sites = sample(model) + assert np.all(sites['E'] == 0.0) + assert count_bounds_warnings(capfd.readouterr().err) == 0 + + +def test_surface_source_sampling(surface_source_file, run_in_tmpdir): + shutil.copy(surface_source_file, 'surface_source.h5') + model, _ = make_model('2g-thermal', plane=openmc.XPlane(0.0)) + model.settings.source = openmc.FileSource('surface_source.h5') + sites = sample(model) + + # The group averages in eV map back to their groups + stored = file_energies('surface_source.h5') + expected = [1.0 if stored[key] == LOW_AVERAGE else 0.0 + for key in site_keys(sites)] + assert np.array_equal(sites['E'], expected) + assert 0 < np.count_nonzero(sites['E']) < sites.size + + +def test_mixed_sources(mg_files, run_in_tmpdir): + # A file of group indices does not change how another source is read + shutil.copy(mg_files('2g-thermal') / 'source.3.h5', 'source.3.h5') + model, _ = make_model('2g-thermal') + point = (4.5, 4.5, 4.5) + model.settings.source = [ + openmc.FileSource('source.3.h5'), + openmc.IndependentSource( + space=openmc.stats.Point(point), + energy=openmc.stats.Discrete([1.0e6], [1.0])), + ] + sites = sample(model) + + from_point = np.all(sites['r'] == point, axis=1) + assert 0 < np.count_nonzero(from_point) < sites.size + assert np.all(sites['E'][from_point] == 0.0) + stored = file_energies('source.3.h5') + file_sites = sites[~from_point] + assert all(key in stored for key in site_keys(file_sites)) + assert np.all(file_sites['E'] == 1.0) From 976bb6a000892623583e941bd72aaac959b19761 Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 20:50:56 -0400 Subject: [PATCH 3/4] Document the energy_mode values of source files State that FileSource also accepts 'continuous-energy' (energies in eV, as without the attribute), that any other value is an error, and that files written by ParticleList.export_to_hdf5 have no energy_mode attribute either. --- docs/source/io_formats/source.rst | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 7d3aca15ebd..6944ecae8c6 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -17,12 +17,18 @@ following the same format. :Attributes: - **filetype** (*char[]*) -- String indicating the type of file. - **version** (*int[2]*) -- Major and minor version of the source file format. - - **energy_mode** (*char[]*) -- 'multi-group', indicating that - ``E`` holds energy group indices. Only present in source point - files written by a multigroup simulation. Surface source files - and files written by a continuous-energy simulation or by - :func:`openmc.write_source_file` do not have this attribute, and - their ``E`` is an energy in [eV]. + - **energy_mode** (*char[]*) -- Whether ``E`` holds energy group + indices ('multi-group') or energies in [eV] + ('continuous-energy'). A file without this attribute is read as + 'continuous-energy'. When the file is read by + :class:`openmc.FileSource`, any other value is an error, and so + is a 'multi-group' file in a continuous-energy simulation. + OpenMC writes this attribute, with the value 'multi-group', only + to source point files written by a multigroup simulation. + Surface source files and files written by a continuous-energy + simulation, by :func:`openmc.write_source_file` or by + :meth:`openmc.ParticleList.export_to_hdf5` do not have this + attribute, so their ``E`` is read as an energy in [eV]. :Datasets: From 93d7766f7664871cec13fa222ea305ed56b0b19c Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 20:53:21 -0400 Subject: [PATCH 4/4] Test that energy bounds are not applied to group indices Run a fixed source simulation from the source file of a two-group simulation, whose sites are all in group 1, with energy bounds of 1 keV to 20 MeV. Read as energies in eV, the stored indices would be rejected; as group indices, the sites start in their group, and one warning says that the bounds are not applied. --- tests/unit_tests/test_source_file_mg.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/unit_tests/test_source_file_mg.py b/tests/unit_tests/test_source_file_mg.py index 530e27e50c1..f6c1b7716e0 100644 --- a/tests/unit_tests/test_source_file_mg.py +++ b/tests/unit_tests/test_source_file_mg.py @@ -269,6 +269,21 @@ def test_fixed_source_from_file(name, mg_files, run_in_tmpdir): assert np.all(np.delete(flux, stored_bin) == 0.0) +def test_energy_bounds_not_applied(mg_files, run_in_tmpdir, capsys): + # Read as an energy in eV, the stored group index 1 would lie outside + # these bounds, which are not applied to group indices + shutil.copy(mg_files('2g-thermal') / 'source.3.h5', 'source.3.h5') + model, _ = make_model('2g-thermal') + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.FileSource( + 'source.3.h5', constraints={'energy_bounds': (1.0e3, 2.0e7)}) + capsys.readouterr() + flux = run_flux(model, '2g-thermal', output=True) + assert flux[0] > 0.0 + assert flux[1] == 0.0 + assert count_bounds_warnings(capsys.readouterr().out) == 1 + + def test_eigenvalue_from_file(mg_files, run_in_tmpdir): # The file holds sites in the low group, the only fissile one shutil.copy(mg_files('2g-thermal') / 'source.3.h5', 'source.3.h5')