Photon transport discrepancies with MCNP and Serpent

Validation and Verification
1

Hello everyone,

I am trying to compare the results for a fixed source photon transport problem between OpenMC, MCNP, and Serpent.

Geometry:
Bare 7cm-radius sphere of U-235.

Source:
The source should emit photons with energies uniformly distributed between 9 and 10 MeV.

Results: (Energies expressed in MeV)
compare2

The integral fluxes are:

  • Serpent: 0.0018058876355745788
  • MCNP: 0.0022937391129849998
  • Integral OpenMC: 0.0020443844006718622

And the ratios are:

  • MCNP/OpenMC: 1.121970561031081
  • Serpent/OpenMC: 0.8833405473946562

Notes:

  • OpenMC is the only one that shows photons below 1keV, which surprises me. The simulation energy cutoff default value is 1keV. Why am I seeing photons below those energies? Serpent and MCNP have the same cutoff and their tallies don’t show photons below the cutoff.
    I tried defining a cutoff for OpenMC as
settings.cutoff = {'energy': 1000.0}

and

settings.cutoff = {'energy_photon': 1000.0}

with no effect on the results.

  • Any ideas why there are discrepancies in the integral fluxes? If I reduce the energy of the source, the results improve considerably. So, I am tempted to think that there might be some discrepancies in the physics models at high energies.

Thank you,
Roberto

2

If I reduce the source energies to 10 - 100 keV, the results improve:
compare2

And the integral fluxes are:

  • Serpent: 6.6645181798415495e-06
  • MCNP: 6.667402350129999e-06
  • OpenMC: 6.716380316594296e-06

And the ratios are:

  • MCNP/OpenMC: 0.9927076841757627
  • Serpent/OpenMC: 0.9922782608625349

Notes:

  • OpenMC result is the only one that shows photons above 100 keV. Any ideas what could be the cause?

Thank you, again.
Roberto

3

Update: I was able to match Serpent and MCNP spectra, I made a mistake with the sampling of the source position.
What would be the proper way to sample the source position uniformly within a cell in OpenMC?

I still don’t know why I see photons below 1 keV and above 100 keV in the second Figure.

4

@froberto Are you using a particle filter in your tally? If not, you may be seeing contributions from electrons/positrons that are generated due to photon interactions. Are you able to share your full model so that I can try to reproduce?

Regarding sampling the source position uniformly, you can use:

lower_left = (-7, -7, -7)
upper_right = (7, 7, 7)
spatial_dist = openmc.stats.Box(lower_left, upper_right)
source = openmc.Source(space=spatial_dist, ...)

This will result in source sites being generated outside your sphere, but those points will be rejected, leaving only points inside the sphere. In our upcoming release of OpenMC, you can properly sample points uniformly in a sphere with:

r = openmc.stats.PowerLaw(0, radius, 1)
theta = openmc.stats.Uniform(0, pi)
phi = openmc.stats.Uniform(0, 2*pi)
spatial_dist = openmc.stats.SphericalIndependent(r, theta, phi)
source = openmc.Source(space=spatial_dist, ...)
5

Thanks, @paulromano for your response.
I added your suggestions to my input file, but I still see undesired peaks below the energy cutoff.
I am probably doing something wrong here.

Here is my input file:

import openmc
import numpy as np


# MATERIALS
openmc.Materials.cross_sections = '/home/roberto/Documents/scratch/openmc/xsdata/endfb71_hdf5/cross_sections.xml'

uranium = openmc.Material(name="uranium")
uranium.add_nuclide('U235', 1.0)  # atom fraction
# uranium.set_density('g/cm3', 20.0)
uranium.set_density('atom/b-cm', 0.0506062)

materials = openmc.Materials([uranium])
materials.export_to_xml()


# GEOMETRY
radius = 7.0
fuel_outer_radius = openmc.Sphere(r=radius, boundary_type='vacuum')
fuel_region = -fuel_outer_radius

fuel = openmc.Cell(name='fuel')
fuel.fill = uranium
fuel.region = fuel_region

root_universe = openmc.Universe(cells=[fuel])
geometry = openmc.Geometry()
geometry.root_universe = root_universe
geometry.export_to_xml()


# SETTINGS
source = openmc.Source()
source.particle = 'photon'
source.energy = openmc.stats.Uniform(9e6, 1e7)

lower_left = (-7, -7, -7)
upper_right = (7, 7, 7)
source.space = openmc.stats.Box(lower_left, upper_right)

settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.photon_transport = True
settings.source = source
settings.batches = 200
settings.particles = 5000
settings.cutoff = {'energy': 1000.0}
settings.export_to_xml()


# TALLIES
cell_filter = openmc.CellFilter(fuel)
energies = np.logspace(2, 7, 501)
energy_filter = openmc.EnergyFilter(energies)  # in eV
particle_filter = openmc.ParticleFilter(['photon'])

tally = openmc.Tally(1)
tally.filters = [cell_filter, energy_filter, particle_filter]
tally.scores = ['flux']

tallies = openmc.Tallies([tally])
tallies.export_to_xml()


# RUN OPENMC
openmc.run()
6

Thanks for sharing your model @froberto. After investigating, it looks like the photons below 1 keV are the result of atomic relaxation (i.e., they are X-rays resulting from electrons being kicked out of the atom). Our cutoffs are not catching these photons right now, and so they end up going on to have a single collision before being killed. I’ll submit a bug fix to remedy this behavior. The spectrum above 1 keV should be OK though – are you seeing agreement above 1 keV for OpenMC vs MCNP vs Serpent?

Edit: bug fix has been submitted.

7

Thanks @paulromano for looking into this.
Yes, I see agreement between the OpenMC, MCNP, and Serpent.
One small caveat is that MCNP doesn’t resolve well the peaks around 10 keV, but OpenMC and Serpent do.

I have another question about sampling the source position.
Can I use the method described above if my geometry consists of two concentric spheres and the source is in one of them only?
Is there a way to choose the cell where the particles are born?
I’ve seen this other post.
I think I need something along those lines.

8

Glad to hear you’re seeing good agreement. Regarding the source question, yes, if we had a general source filtering capability, that would make it easy as you could just ask it to filter on the cell you want the source to be in. In the absence of that, if you wanted to sample a source uniformly over a spherical shell (one of your concentric spheres), you’d have to either write a C++ custom source routine (described here) or use the new PowerLaw distribution as follows:

r = openmc.stats.PowerLaw(r_min, r_max, 1)
theta = openmc.stats.Uniform(0, pi)
phi = openmc.stats.Uniform(0, 2*pi)
spatial_dist = openmc.stats.SphericalIndependent(r, theta, phi)
source = openmc.Source(space=spatial_dist, ...)
9

Great! Thank you for your time.