Different doses between v0.14 and 0.15

[eanutter]

Hi,

I’ve been looking at doses for a shielded reactor, and noticed some big differences between the dose results between v0.14 and v0.15 of OpenMC. The tally filters I was using were a mesh filter, a particle filter and an energy function filter.
I’ve managed to narrow down the difference in results to the effect of energy function filter, but can’t seem to find where the difference is coming from (about 80x as much neutron dose in v0.15 compared to v0.14 and 40x as much photon dose).

Has anyone else had a similar issue/can help diagnose the problem?

Thanks in advance.

[harrytili]

I’m having the same issue. Is anyone else experiencing this please?

[EdP]

On investigation, this issue appears to be related only to the cubic interpolation of the data.dose_coefficients function, when used with the EnergyFunctionFilter. Other interpolation schemes work fine.

It looks to me like there was potentially a bug before 0.15, however I can’t readily see any change in EnergyFunctionFilter between the two versions. Perhaps someone could confirm if this is related to a known PR… is #2887 is related somehow?

For full context, running the simple model below with 0.14 produces a tally value 60x smaller than running with 0.15. Other interpolation schemes are broadly equivalent.

""" DESCRIPTION: A simple cylinder of fuel"""
import openmc
import numpy as np

#### MATERIALS ####
uo2 = openmc.Material(1, "uo2")
uo2.add_nuclide('U235', 0.1)
uo2.add_nuclide('U238', 0.9)
uo2.add_nuclide('O16', 2.0)
uo2.set_density('g/cm3', 10.0)

mats = openmc.Materials([uo2, ])
mats.export_to_xml("materials.xml")

#### GEOMETRY ####
fuel_or = openmc.ZCylinder(r=10.0, boundary_type="vacuum")
top_plane = openmc.ZPlane(z0=10.0, boundary_type="vacuum")
bottom_plane = openmc.ZPlane(z0=-10.0, boundary_type="vacuum")
fuel_region = -fuel_or &+bottom_plane & -top_plane

fuel = openmc.Cell(name="fuel", fill =uo2, region = fuel_region)

geometry = openmc.Geometry([fuel, ])
geometry.export_to_xml("geometry.xml")

#### TALLIES ####
tallies = openmc.Tallies()
# Create dose coefficients 
energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(particle="neutron", geometry="AP", )
# Create a neutron energy filter
energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)

# ****
energy_function_filter_n.interpolation = "cubic"   # TOGGLE cubic/histogram FOR ERROR RE
# ****

# Create neutron tally to score dose
dose_cell_tally_n = openmc.Tally (name="neutron_dose_on_mesh")
dose_cell_tally_n.filters = [energy_function_filter_n, ]  
dose_cell_tally_n.scores = ['flux']
tallies.append(dose_cell_tally_n)
tallies.export_to_xml("tallies.xml")


#### SETTINGS ####
settings = openmc.Settings()
uniform_dist = openmc.stats.Box([-10.0, -10.0, -10.0], [10.0, 10.0, 10.0],)
src =openmc.IndependentSource(space=uniform_dist,)

settings.source = src
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
settings.export_to_xml("settings.xml")

## RUN! ##
openmc.run(threads=40)


## POST-PROCESS ##
sp_name = "statepoint.100.h5"
sp = openmc.StatePoint(sp_name)
tally_output = sp.get_tally(name="neutron_dose_on_mesh")    
tally_val = tally_output.mean.flatten()

print("---tally:")
print(tally_val)
print("---keff:")
print("{:5.4f}".format(unumpy.nominal_values(sp.keff)))

[pshriwise]

Hi all,

This was indeed a bug related to cubit interpolation for the EnergyFunctionFilter. Please see the following issue and fix below:

Cubic interpolation of dose coefficients producing order of magnitude lower doses

Best,

Patrick