Hi all,
I’m trying to adapt the code here (“https://github.com/openmc-dev/openmc-notebooks/blob/main/gamma-detector.ipynb”) to my own study, but I’m getting this error. What could be the reason? I’d appreciate it if you could help me.
Thank you
The error message says that the code is unable to find the statepoint.10.h5 file.
This file is generated by the simulation as an output of the results.
The code then tries to read in the file, but as it can not be found it can not proceed.
If the simulation ran then it will tell you the name of the statepoint output file that was saved, I’ve highlighted it in blue in this screen shot. If you have changed the number of batches then the filename can be different.
Try match the filename of the statepoint printed to the terminal with the statepoint that you are loading.
You might have to change the line of code that looks like this sp = openmc.StatePoint('statepoint.10.h5') so that the filenames are the same
Actually, I wrote it in the code like in the example.
sp = openmc.StatePoint(‘statepoint.10.h5’)
tally = sp.get_tally(name=‘pulse-height’)
pulse_height_values = tally.get_values(scores=[‘pulse-height’]).flatten()
energy_bin_centers = energy_bins[1:] + 0.5 * (energy_bins[1] - energy_bins[0])
Did you run all the cells above in the jupyter notebook. There is a cell above that has an openmc.run() command in an this should produce the statepoint file.
I see the example has a file path so you might want to change
openmc.run("/home/cpf/Desktop/openmc/build/bin/openmc")
to
openmc.run()
I run the this code
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import openmc
from scipy.constants import Avogadro
cesium_iodide = openmc.Material()
cesium_iodide.add_element(‘Cs’, 0.5)
cesium_iodide.add_element(‘I’, 0.5)
cesium_iodide.set_density(‘g/cm3’, 4.51)
air = openmc.Material(name=‘air’)
air.add_element(‘N’, 0.784431)
air.add_element(‘O’, 0.210748)
air.add_element(‘Ar’,0.0046)
air.temperature = 300
air.set_density(‘g/cm3’, 0.0012)
materials = openmc.Materials([cesium_iodide,air])
materials.cross_section=‘C:/Users/Gokcenur/endfb71_hdf5/cross_sections.xml’
materials.export_to_xml()
##
#Create Geometry #
##
Create the detector
-------------------
Create the dectector
--------------------
det_rect_prism = openmc.model.rectangular_prism(1.3, 1.3, origin=(0.0, 0.0))
det_zmin = openmc.ZPlane(z0=-1.0, boundary_type=‘transmission’)
det_zmax = openmc.ZPlane(z0=+1.0, boundary_type=‘transmission’)
detector = det_rect_prism & +det_zmin & -det_zmax
detector_cell = openmc.Cell(fill= cesium_iodide, region=detector)
Create a source sphere just to aid vitualisation but it can be removed later
----------------------------------------------------------------------------
source = openmc.Sphere(x0=+5.0, y0=0.0, z0=0.0, r=2.0, name=‘Source Container’)
source_cell = openmc.Cell(fill=None, region=-source)
Create the external environment
--------------------------------
sphere = openmc.Sphere(x0=30.0, y0=0.0, z0=0.0, r = 80, boundary_type = ‘vacuum’)
environment = -sphere & ~detector & +source
environment_cell = openmc.Cell(fill=air, region=environment)
universe = openmc.Universe(cells=[environment_cell, detector_cell, source_cell])
universe.plot(width = (20.0,20.0),basis = ‘xz’)
geometry = openmc.Geometry(universe)
geometry.export_to_xml()
##
#Plot the Geometry using slice (2D) plot #
##
p = openmc.Plot()
p.filename = ‘detector’
p.width = (250, 250)
p.pixels = (600, 600)
p.color_by = ‘cell’
p.colors = {environment_cell: ‘yellow’, detector_cell: ‘blue’, source_cell: ‘red’}
Create a plot and export to XML
plots = openmc.Plots([p])
plots.export_to_xml()
openmc.plot_geometry()
#Defining Sources
source = openmc.Source()
source.space = openmc.stats.Point(xyz=(5, 0, 0))
source.angle = openmc.stats.Isotropic()
This is a Co60 source, see the task on sources to understand it
source.energy = openmc.stats.Discrete([0.661e6], [1])
source.particle = ‘photon’
source.strength = 1
#Defining Settings
settings = openmc.Settings()
settings.particles = 100000
settings.batches = 20
settings.photon_transport = True
settings.source = [source]
settings.run_mode = ‘fixed source’
settings.export_to_xml()
#Defining Tallies
tallies = openmc.Tallies()
energy_bins = np.linspace(0, 1e6, 1001)
energy_filter = openmc.EnergyFilter(energy_bins)
cell_filter = openmc.CellFilter(detector_cell)
tally = openmc.Tally(name=‘pulse-height’)
tally.filters = [cell_filter, energy_filter]
tally.scores = [‘pulse-height’]
tallies.append(tally)
tallies.export_to_xml()
#Running the Simulation
openmc.run()
#PlotPulseHeightTally
sp = openmc.StatePoint(‘statepoint.10.h5’)
tally = sp.get_tally(name=‘pulse-height’)
pulse_height_values = tally.get_values(scores=[‘pulse-height’]).flatten()
energy_bin_centers = energy_bins[1:] + 0.5 * (energy_bins[1] - energy_bins[0])
plt.figure()
plt.semilogy(energy_bin_centers, pulse_height_values)
plot the strongest sources as vertical lines
plt.axvline(x=800_000, color=“red”, ls=“:”) # source_1
plt.axvline(x=661_700, color=“red”, ls=“:”) # source_2
plt.xlabel(‘Energy [eV]’)
plt.ylabel(‘Counts’)
plt.title(‘Pulse Height Values’)
plt.grid(True)
plt.tight_layout()
#Post Processing Results
a, b, c = 1000, 4, 0.0002
number_broadening_samples = 1e6
def gauss(E, a=a, b=b, c=c):
sigma = (a + b * (E + c * E**2)**0.5) / (2 * (2 * np.log(2))**0.5)
return np.random.normal(loc=E, scale=sigma)
samples = np.random.choice(energy_bin_centers[1:], size=int(number_broadening_samples), p=pulse_height_values[1:]/np.sum(pulse_height_values[1:]))
broaded_pulse_height_values = gauss(samples)
broaded_spectrum, _ = np.histogram(broaded_pulse_height_values, bins=energy_bins)
renormalized_broaded_spectrum = broaded_spectrum / np.sum(broaded_spectrum) * np.sum(pulse_height_values[1:])
plt.figure()
plt.semilogy(energy_bin_centers[1:], pulse_height_values[1:], label=“original simulation result”)
plt.semilogy(energy_bin_centers[1:], renormalized_broaded_spectrum[1:], label=“broaded detector response”)
plot the strongest sources as vertical lines
plt.axvline(x=800_000, color=“red”, ls=“:”, label=“gamma source”) # source_1
plt.axvline(x=661_700, color=“red”, ls=“:”) # source_2
plt.legend()
plt.xlabel(‘Energy [eV]’)
plt.ylabel(‘Counts’)
plt.title(‘Pulse Height Values’)
plt.grid(True)
plt.tight_layout()
Looks like you have got 20 batches there instead of 10, so this will change the statepoint filename
1 Like
Thanks, I didn’t pay attention there.