Tally spectrum of secondary particles

[GuySten]

Description

This PR enable tallying spectrum of secondary particles.

This PR is a different approach to #3453 .

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)

[GuySten]

Currently (in the new approach), the example problem from #3453 gives:

Which looks the same as in #3453.

The script to run this is:

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'])
particleout_filter = openmc.ParticleoutFilter(['photon'])
energyout_filter = openmc.EnergyoutFilter(energy_bins)

# Define tally
tally = openmc.Tally(name='photon spectrum')
tally.filters = [particle_filter, particleout_filter, energyout_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()

[gridley]

So, to accomplish the functionality I was shooting for in the original PR, the goal would be to use a matrix for group-to-group production of different particle types. For example columns might be photon energy groups and rows neutron groups so that the multigroup spectrum of photons given neutrons could be obtained.

With this PR, it's not clear to me how that could be obtained in a tensor product fashion like that, although it may be able to calculate the overall secondary spectrum. You need to capture the correlation between neutron E and photon E, and with this approach where tallying seems to not take place until the secondary is revived and treated as an active particle, I'm not seeing how that correlation could be preserved for output in a tally since the parent particle's energy long gone after the secondary is revived.

[GuySten]

You're right, to do that I will need to store the original particle energy in E_last...

I think afterward you will get your matrices by using a combination of EnergyFilter, EnergyoutFilter, ParticleFilter and ParticleOutFilter.

[gridley]

OK, I see what you're getting at. But suppose in particular we want the photon-from-neutron production matrices' Legendre moments. In that case you have to go and start saving the pre-collision direction of the original particle that made the source site. On the other hand, a SecondaryDirectionCosine could be filtered instead. If other quantities were needed (I can't think of any at the moment), the data saved on the source site (a struct whose size is important to keep small) again has to be increased further.

So, for these reasons, my opinion is that the modified filters on secondaries are superior because they never require storing extra data -- the tallying gets done when the sites are created. It's more code, but we don't have to add any transport-time data tied to source sites.

[GuySten]

Thanks for your feedback @gridley.

Based on your comments:
I think that I should score at sourcesite creation when all data is available (and not after the particle died).
If I want to use the current filters mechanisms I should generate particles for tally purposes and then drop them immediately.
That way the scoring method could access all available data and I won't need to increase sourcesite size.

I will try the new approach and will get back to you.

[GuySten]

In the new approach the tallying gets done when creating secondaries.
A temporary Particle is constucted and then passed to the scoring function.
As a bonus the amount of new code got smaller.

[gridley]

Now that is a very interesting and effective approach! I would counter this by pointing out that putting a Particle on the stack is going to allocate memory in the backend for geometry state, last I checked. Simply calling malloc when particles are in flight can be a pretty tangible performance penalty. I believe it's also doing a malloc for the cached microscopic cross-sections.

At this point, I'd like to leave it to someone else to decide on which implementation, if any, they'd like to merge into the code.

[GuySten]

Fair enough @gridley.

@paulromano, @pshriwise, @jtramm.

What do you think?