This notebook demonstrates some basic post-processing tasks that can be performed with the Python API, such as plotting a 2D mesh tally and plotting neutron source sites from an eigenvalue calculation. The problem we will use is a simple reflected pin-cell.
%matplotlib inline
from IPython.display import Image
import numpy as np
import matplotlib.pyplot as plt
import openmc
First we need to define materials that will be used in the problem. We'll create three materials for the fuel, water, and cladding of the fuel pin.
# 1.6 enriched fuel
fuel = openmc.Material(name='1.6% Fuel')
fuel.set_density('g/cm3', 10.31341)
fuel.add_nuclide('U235', 3.7503e-4)
fuel.add_nuclide('U238', 2.2625e-2)
fuel.add_nuclide('O16', 4.6007e-2)
# borated water
water = openmc.Material(name='Borated Water')
water.set_density('g/cm3', 0.740582)
water.add_nuclide('H1', 4.9457e-2)
water.add_nuclide('O16', 2.4732e-2)
water.add_nuclide('B10', 8.0042e-6)
# zircaloy
zircaloy = openmc.Material(name='Zircaloy')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_nuclide('Zr90', 7.2758e-3)
With our three materials, we can now create a materials file object that can be exported to an actual XML file.
# Instantiate a Materials collection
materials = openmc.Materials([fuel, water, zircaloy])
# Export to "materials.xml"
materials.export_to_xml()
Now let's move on to the geometry. Our problem will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces -- in this case two cylinders and six reflective planes.
# Create cylinders for the fuel and clad
fuel_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.39218)
clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.45720)
# Create boundary planes to surround the geometry
min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective')
max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')
min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')
max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective')
min_z = openmc.ZPlane(z0=-0.63, boundary_type='reflective')
max_z = openmc.ZPlane(z0=+0.63, boundary_type='reflective')
With the surfaces defined, we can now create cells that are defined by intersections of half-spaces created by the surfaces.
# Create a Universe to encapsulate a fuel pin
pin_cell_universe = openmc.Universe(name='1.6% Fuel Pin')
# Create fuel Cell
fuel_cell = openmc.Cell(name='1.6% Fuel')
fuel_cell.fill = fuel
fuel_cell.region = -fuel_outer_radius
pin_cell_universe.add_cell(fuel_cell)
# Create a clad Cell
clad_cell = openmc.Cell(name='1.6% Clad')
clad_cell.fill = zircaloy
clad_cell.region = +fuel_outer_radius & -clad_outer_radius
pin_cell_universe.add_cell(clad_cell)
# Create a moderator Cell
moderator_cell = openmc.Cell(name='1.6% Moderator')
moderator_cell.fill = water
moderator_cell.region = +clad_outer_radius
pin_cell_universe.add_cell(moderator_cell)
OpenMC requires that there is a "root" universe. Let us create a root cell that is filled by the pin cell universe and then assign it to the root universe.
# Create root Cell
root_cell = openmc.Cell(name='root cell')
root_cell.fill = pin_cell_universe
# Add boundary planes
root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z
# Create root Universe
root_universe = openmc.Universe(universe_id=0, name='root universe')
root_universe.add_cell(root_cell)
We now must create a geometry that is assigned a root universe, put the geometry into a geometry file, and export it to XML.
# Create Geometry and set root Universe
geometry = openmc.Geometry(root_universe)
# Export to "geometry.xml"
geometry.export_to_xml()
With the geometry and materials finished, we now just need to define simulation parameters. In this case, we will use 10 inactive batches and 90 active batches each with 5000 particles.
# OpenMC simulation parameters
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 5000
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)
# Export to "settings.xml"
settings.export_to_xml()
Let us also create a plot file that we can use to verify that our pin cell geometry was created successfully.
plot = openmc.Plot.from_geometry(geometry)
plot.pixels = (250, 250)
plot.to_ipython_image()
As we can see from the plot, we have a nice pin cell with fuel, cladding, and water! Before we run our simulation, we need to tell the code what we want to tally. The following code shows how to create a 2D mesh tally.
# Instantiate an empty Tallies object
tallies = openmc.Tallies()
# Create mesh which will be used for tally
mesh = openmc.RegularMesh()
mesh.dimension = [100, 100]
mesh.lower_left = [-0.63, -0.63]
mesh.upper_right = [0.63, 0.63]
# Create mesh filter for tally
mesh_filter = openmc.MeshFilter(mesh)
# Create mesh tally to score flux and fission rate
tally = openmc.Tally(name='flux')
tally.filters = [mesh_filter]
tally.scores = ['flux', 'fission']
tallies.append(tally)
# Export to "tallies.xml"
tallies.export_to_xml()
Now we a have a complete set of inputs, so we can go ahead and run our simulation.
# Run OpenMC!
openmc.run()
%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% ############### %%%%%%%%%%%%%%%%%%%%%%%% ################## %%%%%%%%%%%%%%%%%%%%%%% ################### %%%%%%%%%%%%%%%%%%%%%%% #################### %%%%%%%%%%%%%%%%%%%%%% ##################### %%%%%%%%%%%%%%%%%%%%% ###################### %%%%%%%%%%%%%%%%%%%% ####################### %%%%%%%%%%%%%%%%%% ####################### %%%%%%%%%%%%%%%%% ###################### %%%%%%%%%%%%%%%%% #################### %%%%%%%%%%%%%%%%% ################# %%%%%%%%%%%%%%%%% ############### %%%%%%%%%%%%%%%% ############ %%%%%%%%%%%%%%% ######## %%%%%%%%%%%%%% %%%%%%%%%%% | The OpenMC Monte Carlo Code Copyright | 2011-2021 MIT and OpenMC contributors License | https://docs.openmc.org/en/latest/license.html Version | 0.13.0-dev Git SHA1 | 3dd81a1316ac3b5a0633e4b7a290f3bc97a066d9 Date/Time | 2021-06-26 15:28:06 OpenMP Threads | 2 Reading settings XML file... Reading cross sections XML file... Reading materials XML file... Reading geometry XML file... Reading U235 from /home/shriwise/opt/openmc/xs/nndc_hdf5/U235.h5 Reading U238 from /home/shriwise/opt/openmc/xs/nndc_hdf5/U238.h5 Reading O16 from /home/shriwise/opt/openmc/xs/nndc_hdf5/O16.h5 Reading H1 from /home/shriwise/opt/openmc/xs/nndc_hdf5/H1.h5 Reading B10 from /home/shriwise/opt/openmc/xs/nndc_hdf5/B10.h5 Reading Zr90 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr90.h5 Minimum neutron data temperature: 294.0 K Maximum neutron data temperature: 294.0 K Reading tallies XML file... Preparing distributed cell instances... Writing summary.h5 file... Maximum neutron transport energy: 20000000.0 eV for U235 Initializing source particles... ====================> K EIGENVALUE SIMULATION <==================== Bat./Gen. k Average k ========= ======== ==================== 1/1 1.05252 2/1 1.03787 3/1 1.01943 4/1 1.03989 5/1 1.06679 6/1 1.03713 7/1 1.02400 8/1 1.04289 9/1 1.05130 10/1 1.00878 11/1 1.06773 12/1 1.03922 1.05347 +/- 0.01426 13/1 1.05156 1.05283 +/- 0.00826 14/1 1.06049 1.05475 +/- 0.00614 15/1 1.01018 1.04583 +/- 0.01010 16/1 1.04020 1.04490 +/- 0.00830 17/1 1.05579 1.04645 +/- 0.00719 18/1 1.01592 1.04264 +/- 0.00730 19/1 1.06881 1.04554 +/- 0.00707 20/1 1.02985 1.04397 +/- 0.00651 21/1 1.01496 1.04134 +/- 0.00645 22/1 1.05330 1.04233 +/- 0.00598 23/1 1.05170 1.04305 +/- 0.00554 24/1 1.02888 1.04204 +/- 0.00523 25/1 1.04083 1.04196 +/- 0.00487 26/1 1.01235 1.04011 +/- 0.00492 27/1 1.02785 1.03939 +/- 0.00468 28/1 1.04556 1.03973 +/- 0.00442 29/1 1.05400 1.04048 +/- 0.00425 30/1 1.06213 1.04157 +/- 0.00417 31/1 0.99934 1.03955 +/- 0.00445 32/1 1.04433 1.03977 +/- 0.00425 33/1 1.05184 1.04030 +/- 0.00409 34/1 1.03971 1.04027 +/- 0.00392 35/1 1.05272 1.04077 +/- 0.00379 36/1 1.06881 1.04185 +/- 0.00380 37/1 1.03344 1.04154 +/- 0.00367 38/1 1.04726 1.04174 +/- 0.00354 39/1 1.01440 1.04080 +/- 0.00354 40/1 1.03534 1.04062 +/- 0.00343 41/1 1.04429 1.04073 +/- 0.00332 42/1 1.02142 1.04013 +/- 0.00327 43/1 1.03895 1.04010 +/- 0.00317 44/1 1.05985 1.04068 +/- 0.00313 45/1 1.04737 1.04087 +/- 0.00304 46/1 1.04796 1.04106 +/- 0.00297 47/1 1.06708 1.04177 +/- 0.00297 48/1 1.06523 1.04238 +/- 0.00295 49/1 0.99626 1.04120 +/- 0.00311 50/1 1.04077 1.04119 +/- 0.00303 51/1 1.06327 1.04173 +/- 0.00301 52/1 1.06508 1.04229 +/- 0.00299 53/1 1.03689 1.04216 +/- 0.00292 54/1 1.02899 1.04186 +/- 0.00287 55/1 1.03267 1.04166 +/- 0.00281 56/1 1.05790 1.04201 +/- 0.00277 57/1 1.04353 1.04204 +/- 0.00271 58/1 1.04657 1.04214 +/- 0.00266 59/1 1.02914 1.04187 +/- 0.00261 60/1 1.04882 1.04201 +/- 0.00257 61/1 1.01905 1.04156 +/- 0.00255 62/1 1.03995 1.04153 +/- 0.00251 63/1 1.05377 1.04176 +/- 0.00247 64/1 1.02909 1.04153 +/- 0.00243 65/1 1.06892 1.04202 +/- 0.00244 66/1 1.04216 1.04203 +/- 0.00240 67/1 1.03473 1.04190 +/- 0.00236 68/1 1.04114 1.04188 +/- 0.00232 69/1 1.04955 1.04201 +/- 0.00228 70/1 1.05464 1.04222 +/- 0.00225 71/1 1.02859 1.04200 +/- 0.00223 72/1 1.05387 1.04219 +/- 0.00220 73/1 1.05039 1.04232 +/- 0.00217 74/1 1.04338 1.04234 +/- 0.00213 75/1 1.05838 1.04259 +/- 0.00211 76/1 1.03831 1.04252 +/- 0.00208 77/1 1.03555 1.04242 +/- 0.00205 78/1 1.05684 1.04263 +/- 0.00204 79/1 1.04267 1.04263 +/- 0.00201 80/1 1.05813 1.04285 +/- 0.00199 81/1 1.03512 1.04274 +/- 0.00196 82/1 1.07081 1.04313 +/- 0.00198 83/1 1.04476 1.04315 +/- 0.00195 84/1 1.05153 1.04327 +/- 0.00192 85/1 1.03939 1.04322 +/- 0.00190 86/1 1.04218 1.04320 +/- 0.00187 87/1 1.03688 1.04312 +/- 0.00185 88/1 1.03480 1.04301 +/- 0.00183 89/1 1.05089 1.04311 +/- 0.00181 90/1 1.06251 1.04336 +/- 0.00180 91/1 1.04054 1.04332 +/- 0.00178 92/1 1.05340 1.04344 +/- 0.00176 93/1 1.05938 1.04364 +/- 0.00175 94/1 1.02741 1.04344 +/- 0.00174 95/1 1.08249 1.04390 +/- 0.00178 96/1 1.02858 1.04372 +/- 0.00177 97/1 1.03983 1.04368 +/- 0.00175 98/1 1.04715 1.04372 +/- 0.00173 99/1 1.07443 1.04406 +/- 0.00175 100/1 1.04461 1.04407 +/- 0.00173 Creating state point statepoint.100.h5... =======================> TIMING STATISTICS <======================= Total time for initialization = 2.4568e-01 seconds Reading cross sections = 2.3233e-01 seconds Total time in simulation = 1.1761e+02 seconds Time in transport only = 1.1757e+02 seconds Time in inactive batches = 2.0641e+00 seconds Time in active batches = 1.1554e+02 seconds Time synchronizing fission bank = 2.1808e-02 seconds Sampling source sites = 1.8421e-02 seconds SEND/RECV source sites = 3.3183e-03 seconds Time accumulating tallies = 2.6283e-03 seconds Time writing statepoints = 4.2804e-03 seconds Total time for finalization = 2.2731e-02 seconds Total time elapsed = 1.1789e+02 seconds Calculation Rate (inactive) = 24223.6 particles/second Calculation Rate (active) = 3894.66 particles/second ============================> RESULTS <============================ k-effective (Collision) = 1.04491 +/- 0.00149 k-effective (Track-length) = 1.04407 +/- 0.00173 k-effective (Absorption) = 1.04203 +/- 0.00169 Combined k-effective = 1.04355 +/- 0.00131 Leakage Fraction = 0.00000 +/- 0.00000
Our simulation ran successfully and created a statepoint file with all the tally data in it. We begin our analysis here loading the statepoint file and 'reading' the results. By default, data from the statepoint file is only read into memory when it is requested. This helps keep the memory use to a minimum even when a statepoint file may be huge.
# Load the statepoint file
sp = openmc.StatePoint('statepoint.100.h5')
Next we need to get the tally, which can be done with the StatePoint.get_tally(...)
method.
tally = sp.get_tally(scores=['flux'])
print(tally)
Tally ID = 1 Name = flux Filters = MeshFilter Nuclides = total Scores = ['flux', 'fission'] Estimator = tracklength
The statepoint file actually stores the sum and sum-of-squares for each tally bin from which the mean and variance can be calculated as described here. The sum and sum-of-squares can be accessed using the sum
and sum_sq
properties:
tally.sum
array([[[0.41279586, 0. ]], [[0.41176924, 0. ]], [[0.41096843, 0. ]], ..., [[0.4095409 , 0. ]], [[0.40836217, 0. ]], [[0.40852022, 0. ]]])
However, the mean and standard deviation of the mean are usually what you are more interested in. The Tally class also has properties mean
and std_dev
which automatically calculate these statistics on-the-fly.
print(tally.mean.shape)
(tally.mean, tally.std_dev)
(10000, 1, 2)
(array([[[0.00458662, 0. ]], [[0.00457521, 0. ]], [[0.00456632, 0. ]], ..., [[0.00455045, 0. ]], [[0.00453736, 0. ]], [[0.00453911, 0. ]]]), array([[[1.74741992e-05, 0.00000000e+00]], [[1.68457472e-05, 0.00000000e+00]], [[1.75888801e-05, 0.00000000e+00]], ..., [[1.79971274e-05, 0.00000000e+00]], [[1.89308740e-05, 0.00000000e+00]], [[1.75231302e-05, 0.00000000e+00]]]))
The tally data has three dimensions: one for filter combinations, one for nuclides, and one for scores. We see that there are 10000 filter combinations (corresponding to the 100 x 100 mesh bins), a single nuclide (since none was specified), and two scores. If we only want to look at a single score, we can use the get_slice(...)
method as follows.
flux = tally.get_slice(scores=['flux'])
fission = tally.get_slice(scores=['fission'])
print(flux)
Tally ID = 2 Name = flux Filters = MeshFilter Nuclides = total Scores = ['flux'] Estimator = tracklength
To get the bins into a form that we can plot, we can simply change the shape of the array since it is a numpy array.
flux.std_dev.shape = (100, 100)
flux.mean.shape = (100, 100)
fission.std_dev.shape = (100, 100)
fission.mean.shape = (100, 100)
fig = plt.subplot(121)
fig.imshow(flux.mean)
fig2 = plt.subplot(122)
fig2.imshow(fission.mean)
<matplotlib.image.AxesImage at 0x7f3f6e467e10>
Now let's say we want to look at the distribution of relative errors of our tally bins for flux. First we create a new variable called relative_error
and set it to the ratio of the standard deviation and the mean, being careful not to divide by zero in case some bins were never scored to.
# Determine relative error
relative_error = np.zeros_like(flux.std_dev)
nonzero = flux.mean > 0
relative_error[nonzero] = flux.std_dev[nonzero] / flux.mean[nonzero]
# distribution of relative errors
ret = plt.hist(relative_error[nonzero], bins=50)
Source sites can be accessed from the source
property. As shown below, the source sites are represented as a numpy array with a structured datatype.
sp.source
array([((0.20665803, 0.15081559, -0.57355059), ( 0.49473673, 0.67921184, -0.54213177), 2077978.15846043, 1., 0, 0, 0), ((0.02302023, -0.02944101, -0.45025678), ( 0.53648981, 0.51827967, 0.66600666), 206149.19886773, 1., 0, 0, 0), ((0.19282602, 0.25572118, -0.11262284), ( 0.75853515, 0.55187444, 0.34649535), 1153689.72115824, 1., 0, 0, 0), ..., ((0.14718062, -0.23794414, -0.17253588), (-0.27354594, 0.15713747, 0.94893648), 350211.6847914 , 1., 0, 0, 0), ((0.14718062, -0.23794414, -0.17253588), ( 0.16444666, -0.98360966, 0.0739549 ), 3259134.69914602, 1., 0, 0, 0), ((0.14718062, -0.23794414, -0.17253588), ( 0.16444666, -0.98360966, 0.0739549 ), 3259134.69914602, 1., 0, 0, 0)], dtype={'names':['r','u','E','wgt','delayed_group','surf_id','particle'], 'formats':[[('x', '<f8'), ('y', '<f8'), ('z', '<f8')],[('x', '<f8'), ('y', '<f8'), ('z', '<f8')],'<f8','<f8','<i4','<i4','<i4'], 'offsets':[0,24,48,56,64,68,72], 'itemsize':96})
If we want, say, only the energies from the source sites, we can simply index the source array with the name of the field:
sp.source['E']
array([2077978.15846043, 206149.19886773, 1153689.72115824, ..., 350211.6847914 , 3259134.69914602, 3259134.69914602])
Now, we can look at things like the energy distribution of source sites. Note that we don't directly use the matplotlib.pyplot.hist
method since our binning is logarithmic.
# Create log-spaced energy bins from 1 keV to 10 MeV
energy_bins = np.logspace(3,7)
# Calculate pdf for source energies
probability, bin_edges = np.histogram(sp.source['E'], energy_bins, density=True)
# Make sure integrating the PDF gives us unity
print(sum(probability*np.diff(energy_bins)))
# Plot source energy PDF
plt.semilogx(energy_bins[:-1], probability*np.diff(energy_bins), drawstyle='steps')
plt.xlabel('Energy (eV)')
plt.ylabel('Probability/eV')
1.0
Text(0, 0.5, 'Probability/eV')
Let's also look at the spatial distribution of the sites. To make the plot a little more interesting, we can also include the direction of the particle emitted from the source and color each source by the logarithm of its energy.
plt.quiver(sp.source['r']['x'], sp.source['r']['y'],
sp.source['u']['x'], sp.source['u']['y'],
np.log(sp.source['E']), cmap='jet', scale=20.0)
plt.colorbar()
plt.xlim((-0.5,0.5))
plt.ylim((-0.5,0.5))
(-0.5, 0.5)
# Close the statepoint file as a matter of best practice
sp.close()