in memory cross-section alteration#3979
Open
dave-soliton wants to merge 3 commits into
Open
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
This PR adds support for accessing and modifying continuous-energy nuclide reaction cross sections directly from the Python API through the shared library interface.
Currently, perturbation studies involving nuclear data require generating modified HDF5 cross section libraries and reloading them into OpenMC. This workflow introduces additional file generation and I/O overhead and makes iterative sensitivity or uncertainty studies more cumbersome.
The additions in this PR expose reaction data structures needed for in-memory cross section perturbation workflows, including:
These changes enable nuclear data perturbations to be performed directly in memory without requiring regenerated HDF5 libraries.
Motivation
The primary motivation is support for transport sensitivity and uncertainty analysis workflows where reaction cross sections may be repeatedly perturbed during runtime.
Implementation Notes
The implementation extends the shared-library Python bindings to expose reaction-level nuclide data and provide mechanisms for rebuilding internal nuclide data structures after modification.
Testing
Verified retrieval and modification of reaction cross sections through the Python API
Verified updated cross sections are used during transport
Existing regression tests pass
Example
`import openmc
import openmc.lib
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
def build_thresholded_xs(E, xs_reaction, threshold_energy):
"""
Build an XS array aligned to the unionized grid E,
with zeros below threshold_energy.
"""
xs_full = np.zeros_like(E, dtype=xs_reaction.dtype)
============================================================
Materials
============================================================
metal
metal = openmc.Material(name="Iron")
metal.add_element("Fe", 1.0)
metal.set_density("g/cm3", 7.87)
materials = openmc.Materials([metal])
materials.export_to_xml()
============================================================
Geometry
============================================================
Sphere
sphere = openmc.Sphere(r=40.0, boundary_type='vacuum')
region = -sphere # inside sphere
cell = openmc.Cell(fill=metal, region=region)
geom = openmc.Geometry([cell])
geom.export_to_xml()
============================================================
Fixed Source
============================================================
monoenergetic source at the origin
source = openmc.Source()
source.space = openmc.stats.Point((0.0, 0.0, 0.0))
source.energy = openmc.stats.Discrete([2.4e6], [1.0])
source.angle = openmc.stats.Isotropic()
============================================================
Settings
============================================================
settings = openmc.Settings()
settings.run_mode = "fixed source"
settings.source = source
settings.batches = 1
settings.particles = int(2e5)
settings.export_to_xml()
============================================================
Tallies
============================================================
tally = openmc.Tally(tally_id=1, name="heating")
tally.scores = ["heating"]
tally.filters = [openmc.CellFilter(cell)]
============================================================
Model
============================================================
model = openmc.Model(geom, materials, settings)
model.tallies.append(tally)
model.export_to_xml()
============================================================
Run normal OpenMC (unperturbed)
============================================================
print("\n--- Running base case ---\n")
openmc.lib.init()
openmc.lib.simulation_init()
tally_id = next(iter(openmc.lib.tallies))
openmc.lib.tallies[tally_id].reset()
openmc.lib.run()
result = openmc.lib.tallies[tally_id].mean[0][0]
print(f"heating = {result}, for MT base case")
============================================================
Access cross section data
============================================================
nuclide_name = 'Fe56'
mt = 2 # elastic scattering
T_index = 0
Get index of nuclide
idx = next(i for i, n in enumerate(openmc.lib.nuclides) if n == nuclide_name)
print(f"{nuclide_name} index = {idx}")
mts = openmc.lib.get_nuclide_mt_numbers(idx)
print(f"Nuclide {nuclide_name} has MT numbers: {mts}")
Retrieve xs grid
xs = openmc.lib.get_nuclide_xs(idx, mt, T_index)
Retrieve energy grid
E = openmc.lib.nuclide_energy_grid(idx, T_index)
threshold = openmc.lib.get_reaction_threshold_energy(idx, mt, T_index)
print(f"For {nuclide_name} and MT = {mt} the threshold is {threshold} eV")
xs_lat = build_thresholded_xs(E, xs, threshold)
============================================================
Apply perturbation
============================================================
xs_pass_to_openmc = xs_lat * 1.2
openmc.lib.new_nuclide_xs_with_threshold(
index=idx,
mt=mt,
energy=E,
xs_perturbed=xs_pass_to_openmc,
threshold_energy=threshold,
temperature=T_index
)
Rebuild derived cross sections to apply change properly (manual call)
openmc.lib.nuclide_rebuild_derived_xs(idx)
print("\nApplied a perturbation to cross section in memory and rebuilt derived XS.\n")
Retrieve xs grid
xs_new = openmc.lib.get_nuclide_xs(idx, mt, T_index)
xs_lat_new = build_thresholded_xs(E, xs_new, threshold)
============================================================
Plot before/after comparison
============================================================
plt.figure(figsize=(6,4))
plt.loglog(E, xs_lat, label='Original')
plt.loglog(E, xs_lat_new, label='From memory', ls='--')
plt.loglog(E, xs_pass_to_openmc, label='passed to memory', ls='none', marker='+')
plt.xlabel('Energy [eV]')
plt.ylabel('σ (barns)')
plt.title(f'{nuclide_name} MT={mt} Perturbation')
plt.grid(True, which='both', ls='--', lw=0.5)
plt.legend()
plt.tight_layout()
plt.show()
============================================================
Run with perturbed data
============================================================
tally_id = next(iter(openmc.lib.tallies))
openmc.lib.tallies[tally_id].reset()
openmc.lib.run()
result = openmc.lib.tallies[tally_id].mean[0][0]
print(f"heating = {result}, for MT {mt}")
openmc.lib.simulation_finalize()
openmc.lib.finalize()
print("\n--- Perturbed case complete ---\n")
`
Checklist