Dose discrepancy in Jupyter from OpenMC 14.0 to 15.0

User Support
1

Hello all

I was working with a colleague and noticed a strange discrepenacy between our outputs. We run the same code on a similar setup computationally, but I run OpenMC 14.0 and they run OpenMC 15.0. I was able to replicate this on my machine where a separate environment with 15.0.

I will attach the notebook at the end.
14.0 output: 408.147 +/- 54.948 mSv/yr
15.0 output: 343.0575 +/- 38.614 mSv/yr

Does anyone have ideas for what may be causing this difference? I’m planning on doing a simple attenuation calculation soon and comparing it with the current results but I am curious if anyone else runs into the same issue or if this has to do with a known issue or error on my part. I took a quick look at the 15.0 release notes and didn’t see anything that looked related on first glance.

I also just noticed as I was uploading that this error does not occur in typical .py code, just when put into .ipynb. I will upload the code shortly but want to keep this section relatively clear.

2

Code:

%%

import openmc
import numpy as np
import time
import neutronics_material_maker as nmm

NOTE - CHANGE CROSS SECTION PATH TO YOURS

openmc.config[‘cross_sections’] = ‘../../Libraries/endfb-viii.0-hdf5/cross_sections.xml’

%%

print(openmc.version)

%%

MATERIAL DEFINITONS

concrete_name = ‘Concrete, Regulatory Concrete (developed for U.S. NRC)’
air_mat = nmm.Material.from_library(‘Air (dry, near sea level)’).openmc_material
concrete_mat = nmm.Material.from_library(concrete_name).openmc_material
tissue_mat = nmm.Material.from_library(‘Tissue, Testis (ICRU)’).openmc_material
all_materials = openmc.Materials([concrete_mat,
tissue_mat])

%%

GEOMETRY

adding bounding box/universe

universe_machine = openmc.Universe()
bounding_surface = openmc.model.RectangularParallelepiped(xmin=-100, xmax=100,
ymin=-100, ymax=100,
zmin=-100, zmax=1000,
boundary_type = ‘vacuum’)
room_region = -bounding_surface
room_cell = openmc.Cell(region=room_region,
cell_id=100)
universe_machine.add_cells([room_cell])

Surface

concrete_thickness = 25 ## cm
test_human_radius = 10 ## cm
concrete_start = openmc.ZPlane(z0=100,
surface_id = 1000)
concrete_end = openmc.ZPlane(z0=concrete_start.z0 + concrete_thickness,
surface_id = 1001)
test_human = openmc.YCylinder(x0=0,
z0=concrete_end.z0+test_human_radius,
r=test_human_radius,
surface_id=1002)

Regions

concrete_region = (+concrete_start &
-concrete_end &
-bounding_surface)
test_human_region = (-test_human &
-bounding_surface)

Cells

concrete_cell = openmc.Cell(region = concrete_region,
cell_id = 200,
name = ‘Concrete Slab ala Shake Shack’)
concrete_cell.fill = concrete_mat
universe_machine.add_cells([concrete_cell])
room_region &= ~concrete_region

test_human_cell = openmc.Cell(region = test_human_region,
cell_id = 201,
name = ‘Test Human ala RadMan’)
test_human_cell.fill = tissue_mat
universe_machine.add_cells([test_human_cell])
room_region &= ~test_human_region

Define full geometry

full_geometry = openmc.Geometry(root = universe_machine,
merge_surfaces=True,
surface_precision=1)
full_geometry.root_universe.plot(basis = ‘xz’,
width=(250,1400),
origin= (0,0,400),
pixels=(700,700),
color_by = ‘cell’)

%%

Setting up the source

source_name = ‘DD’

point_source = openmc.IndependentSource()
point_source.particle = ‘neutron’
if source_name == ‘DD’:
point_source.energy = openmc.stats.Discrete([2.45e6],[1])
elif source_name == ‘DT’:
point_source.energy = openmc.stats.Discrete([14.1e6],[1])
else:
print(‘Invalid Source’)

point_source.angle = openmc.stats.Isotropic()
point_source.space = openmc.stats.Point(xyz=(0,0,0))
point_source.strength = 1

%%

Setting up dose within a cell

energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(
particle=“neutron”, geometry=“AP”
)
energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n, interpolation=‘cubic’)
#energy_function_filter_n.interpolation = “cubic” # cubic interpolation is recommended by ICRP

neutron_particle_filter = openmc.ParticleFilter(“neutron”)
cell_filter = openmc.CellFilter(test_human_cell)

dose_cell_tally = openmc.Tally(name=“neutron_dose_on_cell”)

note that the EnergyFunctionFilter is included as a filter

dose_cell_tally.filters = [
cell_filter,
neutron_particle_filter,
energy_function_filter_n,
]
dose_cell_tally.scores = [“flux”]
my_tallies = openmc.Tallies([dose_cell_tally])

%%

Simulation Settings

new_settings = openmc.Settings()
new_settings.run_mode = “fixed source”
new_settings.particles = int(1e4)
new_settings.batches = 3
new_settings.output = {‘tallies’:False,
‘path’:‘Data/’,
‘summary’:False}
new_settings.electron_treatment = ‘led’
new_settings.photon_transport = True
new_settings.source = point_source
new_settings.seed = int(time.time())

model = openmc.Model(full_geometry,
all_materials,
new_settings,
my_tallies)
model.geometry.merge_surfaces = True

model.export_to_xml()

%%

start_time = time.time()

!rm s*.h5

model.run()
end_time = time.time()

# mv

print((end_time - start_time)/60)

%%

with openmc.StatePoint(‘Data/statepoint.3.h5’) as statepoint:

neutron_tally_result = statepoint.get_tally(
    name="neutron_dose_on_cell"
).mean.flatten()[0]
neutron_tally_error = statepoint.get_tally(
    name="neutron_dose_on_cell"
).std_dev.flatten()[0]

print(neutron_tally_result)
print(neutron_tally_error)

neutrons_per_second = 1e8 # units of neutrons per second

tally.mean is in units of pSv-cm3/source neutron

this multiplication changes units to neutron to pSv-cm3/second

total_dose = neutron_tally_result * neutrons_per_second

converts from pSv-cm3/second to pSv/second

total_dose = total_dose / (np.pitest_human_radiustest_human_radius*200)

converts from (pico) pSv/second to (milli) mSv/second

total_dose = total_dose * 1e-9

converts from (milli) mSv/second to (milli) mSv/year

total_dose = total_dose * 60 * 60 * 24 * 365

tally.mean is in units of pSv-cm3/source neutron

this multiplication changes units to neutron to pSv-cm3/second

total_dose_error = neutron_tally_error * neutrons_per_second

converts from pSv-cm3/second to pSv/second

total_dose_error = total_dose_error / (np.pitest_human_radiustest_human_radius*200)

converts from (pico) pSv/second to (milli) mSv/second

total_dose_error = total_dose_error * 1e-9

converts from (milli) mSv/second to (milli) mSv/year

total_dose_error = total_dose_error * 60 * 60 * 24 * 365
print(f’{total_dose} +/- {total_dose_error}')

3

Hi AOlive,
I can’t check directly as I don’t have access to a version with 14.0 installed, nor the model you’re pulling from, but I would guess that either you’re using different cross section databases, or its something to do with the way you export the xmls from a model and run it. You are exporting to XML which is producing xml files for geometry, materials, and settings, and then running mode.run(), which is essentially doing the same thing but then also running it afterwards. I have a feeling that this is where it could be causing issues. I would try to just export the individual XMLs for geometry, materials, settings, and tallies, using their respective commands, and then use openmc.run to see if your results vary. Make sure no model is within the same folder as the XML’s as OpenMC (at least for 15.2 using jupyterlabs which is how I run all my work) will take the model over individual XML’s if it sees both. You can actually see this at the beginning of the of the openmc.run, it will say it found a model and is running it instead. Let me know if this helped you!

4

Thank you for your quick response!

I will say the cross section path is set at the top of the script so there should be consistency between the 14.0 and 15.0 with that. I understand it is very hard to read in the format I put it in, but as a new user could not put it in a more readable format, apologies for that. I tried rerunning it with using export_to_xml() and openmc.run(), which yielded the same result (I did remove all the XML files from the directory before both runs.)
Aside from changing the cross section path at the top of the code, I think it should run and generate the model without any new files needed so it is (hopefully) easy to run. I am still unable to attach files, but if this thread continues to be active I will upload it in .ipynb format when possible.

5

No problem, it’s not intuitive, to include the code in a more readable way you can paste the text and then make sure its all on its own line by pushing enter, and then drag over it all and push tab so that its all indented. If it works you should see on the preview on the right that it’s working and all in a special format. You are saying error, which result do you think is the correct one? Did you try the export to xml using both 14.0 and 15.0? as well as in .py and .ipynb? and still get different results?

6

I am unable to go back and edit the first code snippet but here is the fixed version along with the dose total and error for each:

import openmc
import numpy as np
import time
import neutronics_material_maker as nmm

# NOTE - CHANGE CROSS SECTION PATH TO YOURS
openmc.config['cross_sections'] = '../../Libraries/endfb-viii.0-hdf5/cross_sections.xml'

print(openmc.__version__)

## MATERIAL DEFINITONS
concrete_name = 'Concrete, Regulatory Concrete (developed for U.S. NRC)'
air_mat = nmm.Material.from_library('Air (dry, near sea level)').openmc_material
concrete_mat = nmm.Material.from_library(concrete_name).openmc_material
tissue_mat = nmm.Material.from_library('Tissue, Testis (ICRU)').openmc_material
all_materials = openmc.Materials([concrete_mat,
                                 tissue_mat])

## GEOMETRY

# adding bounding box/universe
universe_machine = openmc.Universe()
bounding_surface = openmc.model.RectangularParallelepiped(xmin=-100, xmax=100, 
                                                          ymin=-100, ymax=100, 
                                                          zmin=-100, zmax=1000, 
                                                          boundary_type = 'vacuum')
room_region = -bounding_surface
room_cell = openmc.Cell(region=room_region,
                       cell_id=100)
universe_machine.add_cells([room_cell])

# Surface
concrete_thickness = 25 ## cm
test_human_radius = 10 ## cm
concrete_start = openmc.ZPlane(z0=100,
                              surface_id = 1000)
concrete_end = openmc.ZPlane(z0=concrete_start.z0 + concrete_thickness,
                            surface_id = 1001)
test_human = openmc.YCylinder(x0=0,
                             z0=concrete_end.z0+test_human_radius,
                             r=test_human_radius,
                             surface_id=1002)

# Regions
concrete_region = (+concrete_start & 
                  -concrete_end &
                   -bounding_surface)
test_human_region = (-test_human &
                    -bounding_surface)

# Cells
concrete_cell = openmc.Cell(region = concrete_region,
                           cell_id = 200,
                           name = 'Concrete Slab ala Shake Shack')
concrete_cell.fill = concrete_mat
universe_machine.add_cells([concrete_cell])
room_region &= ~concrete_region

test_human_cell = openmc.Cell(region = test_human_region,
                             cell_id = 201,
                             name = 'Test Human ala RadMan')
test_human_cell.fill = tissue_mat
universe_machine.add_cells([test_human_cell])
room_region &= ~test_human_region

## Define full geometry
full_geometry = openmc.Geometry(root = universe_machine,
                               merge_surfaces=True,
                               surface_precision=1)
full_geometry.root_universe.plot(basis = 'xz', 
                                 width=(250,1400),
                                 origin= (0,0,400),
                                 pixels=(700,700), 
                                 color_by = 'cell')

## Setting up the source

source_name = 'DD'

point_source = openmc.IndependentSource()
point_source.particle = 'neutron'
if source_name == 'DD':
    point_source.energy = openmc.stats.Discrete([2.45e6],[1])
elif source_name == 'DT':
    point_source.energy = openmc.stats.Discrete([14.1e6],[1])
else:
    print('Invalid Source')
    
point_source.angle = openmc.stats.Isotropic()
point_source.space = openmc.stats.Point(xyz=(0,0,0))
point_source.strength = 1

## Setting up dose within a cell

energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(
    particle="neutron", geometry="AP"
)
energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n, interpolation='cubic')
#energy_function_filter_n.interpolation = "cubic"  # cubic interpolation is recommended by ICRP

neutron_particle_filter = openmc.ParticleFilter("neutron")
cell_filter = openmc.CellFilter(test_human_cell)

dose_cell_tally = openmc.Tally(name="neutron_dose_on_cell")
# note that the EnergyFunctionFilter is included as a filter
dose_cell_tally.filters = [
    cell_filter,
    neutron_particle_filter,
    energy_function_filter_n,
]
dose_cell_tally.scores = ["flux"]
my_tallies = openmc.Tallies([dose_cell_tally])

## Simulation Settings
new_settings = openmc.Settings()
new_settings.run_mode = "fixed source"
new_settings.particles = int(1e5)
new_settings.batches = 10
new_settings.output = {'tallies':False,
                      'summary':False}
new_settings.electron_treatment = 'led'
new_settings.photon_transport = True
new_settings.source = point_source
new_settings.seed = int(time.time())

model = openmc.Model(full_geometry,
                    all_materials,
                    new_settings,
                    my_tallies)
model.geometry.merge_surfaces = True

model.export_to_xml()

start_time = time.time()
# !rm s*.h5

model.run()
end_time = time.time()

# # mv
print((end_time - start_time)/60)

with openmc.StatePoint(f'statepoint.{new_settings.batches}.h5') as statepoint:

    neutron_tally_result = statepoint.get_tally(
        name="neutron_dose_on_cell"
    ).mean.flatten()[0]
    neutron_tally_error = statepoint.get_tally(
        name="neutron_dose_on_cell"
    ).std_dev.flatten()[0]

print(neutron_tally_result)
print(neutron_tally_error)

neutrons_per_second = 1e8  # units of neutrons per second

# tally.mean is in units of pSv-cm3/source neutron
# this multiplication changes units to neutron to pSv-cm3/second
total_dose = neutron_tally_result * neutrons_per_second
# converts from pSv-cm3/second to pSv/second
total_dose = total_dose / (np.pi*test_human_radius*test_human_radius*200)
# converts from (pico) pSv/second to (milli) mSv/second
total_dose = total_dose * 1e-9
# converts from (milli) mSv/second to (milli) mSv/year
total_dose = total_dose * 60 * 60 * 24 * 365

# tally.mean is in units of pSv-cm3/source neutron
# this multiplication changes units to neutron to pSv-cm3/second
total_dose_error = neutron_tally_error * neutrons_per_second
# converts from pSv-cm3/second to pSv/second
total_dose_error = total_dose_error / (np.pi*test_human_radius*test_human_radius*200)
# converts from (pico) pSv/second to (milli) mSv/second
total_dose_error = total_dose_error * 1e-9
# converts from (milli) mSv/second to (milli) mSv/year
total_dose_error = total_dose_error * 60 * 60 * 24 * 365
print(f'{total_dose} +/- {total_dose_error}')

15.0 and .py: 458.7325747958416 +/- 12.754034265550862
14.0 and .py: 482.30560557190455 +/- 9.111298901804613
15.0 and .ipynb: 463.94745477792696 +/- 8.028673977084425
14.0 and .ipynb: 30.88386888500472 +/- 0.5119372125286035
Note: The only things to change between these runs are an environment change (to move between the 14.0 and 15.0 versions) and removing XML files after each run.
Also, apologies for saying ‘error’ since I have not done the attenuation calculation yet but it seems like the only tally that disagrees to a large extent is 14.0 when run in .ipynb format. Due to two of my colleagues also having this error on their separate computers (of the same model) I am curious if it is a common environment error or an error with some segment of how it executes.

7

So I ran your simulation with better counting statistics and got the following result,

da558315b26554e368bffee7e45b3fe6
da558315b26554e368bffee7e45b3fe6729×388 45.1 KB

I ran using Py as well and got the same result. I would suggest increasing your counting statistics because the difference between your 15.0 using .py, 14.0 py, and 15.0 and ipynb are somewhat close enough and possibly due to low counting statistics and propagation of error. As for why your 14.0 ipynb is different I do not know, could be an environment setting that is present that you are unaware of.