Secondary energy filter

[gridley]

Description

Towards being able to generate multigroup photon/neutron libraries for shielding purposes, with the eventual goal of random ray adjoints for dose in shielding contexts, I have implemented a filter on secondary particle energies. The first step is to be able to generate the multigroup photon source vector from a neutron flux vector for a given MGXS material. This can be done by tallying secondary particle energies generated during a collision.

To implement this, I added a counter on the secondaries created in a collision. Analog estimation can then be used to bin the secondary particle energies. Naturally, to ensure that secondary photons from photons are not captured, this requires specifying the primary and secondary particle type on the energy filter.

I want to check this by comparing the photon source from an n/p simulation, and using it as the source in a plain photon spectrum to check they match up.

Here is an example of using this feature to find the secondary photon specrum emitted by 2 MeV source neutrons in U238. There are about 20 photons generated per source neutron.

import openmc
import numpy as np
import matplotlib.pyplot as plt

# Define material: pure U-238
u238 = openmc.Material(name='U-238')
u238.add_nuclide('U238', 1.0)
u238.set_density('g/cm3', 19.1)

materials = openmc.Materials([u238])

# Define geometry: large U-238 sphere (radius = 1e90 cm)
sphere = openmc.Sphere(r=1e90, boundary_type='vacuum')
cell = openmc.Cell(fill=u238, region=-sphere)
universe = openmc.Universe(cells=[cell])
geometry = openmc.Geometry(universe)

# Define fixed 2 MeV neutron source at origin
source = openmc.Source()
source.space = openmc.stats.Point()
source.energy = openmc.stats.Discrete([2.0e6], [1.0])
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.particles = 10_000
settings.batches = 10
settings.source = source
settings.photon_transport = True

# Define secondary energy filter for photons
energy_bins = np.logspace(4, 8, 100)  # 10 keV to 100 MeV
particle_filter = openmc.ParticleFilter(['neutron'])
sec_energy_filter = openmc.SecondaryEnergyFilter(energy_bins, particle='photon')

# Define tally
tally = openmc.Tally(name='photon spectrum')
tally.filters = [particle_filter, sec_energy_filter]
tally.scores = ['total']
tally.estimator = 'analog'

tallies = openmc.Tallies([tally])

# Create model
model = openmc.model.Model(geometry, materials, settings, tallies)

# Run simulation
sp_filename = model.run()

# Postprocess results
with openmc.StatePoint(sp_filename) as sp:
    t = sp.get_tally(name='photon spectrum')
    spectrum = t.mean.ravel()

# Plotting
bin_centers = 0.5 * (energy_bins[:-1] + energy_bins[1:])
plt.figure(figsize=(8, 5))
plt.loglog(bin_centers, spectrum, drawstyle='steps-mid', label="Secondary photon spectrum")

# Superimpose the excitation energies for inelastic scattering
u238 = openmc.data.IncidentNeutron.from_hdf5('/Users/gavin/Code/nndc_hdf5/U238.h5')
q_values = [-u238.reactions[mt].q_value for mt in range(51, 91)]
plt.vlines(q_values, 1e-2, 1, label="Inelastic -Q Values", color="orange")
plt.ylim([1e-2, 1])
plt.xlim([energy_bins.min(), energy_bins.max()])
plt.xlabel('Photon Energy [eV]')
plt.ylabel('Photons per Source Neutron')
plt.title('Secondary Photon Energy Spectrum (2 MeV n on U-238)')
plt.grid(True, which='both')
plt.legend()
plt.tight_layout()
plt.show()

Which gave me this result:

For now this is just a first pass to open up the implementation to feedback. I plan to add a regression test and further validation of the idea soon.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)