Cubic interpolation of dose coefficients producing order of magnitude lower doses

[rlbarker]

Bug Description

When running a prompt photon dose simulation I have noticed very different results when using cubic interpolation of dose coefficients compared to all other types of interpolation. This effect is consistent whichever dose coefficients you use. The cubic interpolation provides a dose that is an order of magnitude lower than all other interpolation methods.

For example, I have attached a plot produced from a shutdown dose simulation of a steel sphere that was first irradiated with neutrons and the photon dose from the activated steel then plotted on a mesh. Results are consistent between different dose coefficients and interpolation methods (I used ICRP116 and ICRU57 + Pellicioni which Fluka uses) except for cubic interpolation.

Steps to Reproduce

I have put below a python script with a simple (2 spheres) geometry and a photon source. Dose tallies are performed for all 5 interpolation methods. All tally results agree the dose is roughly 18 pSvcm2, except the cubic interpolation tally which calculates 2 pSvcm2.

import openmc

# Geometry
material1 = openmc.Material(material_id=1)
material1.add_element(element="H", percent=1)

my_materials = openmc.Materials([material1])

sphere1 = openmc.Sphere(r=10)
sphere2 = openmc.Sphere(r=20, boundary_type='vacuum')

region1 = -sphere1
region2 = -sphere2 & +sphere1

cell1 = openmc.Cell(region=region1, fill=material1)
cell2 = openmc.Cell(region=region2)

my_geometry = openmc.Geometry([cell1, cell2])

# Source
source1 = openmc.Source()
source1.space = openmc.stats.Point((0, 0, 0))
source1.angle = openmc.stats.Isotropic()
source1.energy = openmc.stats.Discrete([0.3e6, 0.9e6, 2e6], [0.3, 0.4, 0.3])
source1.particle = 'photon'

# Settings
my_settings = openmc.Settings(
    batches=100,
    particles=10_000,
    run_mode="fixed source",
    source=source1
    )

# Tallies: one dose tally for each interpolation method
energies_p, pSv_cm2_p = openmc.data.dose_coefficients('photon', geometry="AP")
energy_function_filter_cubic = openmc.EnergyFunctionFilter(energies_p, pSv_cm2_p)
energy_function_filter_cubic.interpolation = "cubic"

energy_function_filter_log_log = openmc.EnergyFunctionFilter(energies_p, pSv_cm2_p)
energy_function_filter_log_log.interpolation = "log-log"

energy_function_filter_linear_linear = openmc.EnergyFunctionFilter(energies_p, pSv_cm2_p)
energy_function_filter_linear_linear.interpolation = "linear-linear"

energy_function_filter_linear_log = openmc.EnergyFunctionFilter(energies_p, pSv_cm2_p)
energy_function_filter_linear_log.interpolation = "linear-log"

energy_function_filter_log_linear = openmc.EnergyFunctionFilter(energies_p, pSv_cm2_p)
energy_function_filter_log_linear.interpolation = "log-linear"

cell_filter = openmc.CellFilter(cell2)
particle_filter = openmc.ParticleFilter('photon')

cubic_dose_tally = openmc.Tally(name='cubic_dose_tally')
cubic_dose_tally.filters = [
    cell_filter,
    particle_filter,
    energy_function_filter_cubic]
cubic_dose_tally.scores = ['flux']

log_log_dose_tally = openmc.Tally(name='log_log_dose_tally')
log_log_dose_tally.filters = [
    cell_filter,
    particle_filter,
    energy_function_filter_log_log]
log_log_dose_tally.scores = ['flux']

linear_linear_dose_tally = openmc.Tally(name='linear_linear_dose_tally')
linear_linear_dose_tally.filters = [
    cell_filter,
    particle_filter,
    energy_function_filter_linear_linear]
linear_linear_dose_tally.scores = ['flux']

linear_log_dose_tally = openmc.Tally(name='linear_log_dose_tally')
linear_log_dose_tally.filters = [
    cell_filter,
    particle_filter,
    energy_function_filter_linear_log]
linear_log_dose_tally.scores = ['flux']

log_linear_dose_tally = openmc.Tally(name='log_linear_dose_tally')
log_linear_dose_tally.filters = [
    cell_filter,
    particle_filter,
    energy_function_filter_log_linear]
log_linear_dose_tally.scores = ['flux']

my_tallies = openmc.Tallies([
    cubic_dose_tally,
    log_log_dose_tally,
    linear_linear_dose_tally,
    linear_log_dose_tally,
    log_linear_dose_tally
    ])

model = openmc.model.Model(
    geometry=my_geometry, 
    materials=my_materials, 
    settings=my_settings, 
    tallies=my_tallies
    )

output_file = model.run()

with openmc.StatePoint(output_file) as results:
    cubic_dose_tally = results.get_tally(
        name="cubic_dose_tally"
    )
    cubic_dose = cubic_dose_tally.get_pandas_dataframe()["mean"].sum()

    log_log_dose_tally = results.get_tally(
        name="log_log_dose_tally"
    )
    log_log_dose = log_log_dose_tally.get_pandas_dataframe()["mean"].sum()

    linear_linear_dose_tally = results.get_tally(
        name="linear_linear_dose_tally"
    )
    linear_linear_dose = linear_linear_dose_tally.get_pandas_dataframe()["mean"].sum()

    linear_log_dose_tally = results.get_tally(
        name="linear_log_dose_tally"
    )
    linear_log_dose = linear_log_dose_tally.get_pandas_dataframe()["mean"].sum()

    log_linear_dose_tally = results.get_tally(
        name="log_linear_dose_tally"
    )
    log_linear_dose = log_linear_dose_tally.get_pandas_dataframe()["mean"].sum()

print(f"cubic: {cubic_dose} pSv cm2")
print(f"log_log: {log_log_dose} pSv cm2") 
print(f"linear_linear: {linear_linear_dose} pSv cm2")
print(f"linear_log: {linear_log_dose} pSv cm2")
print(f"log_linear: {log_linear_dose} pSv cm2")

Environment

I am using the development branch of openmc version 0.13.4-dev. ENDF VIII data. Ubuntu OS.

[giovanni-mariano]

Hello, I was facing the same issue some days ago. I think the bug is in the interpolate_lagrangian function:

openmc/include/openmc/interpolate.h

Line 50 in a3695c7

output += (numerator / denominator) * ys[i];

in particular, I think that "ys[i]" should be ys[idx + i] but I still haven't done any test.

Giovanni

[pshriwise]

Thanks for reporting this @rlbarker! And @giovanni-mariano I believe your fix is correct. 👍🏻

I'll submit a PR and add a test to our C++ test suite to ensure the behavior of the interpolate_lagrangian is correct going forward.

[rlbarker]

Thanks everyone!

[pshriwise]

Thanks for the example above @rlbarker. After the fix in #2775 things look better. If those plots above are relatively easy to reproduce, mind posting them after a re-run with that fix?

cubic: 18.110131567651347 pSv cm2
log_log: 18.098182424892507 pSv cm2
linear_linear: 18.09241280539433 pSv cm2
linear_log: 18.155939740734116 pSv cm2
log_linear: 18.034918963809766 pSv cm2

[rlbarker]

Here are the updated plots - looks much better! Thank you

[pshriwise]

Fantastic. Thanks @rlbarker!