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)

[jtramm]

This is very cool! The general idea for the implementation looks great as well -- seems like an elegant way of getting this done.

I'm also interested in hearing more about the workflows you have in mind:

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

[gridley]

Hey John! As for the workflow in mind, I have been working on some reactor shielding problems that have both gammas and neutrons to shield both coming off a surface source.

Maybe I'm missing something about the mainstream workflows for shielding, but I got the impression that if the quantity of interest is a dose response function at some point beyond the shield, it's not known ahead of time whether neutrons or gammas are the dominant contributor to dose. Since the overall dose rate is a functional over both the neutron and gamma flux spectrum, then to get a flat estimate of that via FW-CADIS you would need to have a multigroup formulation of neutrons and gammas put together. I think that they can be decoupled, since gammas don't make neutrons (to first order).

So in the forward problem you can solve for neutron transport first, then get the gamma transport problem solved including the neutron-to-gamma secondary source. In the adjoint version, it should be possible to first solve the adjoint gamma transport problem, then use that as a source contribution in the multigroup neutron transport problem. As far as I can tell there is no capability like this in the current weight window generator.

So, in conclusion, getting that off-diagonal block for gamma production from neutrons in multigroup is what this PR is meant to feed into.

[gridley]

Alright, @jtramm, I think this is good to go. I can add more to this if you like, but am just working towards that dream of automatic weight windows for coupled n/gamma problems bit by bit. The next step would be generating the multigroup photon library. With this filter in place, I think all the C++ code is there to do that.

[shimwell]

Looks great, I think one of the regression tests is failing due to the change of an input file which now includes more filters.

I am also interested in seeing the FW-Cadis weight window generator used for photons. In my case the prompt photons contribute less that 1 percent of the neutron dose but I would be interested in weight windows for decay photons (which I think we can't currently do in openmc)

[gridley]

@shimwell do you understand why that test would fail? I updated the inputs_true.dat and results_true.dat files, and have the test passing locally. Is there some other location where inputs_true files are kept which also needs to be updated?

[gridley]

Oh, I see now! There is another test using that same photon production test on model.xml styled input. The new commit has updated that one as well.