Neutron flux in a graphite stack

User Support
1

Hi All,

I’m quite new to OpenMC and I’m attempting to replicate an experiment I have carried out in OpenMC. The experiment is an AmBe neutron source in middle-bottom of a large graphite stack.

import openmc

# create graphite - 100% carbon, density 2.26g/cm^3
graphite = openmc.Material()
graphite.add_element('C', 1.0)
graphite.set_density('g/cm3', 2.26)
materials = openmc.Materials([graphite])

# cuboid - simple graphite cube for now
graphite_surface = openmc.model.RectangularParallelepiped(-70.0, 70.0, -70.0, 70.0, 0, 210.0, boundary_type='vacuum')
graphite_cell = openmc.Cell(fill=graphite, region=-graphite_surface)

# universe with graphite cell - i think this is right
root_universe = openmc.Universe(cells=[graphite_cell])

# geometry - i think this is right
geometry = openmc.Geometry(root_universe)

#source (AmBe) - confident on this part
source = openmc.IndependentSource()
source.space = openmc.stats.Point((0,0,10))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([4.2e6], [1]) #average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually


# settings - this works but not sure if its ideal
settings = openmc.Settings()
settings.source = source
settings.particles = 2500
settings.batches = 20
settings.inactive = 5
settings.output= {'tallies': True}

# figure out how to tally neutrons?

# export to xml
geometry.export_to_xml()
materials.export_to_xml()
settings.export_to_xml()

# run - is this what i want or do I want something else if i'm just trying to get neutron flux/counts
openmc.run()

Above is what I’ve made so far but not sure how to progress to what I want to get which is ideally a flux measurement throughout the graphite cube

Thanks,
Nooradin

2

Welcome to the community Nooradin

It looks like you would just need to add a tally to that model, please check for typos but this is roughly what would be added

here are some other flux tally examples

import openmc

# create graphite - 100% carbon, density 2.26g/cm^3
graphite = openmc.Material()
graphite.add_element('C', 1.0)
graphite.set_density('g/cm3', 2.26)
materials = openmc.Materials([graphite])

# cuboid - simple graphite cube for now
graphite_surface = openmc.model.RectangularParallelepiped(-70.0, 70.0, -70.0, 70.0, 0, 210.0, boundary_type='vacuum')
graphite_cell = openmc.Cell(fill=graphite, region=-graphite_surface)

# universe with graphite cell - i think this is right
root_universe = openmc.Universe(cells=[graphite_cell])

# geometry - i think this is right
geometry = openmc.Geometry(root_universe)

#source (AmBe) - confident on this part
source = openmc.IndependentSource()
source.space = openmc.stats.Point((0,0,10))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([4.2e6], [1]) #average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually


# settings - this works but not sure if its ideal
settings = openmc.Settings()
settings.source = source
settings.particles = 2500
settings.batches = 20
settings.inactive = 5
settings.output= {'tallies': True}

# figure out how to tally neutrons?
tally = openmc.Tally()
tally.score = ['flux']
mat_filter = openmc.MaterialFilter(graphite)
tally.filters = openmc.Filters = [mat_filter]
tallies = openmc.Tallies([tally])

# export to xml
geometry.export_to_xml()
materials.export_to_xml()
settings.export_to_xml()
tallies.export_to_xml()

# run - is this what i want or do I want something else if i'm just trying to get neutron flux/counts
openmc.run()
3

I think it needed ‘tally.scores =[‘flux’]’ but I seem to be getting the error:

RuntimeError: No fission sites banked on MPI rank 0 Abort(-1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0

Heres the updated code just in case

import openmc

# create graphite - 100% carbon, density 2.26g/cm^3
graphite = openmc.Material()
graphite.add_element('C', 1.0)
graphite.set_density('g/cm3', 2.26)
materials = openmc.Materials([graphite])

# cuboid - simple graphite cube for now
graphite_surface = openmc.model.RectangularParallelepiped(-70.0, 70.0, -70.0, 70.0, 0, 210.0, boundary_type='vacuum')
graphite_cell = openmc.Cell(fill=graphite, region=-graphite_surface)

# universe with graphite cell - i think this is right
root_universe = openmc.Universe(cells=[graphite_cell])

# geometry - i think this is right
geometry = openmc.Geometry(root_universe)

#source (AmBe) - confident on this part
source = openmc.IndependentSource()
source.space = openmc.stats.Point((0,0,10))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([4.2e6], [1]) #average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually


# settings - this works but not sure if its ideal
settings = openmc.Settings()
settings.source = source
settings.particles = 10000
settings.batches = 20
settings.inactive = 5
settings.output= {'tallies': True}

# tally
tally = openmc.Tally()
tally.score = ['flux']
mat_filter = openmc.MaterialFilter(graphite)
tally.filters = openmc.Filters = [mat_filter]
tallies = openmc.Tallies([tally])
tally.scores =['flux']

# export to xml
geometry.export_to_xml()
materials.export_to_xml()
settings.export_to_xml()
tallies.export_to_xml()

# run - is this what i want or do I want something else if i'm just trying to get neutron flux/counts
openmc.run()

and thanks for the help so far :smiley:

4

Oh I was missing

settings.run_mode = ‘fixed source’

seems to be running now.

Thank you for the help with the tallies section but it seems my tallies.xml file is completely empty.

Heres the current version of the code

import openmc

# create graphite - 100% carbon, density 2.26g/cm^3
graphite = openmc.Material()
graphite.add_element('C', 1.0)
graphite.set_density('g/cm3', 2.26)
materials = openmc.Materials([graphite])

# cuboid - simple graphite cube for now
graphite_surface = openmc.model.RectangularParallelepiped(-70.0, 70.0, -70.0, 70.0, 0, 210.0, boundary_type='vacuum')
graphite_cell = openmc.Cell(fill=graphite, region=-graphite_surface)

# universe with graphite cell - i think this is right
root_universe = openmc.Universe(cells=[graphite_cell])

# geometry - i think this is right
geometry = openmc.Geometry(root_universe)

#source (AmBe) - confident on this part
source = openmc.IndependentSource()
source.space = openmc.stats.Point((0,0,10))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([4.2e6], [1]) #average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually

# settings - this works but not sure if its ideal
settings = openmc.Settings()
settings.source = source
settings.run_mode = 'fixed source'
settings.particles = 10000
settings.batches = 20
settings.inactive = 5
settings.output= {'tallies': True}

# tally
tally = openmc.Tally()
tally.score = ['flux']
mat_filter = openmc.MaterialFilter(graphite)
tally.filters = openmc.Filters = [mat_filter]
tallies = openmc.Tallies([tally])
tally.scores =['flux']

#plot
plots = openmc.Plots()

# export to xml
geometry.export_to_xml()
materials.export_to_xml()
settings.export_to_xml()
tallies.export_to_xml()
plots.export_to_xml()

# run - is this what i want or do I want something else if i'm just trying to get neutron flux/counts
openmc.run()


#post process
sp = openmc.StatePoint('statepoint.20.h5')

flux = tally.get_slice(scores=['flux'])
print(flux)
5

I’ve got the basics of the code down, I’m still struggling to understand why my tallies are completely empty. I’ve added filters now and some other things to clean it up a bit

import openmc
import numpy as np

# create graphite - 100% carbon, density 2.26g/cm^3
graphite = openmc.Material()
graphite.add_element('C', 1.0)
graphite.set_density('g/cm3', 2.26)
materials = openmc.Materials([graphite])

# cuboid - simple graphite cube for now
graphite_surface = openmc.model.RectangularParallelepiped(-70.0, 70.0, -70.0, 70.0, 0, 210.0, boundary_type='vacuum')
graphite_cell = openmc.Cell(fill=graphite, region=-graphite_surface)

# geometry 
geometry = openmc.Geometry([graphite_cell])

#source (AmBe) - confident on this part
source = openmc.IndependentSource()
source.space = openmc.stats.Point((0,0,10))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([4.2e6], [1]) #average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually

# settings - this works but not sure if its ideal
settings = openmc.Settings()
settings.source = source
settings.run_mode = 'fixed source'
settings.particles = 10000
settings.batches = 10
settings.inactive = 5
settings.output = {'tallies': True}

# filters
neutron_filter = openmc.ParticleFilter(['neutron']) #selects neutrons
energy_filter = openmc.EnergyFilter(np.linspace(0, 0.025, 200))  #select thermal neutrons
cell_filter = openmc.CellFilter([graphite_cell]) #select cell

# tally
tally = openmc.Tally(name='FluxTally')
tally.scores =['flux']
tally.filters = [cell_filter, neutron_filter, energy_filter]
tallies = openmc.Tallies([tally])

# export to xml
#geometry.export_to_xml()
#materials.export_to_xml()
#settings.export_to_xml()
#tallies.export_to_xml()
#plots.export_to_xml()

# run - is this what i want or do I want something else if i'm just trying to get neutron flux/counts
model = openmc.model.Model(geometry, materials, settings, tallies)

results_filename = model.run(tracks=True)


#post process
results = openmc.StatePoint(results_filename)

cell_tally = results.get_tally(name='FluxTally')

Any help on why my tallies are completely empty would be great, thanks!

6

Are you accessing the cell_tally.mean attribute

7

I have a tally with results now, however I want to be able to plot a flux “heat map” of a x-z or y-z slice with my data, do you know how I could go about doing that?

iimport openmc
import numpy as np
import matplotlib.pyplot as plt
import os

# create graphite - 100% carbon, density 2.26g/cm^3
graphite = openmc.Material()
graphite.add_element('C', 1.0)
graphite.set_density('g/cm3', 2.26)
materials = openmc.Materials([graphite])

# cuboid - simple graphite cube for now
graphite_surface = openmc.model.RectangularParallelepiped(-70.0, 70.0, -70.0, 70.0, 0, 210.0, boundary_type='vacuum')
graphite_cell = openmc.Cell(fill=graphite, region=-graphite_surface)

# geometry 
geometry = openmc.Geometry([graphite_cell])

# source (AmBe)
source = openmc.IndependentSource()
source.space = openmc.stats.Point((0, 0, 10))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([4.2e6], [1])  # average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually

# settings 
settings = openmc.Settings()
settings.source = source
settings.run_mode = 'fixed source'
settings.particles = 1000
settings.batches = 10
settings.inactive = 0
settings.output = {'tallies': True}

# filters
neutron_filter = openmc.ParticleFilter(['neutron'])  # selects neutrons
energy_filter = openmc.EnergyFilter(np.linspace(0, 0.025, 201))  # select thermal neutrons
energy_filter2 = openmc.EnergyFilter(np.linspace(0, 4.2e6, 201))  # selects all
cell_filter = openmc.CellFilter([graphite_cell])  # select cell

# mesh filter
mesh = openmc.RegularMesh()
mesh.dimension = [100, 100]
mesh.lower_left = [-70.0, 0]
mesh.upper_right = [70.0, 210.0]
mesh_filter = openmc.MeshFilter(mesh)
mesh_tally = openmc.Tally(name='tallies_on_mesh')
mesh_tally.filters = [mesh_filter, energy_filter2, neutron_filter, cell_filter]
mesh_tally.scores = ['flux']
tallies = openmc.Tallies([mesh_tally])

if os.path.exists('tallies.h5'):
    os.remove('tallies.h5')
if os.path.exists('statepoint.10.h5'):
    os.remove('statepoint.10.h5')
if os.path.exists('summary.h5'):
    os.remove('summary.h5')
if os.path.exists('tallies.out'):
    os.remove('tallies.out')
if os.path.exists('tracks.h5'):
    os.remove('tracks.h5')


# run
model = openmc.model.Model(geometry, materials, settings, tallies)
results_filename = model.run()

# post process - not sure
results = openmc.StatePoint(results_filename)

# Get the tally results
my_tally = results.get_tally(scores=['flux'])
my_slice = my_tally.get_slice(scores=['flux'])
print(my_tally)
# Get the mean of the tally results
mean = my_slice.mean.shape = (mesh.dimension[0], mesh.dimension[1])

# Plot the reshaped tally results
plt.imshow(mean, extent=(mesh.lower_left[0], mesh.upper_right[0], mesh.lower_left[1], mesh.upper_right[1]))
plt.show()

This is what i’ve tried so far but the output is far from ideal. I think the main issue is that my output is one dimensional and I’m not sure what is going wrong with my mesh to make that happen.

8

I have written a small package for plotting regular mesh tallies that might be of interest

The plotter supports plotting mesh tallies

9

Hi, thank you for the packages. I used your first one and I have what seems to be a good slice of data.
graphite_flux
graphite_flux_log

Is there anyway I could extract the diffusion length from the data produced now. I’ve looked up some sort but it was using multigroup diffusion coefficient but I’m not sure that would be appropiate.

Once again really appreciate the help

10

L^2 should equal 1/6-th of the mean distance between birth and absorption of a thermal neutron. See https://www.nuclear-power.com/nuclear-power/reactor-physics/neutron-diffusion-theory/diffusion-length/

since L=53.5 cm for graphite, i dont think your geometry is big enough to sample the mean ‘life’ distance

11

I am happy to vary the dimension of the graphite, how do I go about sampling the lifetime of a thermal neutron?

12

Hi again,

Just another question regarding mesh tallies, right now my mesh filter looks like this

# mesh filter
    mesh = openmc.RegularMesh()
    mesh.dimension = [140, 1, 210]
    mesh.lower_left = [-70.0, -70.0, 0]
    mesh.upper_right = [70.0, 70.0, 210.0]
    mesh_filter = openmc.MeshFilter(mesh)
    mesh_tally = openmc.Tally(name='mesh_tally')
    mesh_tally.filters = [mesh_filter, neutron_filter, cell_filter]
    mesh_tally.scores = ['flux']
    tallies = openmc.Tallies([mesh_tally])

Is there a way for the mesh tally to only tally thermal neutrons, I’ve tried adding an energy filter but I get ValueErrors when trying to plot.

Thanks

13

I would have gobe with an energy filter with a couple of low energy values like you mentioned. Not sure why this would cause an error.

14

Oh a value error with my plotter, oh right yes i can believe that. Can you poster the value error and i shall fix it

15 
ValueError: cannot select an axis to squeeze out which has size not equal to one

I’ve tried plotting through another means but it looks quite noisey.
graphite_flux_log_thermal

Had to do some weird resizing with my arrays to get it to work though

    my_slice = my_tally.get_slice(scores=['flux'])
    print(my_slice.mean)
    print(my_slice.mean.shape)
    new_dim1 = 140
    new_dim2 = int(my_slice.mean.shape[0] / new_dim1)
    reshaped = my_slice.mean.reshape(new_dim1, new_dim2)
    plt.imshow(reshaped, extent=mesh.bounding_box.extent['xz'], norm=LogNorm())
    plt.savefig(f'{material_choice}_flux_log_thermal.png')
    plt.show()

I think this is probably due to do with my energy filters number of bins, not sure what kind of number I should use.

    energy_filter_thermal = openmc.EnergyFilter(np.linspace(0, 0.625, 200))  # select thermal neutrons
    energy_filter_all = openmc.EnergyFilter(np.linspace(0, 4.2e6, 200))  # selects all
16

Try an energy function filter with just one energy bin

instead of this

energy_filter_all = openmc.EnergyFilter(np.linspace(0, 4.2e6, 200))

try this

energy_filter_all = openmc.EnergyFilter([0, 4.2e6])
17

Even with one energy bin I still get

ValueError: cannot select an axis to squeeze out which has size not equal to one
18

I was able to reproduce that error ValueError: cannot select an axis to squeeze out which has size not equal to one

I’ve put a check in the openmc_regular_mesh_plotter with this PR Adding energy filter num_dims check to avoid crash by shimwell · Pull Request #64 · fusion-energy/openmc_regular_mesh_plotter · GitHub

I’ve added tests that check the valueerror gets raised for EnergyFilters with num dim above 1

I’ve also released a new version of openmc-regular-mesh-plotter that has the check in it.
pip install openmc-regular-mesh-plotter==1.2.1

19

I still seem to be getting an error even after updating of

ValueError: cannot select an axis to squeeze out which has size not equal to one

This is my code excluding the materials section if it helps

    # source (AmBe)
    source = openmc.IndependentSource()
    source.space = openmc.stats.Point((0, 0, 20.32))
    source.angle = openmc.stats.Isotropic()
    source.energy = openmc.stats.Discrete([4.2e6], [1])  # average energy of neutrons from AmBe source is 4.2MeV, want to have a spectrum eventually

    # settings
    settings = openmc.Settings()
    settings.source = source
    settings.run_mode = 'fixed source'
    settings.particles = 10000
    settings.batches = 20
    settings.inactive = 0
    settings.output = {'tallies': True}

    # filters
    neutron_filter = openmc.ParticleFilter(['neutron'])  # selects neutrons
    energy_filter_thermal = openmc.EnergyFilter(np.linspace(0, 0.030, 1))  # select thermal neutrons
    energy_filter_all = openmc.EnergyFilter(np.linspace(0, 4.2e6))  # selects all - still doesn't work?
    # using energy filters breaks it, not sure why

    # mesh filter

    mesh = openmc.RegularMesh()
    mesh.dimension = [dim_xy[1]*2, 1, dim_z[1]]
    mesh.lower_left = [dim_xy[0], dim_xy[0], dim_z[0]]
    mesh.upper_right = [dim_xy[1], dim_xy[1], dim_z[1]]
    mesh_filter = openmc.MeshFilter(mesh)
    mesh_tally = openmc.Tally(name='mesh_tally')
    mesh_tally.filters = [mesh_filter, cell_filter, energy_filter_thermal]
    mesh_tally.scores = ['flux']
    tallies = openmc.Tallies([mesh_tally])

    # remove old files
    files_to_remove = ['tallies.h5', f'statepoint.{settings.batches}.h5', 'summary.h5', 'tallies.out', 'tracks.h5', 'geometry.xml',
                       'materials.xml', 'settings.xml', 'tallies.xml', 'plots.xml']
    for file in files_to_remove:
        if os.path.exists(file):
            os.remove(file)

    # run
    model = openmc.model.Model(geometry, materials, settings, tallies)
    results_filename = model.run(tracks=True)

    # post process - not sure
    results = openmc.StatePoint(results_filename)

    # Get the tally results
    my_tally = results.get_tally(name='mesh_tally')

    # plotting with matplotlib
    # my_slice = my_tally.get_slice(scores=['flux'])
    # print(my_slice.mean)
    # print(my_slice.mean.shape)
    # new_dim1 = (dim_xy[1]*2)
    # new_dim2 = int(my_slice.mean.shape[0] / new_dim1)
    # reshaped = my_slice.mean.reshape(new_dim1, new_dim2)
    # plt.imshow(reshaped, extent=mesh.bounding_box.extent['xz'], norm=LogNorm())
    # plt.xlabel('x [cm]')
    # plt.ylabel('z [cm]')
    # plt.title(f'Flux in {material_choice} (log scale)')
    # plt.savefig(f'{material_choice}_flux_log_thermal.png')
    # plt.show()

    # plotting with openmc_regular_mesh_plotter
    plot = plot_mesh_tally(basis='xz',
                           tally=my_tally,
                           outline=True,
                           geometry=geometry,
                           norm=LogNorm(),
                           colorbar=True)

    plot.title.set_text(f'Flux in {material_choice} (log scale)')
    plot.figure.savefig(f'{material_choice}_flux_log.png')

    plot = plot_mesh_tally(basis='xz',
                           tally=my_tally,
                           outline=True,
                           geometry=geometry,
                           colorbar=True)

    plot.title.set_text(f'Flux in {material_choice}')
    plot.figure.savefig(f'{material_choice}_flux.png')
20

Not sure but I wonder if you need a 2nd value in that EnergyFilter, as it is only one value long at the moment. How about putting a 2 at the end

openmc.EnergyFilter(np.linspace(0, 0.030, 2)