This notebook shows the how tallies can be combined (added, subtracted, multiplied, etc.) using the Python API in order to create derived tallies. Since no covariance information is obtained, it is assumed that tallies are completely independent of one another when propagating uncertainties. The target problem is a simple pin cell.

Note: that this Notebook was created using the latest Pandas v0.16.1. Everything in the Notebook will wun with older versions of Pandas, but the multi-indexing option in >v0.15.0 makes the tables look prettier.


In [1]:
%load_ext autoreload
%autoreload 2

In [2]:
import glob
from IPython.display import Image
import numpy as np

import openmc
from openmc.statepoint import StatePoint
from openmc.summary import Summary
from openmc.source import Source
from openmc.stats import Box

%matplotlib inline

Generate Input Files

First we need to define materials that will be used in the problem. Before defining a material, we must create nuclides that are used in the material.


In [3]:
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
b10 = openmc.Nuclide('B-10')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
u238 = openmc.Nuclide('U-238')
zr90 = openmc.Nuclide('Zr-90')

With the nuclides we defined, we will now create three materials for the fuel, water, and cladding of the fuel pin.


In [4]:
# 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.


In [5]:
# Instantiate a MaterialsFile, add Materials
materials_file = openmc.MaterialsFile()
materials_file.add_material(fuel)
materials_file.add_material(water)
materials_file.add_material(zircaloy)
materials_file.default_xs = '71c'

# Export to "materials.xml"
materials_file.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.


In [6]:
# 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
# Use both reflective and vacuum boundaries to make life interesting
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.


In [7]:
# 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.


In [8]:
# 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.


In [9]:
# Create Geometry and set root Universe
geometry = openmc.Geometry()
geometry.root_universe = root_universe

In [10]:
# Instantiate a GeometryFile
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry

# Export to "geometry.xml"
geometry_file.export_to_xml()

With the geometry and materials finished, we now just need to define simulation parameters. In this case, we will use 5 inactive batches and 15 active batches each with 2500 particles.


In [11]:
# OpenMC simulation parameters
batches = 20
inactive = 5
particles = 2500

# Instantiate a SettingsFile
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.output = {'tallies': True, 'summary': True}
source_bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]
settings_file.source = Source(space=Box(
        source_bounds[:3], source_bounds[3:]))

# Export to "settings.xml"
settings_file.export_to_xml()

Let us also create a plot file that we can use to verify that our pin cell geometry was created successfully.


In [12]:
# Instantiate a Plot
plot = openmc.Plot(plot_id=1)
plot.filename = 'materials-xy'
plot.origin = [0, 0, 0]
plot.width = [1.26, 1.26]
plot.pixels = [250, 250]
plot.color = 'mat'

# Instantiate a PlotsFile, add Plot, and export to "plots.xml"
plot_file = openmc.PlotsFile()
plot_file.add_plot(plot)
plot_file.export_to_xml()

With the plots.xml file, we can now generate and view the plot. OpenMC outputs plots in .ppm format, which can be converted into a compressed format like .png with the convert utility.


In [13]:
# Run openmc in plotting mode
executor = openmc.Executor()
executor.plot_geometry(output=False)


Out[13]:
0

In [14]:
# Convert OpenMC's funky ppm to png
!convert materials-xy.ppm materials-xy.png

# Display the materials plot inline
Image(filename='materials-xy.png')


Out[14]:

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 variety of tallies.


In [15]:
# Instantiate an empty TalliesFile
tallies_file = openmc.TalliesFile()

In [16]:
# Create Tallies to compute microscopic multi-group cross-sections

# Instantiate energy filter for multi-group cross-section Tallies
energy_filter = openmc.Filter(type='energy', bins=[0., 0.625e-6, 20.])

# Instantiate flux Tally in moderator and fuel
tally = openmc.Tally(name='flux')
tally.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id, moderator_cell.id]))
tally.add_filter(energy_filter)
tally.add_score('flux')
tallies_file.add_tally(tally)

# Instantiate reaction rate Tally in fuel
tally = openmc.Tally(name='fuel rxn rates')
tally.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id]))
tally.add_filter(energy_filter)
tally.add_score('nu-fission')
tally.add_score('scatter')
tally.add_nuclide(u238)
tally.add_nuclide(u235)
tallies_file.add_tally(tally)

# Instantiate reaction rate Tally in moderator
tally = openmc.Tally(name='moderator rxn rates')
tally.add_filter(openmc.Filter(type='cell', bins=[moderator_cell.id]))
tally.add_filter(energy_filter)
tally.add_score('absorption')
tally.add_score('total')
tally.add_nuclide(o16)
tally.add_nuclide(h1)
tallies_file.add_tally(tally)

In [17]:
# K-Eigenvalue (infinity) tallies
fiss_rate = openmc.Tally(name='fiss. rate')
abs_rate = openmc.Tally(name='abs. rate')
fiss_rate.add_score('nu-fission')
abs_rate.add_score('absorption')
tallies_file.add_tally(fiss_rate)
tallies_file.add_tally(abs_rate)

In [18]:
# Resonance Escape Probability tallies
therm_abs_rate = openmc.Tally(name='therm. abs. rate')
therm_abs_rate.add_score('absorption')
therm_abs_rate.add_filter(openmc.Filter(type='energy', bins=[0., 0.625]))
tallies_file.add_tally(therm_abs_rate)

In [19]:
# Thermal Flux Utilization tallies
fuel_therm_abs_rate = openmc.Tally(name='fuel therm. abs. rate')
fuel_therm_abs_rate.add_score('absorption')
fuel_therm_abs_rate.add_filter(openmc.Filter(type='energy', bins=[0., 0.625]))
fuel_therm_abs_rate.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id]))
tallies_file.add_tally(fuel_therm_abs_rate)

In [20]:
# Fast Fission Factor tallies
therm_fiss_rate = openmc.Tally(name='therm. fiss. rate')
therm_fiss_rate.add_score('nu-fission')
therm_fiss_rate.add_filter(openmc.Filter(type='energy', bins=[0., 0.625]))
tallies_file.add_tally(therm_fiss_rate)

In [21]:
# Instantiate energy filter to illustrate Tally slicing
energy_filter = openmc.Filter(type='energy', bins=np.logspace(np.log10(1e-8), np.log10(20), 10))

# Instantiate flux Tally in moderator and fuel
tally = openmc.Tally(name='need-to-slice')
tally.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id, moderator_cell.id]))
tally.add_filter(energy_filter)
tally.add_score('nu-fission')
tally.add_score('scatter')
tally.add_nuclide(h1)
tally.add_nuclide(u238)
tallies_file.add_tally(tally)

In [22]:
# Export to "tallies.xml"
tallies_file.export_to_xml()

Now we a have a complete set of inputs, so we can go ahead and run our simulation.


In [23]:
# Remove old HDF5 (summary, statepoint) files
!rm statepoint.*

# Run OpenMC with MPI!
executor.run_simulation()


       .d88888b.                             888b     d888  .d8888b.
      d88P" "Y88b                            8888b   d8888 d88P  Y88b
      888     888                            88888b.d88888 888    888
      888     888 88888b.   .d88b.  88888b.  888Y88888P888 888       
      888     888 888 "88b d8P  Y8b 888 "88b 888 Y888P 888 888       
      888     888 888  888 88888888 888  888 888  Y8P  888 888    888
      Y88b. .d88P 888 d88P Y8b.     888  888 888   "   888 Y88b  d88P
       "Y88888P"  88888P"   "Y8888  888  888 888       888  "Y8888P"
__________________888______________________________________________________
                  888
                  888

      Copyright:      2011-2015 Massachusetts Institute of Technology
      License:        http://mit-crpg.github.io/openmc/license.html
      Version:        0.7.1
      Git SHA1:       ea9fb637f63f9374c7436456141afa850b84acf9
      Date/Time:      2016-01-14 07:00:14

 ===========================================================================
 ========================>     INITIALIZATION     <=========================
 ===========================================================================

 Reading settings XML file...
 Reading cross sections XML file...
 Reading geometry XML file...
 Reading materials XML file...
 Reading tallies XML file...
 Building neighboring cells lists for each surface...
 Loading ACE cross section table: 92235.71c
 Loading ACE cross section table: 92238.71c
 Loading ACE cross section table: 8016.71c
 Loading ACE cross section table: 1001.71c
 Loading ACE cross section table: 5010.71c
 Loading ACE cross section table: 40090.71c
 Maximum neutron transport energy: 20.0000 MeV for 92235.71c
 Initializing source particles...

 ===========================================================================
 ====================>     K EIGENVALUE SIMULATION     <====================
 ===========================================================================

  Bat./Gen.      k            Average k         
  =========   ========   ====================   
        1/1    1.05992                       
        2/1    1.05251                       
        3/1    1.05204                       
        4/1    1.02100                       
        5/1    1.07784                       
        6/1    1.04814                       
        7/1    1.02335    1.03574 +/- 0.01239
        8/1    1.02415    1.03188 +/- 0.00813
        9/1    1.10331    1.04974 +/- 0.01876
       10/1    1.05452    1.05069 +/- 0.01456
       11/1    1.07867    1.05536 +/- 0.01277
       12/1    1.04203    1.05345 +/- 0.01096
       13/1    1.04482    1.05237 +/- 0.00955
       14/1    1.04117    1.05113 +/- 0.00852
       15/1    1.07581    1.05360 +/- 0.00801
       16/1    1.04235    1.05257 +/- 0.00731
       17/1    1.02710    1.05045 +/- 0.00701
       18/1    1.01970    1.04809 +/- 0.00687
       19/1    1.01022    1.04538 +/- 0.00691
       20/1    1.01449    1.04332 +/- 0.00675
 Creating state point statepoint.20.h5...

 ===========================================================================
 ======================>     SIMULATION FINISHED     <======================
 ===========================================================================


 =======================>     TIMING STATISTICS     <=======================

 Total time for initialization     =  1.2510E+00 seconds
   Reading cross sections          =  9.7600E-01 seconds
 Total time in simulation          =  1.5844E+01 seconds
   Time in transport only          =  1.5834E+01 seconds
   Time in inactive batches        =  2.2840E+00 seconds
   Time in active batches          =  1.3560E+01 seconds
   Time synchronizing fission bank =  3.0000E-03 seconds
     Sampling source sites         =  2.0000E-03 seconds
     SEND/RECV source sites        =  1.0000E-03 seconds
   Time accumulating tallies       =  0.0000E+00 seconds
 Total time for finalization       =  1.0000E-03 seconds
 Total time elapsed                =  1.7110E+01 seconds
 Calculation Rate (inactive)       =  5472.85 neutrons/second
 Calculation Rate (active)         =  2765.49 neutrons/second

 ============================>     RESULTS     <============================

 k-effective (Collision)     =  1.03935 +/-  0.00682
 k-effective (Track-length)  =  1.04332 +/-  0.00675
 k-effective (Absorption)    =  1.03845 +/-  0.00598
 Combined k-effective        =  1.04024 +/-  0.00523
 Leakage Fraction            =  0.00000 +/-  0.00000

Out[23]:
0

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, the tally results are not read into memory because they might be large, even large enough to exceed the available memory on a computer.


In [24]:
# Load the statepoint file
sp = StatePoint('statepoint.20.h5')

You may have also noticed we instructed OpenMC to create a summary file with lots of geometry information in it. This can help to produce more sensible output from the Python API, so we will use the summary file to link against.


In [25]:
# Load the summary file and link with statepoint
su = Summary('summary.h5')
sp.link_with_summary(su)

We have a tally of the total fission rate and the total absorption rate, so we can calculate k-infinity as: $$k_\infty = \frac{\langle \nu \Sigma_f \phi \rangle}{\langle \Sigma_a \phi \rangle}$$ In this notation, $\langle \cdot \rangle^a_b$ represents an OpenMC that is integrated over region $a$ and energy range $b$. If $a$ or $b$ is not reported, it means the value represents an integral over all space or all energy, respectively.


In [26]:
# Compute k-infinity using tally arithmetic
fiss_rate = sp.get_tally(name='fiss. rate')
abs_rate = sp.get_tally(name='abs. rate')
keff = fiss_rate / abs_rate
keff.get_pandas_dataframe()


Out[26]:
nuclide score mean std. dev.
0 total (nu-fission / absorption) 1.040166 0.009069

Notice that even though the neutron production rate and absorption rate are separate tallies, we still get a first-order estimate of the uncertainty on the quotient of them automatically!

Often in textbooks you'll see k-infinity represented using the four-factor formula $$k_\infty = p \epsilon f \eta.$$ Let's analyze each of these factors, starting with the resonance escape probability which is defined as $$p=\frac{\langle\Sigma_a\phi\rangle_T}{\langle\Sigma_a\phi\rangle}$$ where the subscript $T$ means thermal energies.


In [27]:
# Compute resonance escape probability using tally arithmetic
therm_abs_rate = sp.get_tally(name='therm. abs. rate')
res_esc = therm_abs_rate / abs_rate
res_esc.get_pandas_dataframe()


Out[27]:
energy [MeV] nuclide score mean std. dev.
0 (0.0e+00 - 6.2e-01) total absorption 0.95938 0.008187

The fast fission factor can be calculated as $$\epsilon=\frac{\langle\nu\Sigma_f\phi\rangle}{\langle\nu\Sigma_f\phi\rangle_T}$$


In [28]:
# Compute fast fission factor factor using tally arithmetic
therm_fiss_rate = sp.get_tally(name='therm. fiss. rate')
fast_fiss = fiss_rate / therm_fiss_rate
fast_fiss.get_pandas_dataframe()


Out[28]:
energy [MeV] nuclide score mean std. dev.
0 (0.0e+00 - 6.2e-01) total nu-fission 1.090899 0.010602

The thermal flux utilization is calculated as $$f=\frac{\langle\Sigma_a\phi\rangle^F_T}{\langle\Sigma_a\phi\rangle_T}$$ where the superscript $F$ denotes fuel.


In [29]:
# Compute thermal flux utilization factor using tally arithmetic
fuel_therm_abs_rate = sp.get_tally(name='fuel therm. abs. rate')
therm_util = fuel_therm_abs_rate / therm_abs_rate
therm_util.get_pandas_dataframe()


Out[29]:
energy [MeV] cell nuclide score mean std. dev.
0 (0.0e+00 - 6.2e-01) 10000 total absorption 0.803413 0.007031

The final factor is the number of fission neutrons produced per absorption in fuel, calculated as $$\eta = \frac{\langle \nu\Sigma_f\phi \rangle_T}{\langle \Sigma_a \phi \rangle^F_T}$$


In [30]:
# Compute neutrons produced per absorption (eta) using tally arithmetic
eta = therm_fiss_rate / fuel_therm_abs_rate
eta.get_pandas_dataframe()


Out[30]:
energy [MeV] cell nuclide score mean std. dev.
0 (0.0e+00 - 6.2e-01) 10000 total (nu-fission / absorption) 1.237053 0.011765

Now we can calculate $k_\infty$ using the product of the factors form the four-factor formula.


In [31]:
keff = res_esc * fast_fiss * therm_util * eta
keff.get_pandas_dataframe()


Out[31]:
energy [MeV] cell nuclide score mean std. dev.
0 (0.0e+00 - 6.2e-01) 10000 total (((absorption * nu-fission) * absorption) * (n... 1.040166 0.019018

We see that the value we've obtained here has exactly the same mean as before. However, because of the way it was calculated, the standard deviation appears to be larger.

Let's move on to a more complicated example now. Before we set up tallies to get reaction rates in the fuel and moderator in two energy groups for two different nuclides. We can use tally arithmetic to divide each of these reaction rates by the flux to get microscopic multi-group cross sections.


In [32]:
# Compute microscopic multi-group cross-sections
flux = sp.get_tally(name='flux')
flux = flux.get_slice(filters=['cell'], filter_bins=[(fuel_cell.id,)])
fuel_rxn_rates = sp.get_tally(name='fuel rxn rates')
mod_rxn_rates = sp.get_tally(name='moderator rxn rates')

In [33]:
fuel_xs = fuel_rxn_rates / flux
fuel_xs.get_pandas_dataframe()


Out[33]:
cell energy [MeV] nuclide score mean std. dev.
0 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (nu-fission / flux) 0.000001 7.377419e-09
1 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (scatter / flux) 0.209989 2.303838e-03
2 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (nu-fission / flux) 0.356420 3.951669e-03
3 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (scatter / flux) 0.005555 6.101004e-05
4 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (nu-fission / flux) 0.007155 8.053460e-05
5 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (scatter / flux) 0.227770 1.079289e-03
6 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (nu-fission / flux) 0.008067 5.254797e-05
7 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (scatter / flux) 0.003367 1.647058e-05

We see that when the two tallies with multiple bins were divided, the derived tally contains the outer product of the combinations. If the filters/scores are the same, no outer product is needed. The get_values(...) method allows us to obtain a subset of tally scores. In the following example, we obtain just the neutron production microscopic cross sections.


In [34]:
# Show how to use Tally.get_values(...) with a CrossScore
nu_fiss_xs = fuel_xs.get_values(scores=['(nu-fission / flux)'])
print(nu_fiss_xs)


[[[  6.65702880e-07]
  [  3.56420449e-01]]

 [[  7.15488656e-03]
  [  8.06673774e-03]]]

The same idea can be used not only for scores but also for filters and nuclides.


In [35]:
# Show how to use Tally.get_values(...) with a CrossScore and CrossNuclide
u235_scatter_xs = fuel_xs.get_values(nuclides=['(U-235 / total)'], 
                                scores=['(scatter / flux)'])
print(u235_scatter_xs)


[[[ 0.00555533]]

 [[ 0.0033668 ]]]

In [36]:
# Show how to use Tally.get_values(...) with a CrossFilter and CrossScore
fast_scatter_xs = fuel_xs.get_values(filters=['energy'], 
                                     filter_bins=[((0.625e-6, 20.),)], 
                                     scores=['(scatter / flux)'])
print(fast_scatter_xs)


[[[ 0.22777006]
  [ 0.0033668 ]]]

A more advanced method is to use get_slice(...) to create a new derived tally that is a subset of an existing tally. This has the benefit that we can use get_pandas_dataframe() to see the tallies in a more human-readable format.


In [37]:
# "Slice" the nu-fission data into a new derived Tally
nu_fission_rates = fuel_rxn_rates.get_slice(scores=['nu-fission'])
nu_fission_rates.get_pandas_dataframe()


Out[37]:
cell energy [MeV] nuclide score mean std. dev.
0 10000 (0.0e+00 - 6.3e-07) U-238 nu-fission 0.000002 1.283958e-08
1 10000 (0.0e+00 - 6.3e-07) U-235 nu-fission 0.868553 6.880390e-03
2 10000 (6.3e-07 - 2.0e+01) U-238 nu-fission 0.082149 8.837250e-04
3 10000 (6.3e-07 - 2.0e+01) U-235 nu-fission 0.092618 5.195308e-04

In [38]:
# "Slice" the H-1 scatter data in the moderator Cell into a new derived Tally
need_to_slice = sp.get_tally(name='need-to-slice')
slice_test = need_to_slice.get_slice(scores=['scatter'], nuclides=['H-1'],
                                 filters=['cell'], filter_bins=[(moderator_cell.id,)])
slice_test.get_pandas_dataframe()


Out[38]:
cell energy [MeV] nuclide score mean std. dev.
0 10002 (1.0e-08 - 1.1e-07) H-1 scatter 4.619398 0.040124
1 10002 (1.1e-07 - 1.2e-06) H-1 scatter 2.030757 0.011239
2 10002 (1.2e-06 - 1.3e-05) H-1 scatter 1.658488 0.009777
3 10002 (1.3e-05 - 1.4e-04) H-1 scatter 1.853002 0.007378
4 10002 (1.4e-04 - 1.5e-03) H-1 scatter 2.050773 0.012484
5 10002 (1.5e-03 - 1.6e-02) H-1 scatter 2.131759 0.007821
6 10002 (1.6e-02 - 1.7e-01) H-1 scatter 2.213710 0.015159
7 10002 (1.7e-01 - 1.9e+00) H-1 scatter 2.011925 0.009406
8 10002 (1.9e+00 - 2.0e+01) H-1 scatter 0.371280 0.003949