Post Processing¶
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
Generate Input Files¶
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.IndependentSource(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
Tally Data Processing¶
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¶
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()