This IPython Notebook illustrates the use of the openmc.mgxs.Library class. The Library class is designed to automate the calculation of multi-group cross sections for use cases with one or more domains, cross section types, and/or nuclides. In particular, this Notebook illustrates the following features:

  • Calculation of multi-energy-group and multi-delayed-group cross sections for a fuel assembly
  • Automated creation, manipulation and storage of MGXS with openmc.mgxs.Library
  • Steady-state pin-by-pin delayed neutron fractions (beta) for each delayed group.
  • Generation of surface currents on the interfaces and surfaces of a Mesh.

Generate Input Files


In [1]:
%matplotlib inline
import math

import matplotlib.pyplot as plt
import numpy as np

import openmc
import openmc.mgxs

First we need to define materials that will be used in the problem: fuel, water, and cladding.


In [2]:
# 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 object that can be exported to an actual XML file.


In [3]:
# Create a materials collection and export to XML
materials = openmc.Materials((fuel, water, zircaloy))
materials.export_to_xml()

Now let's move on to the geometry. This problem will be a square array of fuel pins and control rod guide tubes for which we can use OpenMC's lattice/universe feature. The basic universe will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces for fuel and clad, as well as the outer bounding surfaces of the problem.


In [4]:
# Create cylinders for the fuel and clad
fuel_outer_radius = openmc.ZCylinder(R=0.39218)
clad_outer_radius = openmc.ZCylinder(R=0.45720)

# Create boundary planes to surround the geometry
min_x = openmc.XPlane(x0=-10.71, boundary_type='reflective')
max_x = openmc.XPlane(x0=+10.71, boundary_type='reflective')
min_y = openmc.YPlane(y0=-10.71, boundary_type='reflective')
max_y = openmc.YPlane(y0=+10.71, boundary_type='reflective')
min_z = openmc.ZPlane(z0=-10., boundary_type='reflective')
max_z = openmc.ZPlane(z0=+10., boundary_type='reflective')

With the surfaces defined, we can now construct a fuel pin cell from cells that are defined by intersections of half-spaces created by the surfaces.


In [5]:
# Create a Universe to encapsulate a fuel pin
fuel_pin_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
fuel_pin_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
fuel_pin_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
fuel_pin_universe.add_cell(moderator_cell)

Likewise, we can construct a control rod guide tube with the same surfaces.


In [6]:
# Create a Universe to encapsulate a control rod guide tube
guide_tube_universe = openmc.Universe(name='Guide Tube')

# Create guide tube Cell
guide_tube_cell = openmc.Cell(name='Guide Tube Water')
guide_tube_cell.fill = water
guide_tube_cell.region = -fuel_outer_radius
guide_tube_universe.add_cell(guide_tube_cell)

# Create a clad Cell
clad_cell = openmc.Cell(name='Guide Clad')
clad_cell.fill = zircaloy
clad_cell.region = +fuel_outer_radius & -clad_outer_radius
guide_tube_universe.add_cell(clad_cell)

# Create a moderator Cell
moderator_cell = openmc.Cell(name='Guide Tube Moderator')
moderator_cell.fill = water
moderator_cell.region = +clad_outer_radius
guide_tube_universe.add_cell(moderator_cell)

Using the pin cell universe, we can construct a 17x17 rectangular lattice with a 1.26 cm pitch.


In [7]:
# Create fuel assembly Lattice
assembly = openmc.RectLattice(name='1.6% Fuel Assembly')
assembly.pitch = (1.26, 1.26)
assembly.lower_left = [-1.26 * 17. / 2.0] * 2

Next, we create a NumPy array of fuel pin and guide tube universes for the lattice.


In [8]:
# Create array indices for guide tube locations in lattice
template_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8,
                       11, 14, 2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
template_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8,
                       8, 11, 11, 11, 11, 11, 13, 13, 14, 14, 14])

# Create universes array with the fuel pin and guide tube universes
universes = np.tile(fuel_pin_universe, (17,17))
universes[template_x, template_y] = guide_tube_universe

# Store the array of universes in the lattice
assembly.universes = universes

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 [9]:
# Create root Cell
root_cell = openmc.Cell(name='root cell', fill=assembly)

# 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 and export it to XML.


In [10]:
# Create Geometry and export to XML
geometry = openmc.Geometry(root_universe)
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 40 active batches each with 2500 particles.


In [11]:
# OpenMC simulation parameters
batches = 50
inactive = 10
particles = 2500

# Instantiate a Settings object
settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.output = {'tallies': False}

# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]
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 to verify that our fuel assembly geometry was created successfully.


In [12]:
# Plot our geometry
plot = openmc.Plot.from_geometry(geometry)
plot.pixels = (250, 250)
plot.color_by = 'material'
openmc.plot_inline(plot)


As we can see from the plot, we have a nice array of fuel and guide tube pin cells with fuel, cladding, and water!

Create an MGXS Library

Now we are ready to generate multi-group cross sections! First, let's define a 20-energy-group and 1-energy-group.


In [13]:
# Instantiate a 20-group EnergyGroups object
energy_groups = openmc.mgxs.EnergyGroups()
energy_groups.group_edges = np.logspace(-3, 7.3, 21)

# Instantiate a 1-group EnergyGroups object
one_group = openmc.mgxs.EnergyGroups()
one_group.group_edges = np.array([energy_groups.group_edges[0], energy_groups.group_edges[-1]])

Next, we will instantiate an openmc.mgxs.Library for the energy and delayed groups with our the fuel assembly geometry.


In [14]:
# Instantiate a tally mesh                                                                                                                                                                  
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [17, 17, 1]
mesh.lower_left = [-10.71, -10.71, -10000.]
mesh.width = [1.26, 1.26, 20000.]

# Initialize an 20-energy-group and 6-delayed-group MGXS Library
mgxs_lib = openmc.mgxs.Library(geometry)
mgxs_lib.energy_groups  = energy_groups
mgxs_lib.num_delayed_groups = 6

# Specify multi-group cross section types to compute
mgxs_lib.mgxs_types = ['total', 'transport', 'nu-scatter matrix', 'kappa-fission', 'inverse-velocity', 'chi-prompt',
                      'prompt-nu-fission', 'chi-delayed', 'delayed-nu-fission', 'beta']

# Specify a "mesh" domain type for the cross section tally filters
mgxs_lib.domain_type = 'mesh'

# Specify the mesh domain over which to compute multi-group cross sections
mgxs_lib.domains = [mesh]

# Construct all tallies needed for the multi-group cross section library
mgxs_lib.build_library()

# Create a "tallies.xml" file for the MGXS Library
tallies_file = openmc.Tallies()
mgxs_lib.add_to_tallies_file(tallies_file, merge=True)

# Instantiate a current tally
mesh_filter = openmc.MeshFilter(mesh)
current_tally = openmc.Tally(name='current tally')
current_tally.scores = ['current']
current_tally.filters = [mesh_filter]

# Add current tally to the tallies file
tallies_file.append(current_tally)

# Export to "tallies.xml"
tallies_file.export_to_xml()

Now, we can run OpenMC to generate the cross sections.


In [15]:
# Run OpenMC
openmc.run()


                               %%%%%%%%%%%%%%%
                          %%%%%%%%%%%%%%%%%%%%%%%%
                       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                   %%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                ###############      %%%%%%%%%%%%%%%%%%%%%%%%
               ##################     %%%%%%%%%%%%%%%%%%%%%%%
               ###################     %%%%%%%%%%%%%%%%%%%%%%%
               ####################     %%%%%%%%%%%%%%%%%%%%%%
               #####################     %%%%%%%%%%%%%%%%%%%%%
               ######################     %%%%%%%%%%%%%%%%%%%%
               #######################     %%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%
                 ####################     %%%%%%%%%%%%%%%%%
                   #################     %%%%%%%%%%%%%%%%%
                    ###############     %%%%%%%%%%%%%%%%
                      ############     %%%%%%%%%%%%%%%
                         ########     %%%%%%%%%%%%%%
                                     %%%%%%%%%%%

                   | The OpenMC Monte Carlo Code
         Copyright | 2011-2017 Massachusetts Institute of Technology
           License | http://openmc.readthedocs.io/en/latest/license.html
           Version | 0.8.0
          Git SHA1 | f7edad68f0654d775ed363bfdcbe4aa5d3cfba23
         Date/Time | 2017-04-03 21:15:56

 Reading settings XML file...
 Reading geometry XML file...
 Reading materials XML file...
 Reading cross sections XML file...
 Reading U235 from /home/romano/openmc/scripts/nndc_hdf5/U235.h5
 Reading U238 from /home/romano/openmc/scripts/nndc_hdf5/U238.h5
 Reading O16 from /home/romano/openmc/scripts/nndc_hdf5/O16.h5
 Reading H1 from /home/romano/openmc/scripts/nndc_hdf5/H1.h5
 Reading B10 from /home/romano/openmc/scripts/nndc_hdf5/B10.h5
 Reading Zr90 from /home/romano/openmc/scripts/nndc_hdf5/Zr90.h5
 Maximum neutron transport energy: 2.00000E+07 eV for U235
 Reading tallies XML file...
 Building neighboring cells lists for each surface...
 Initializing source particles...

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

  Bat./Gen.      k            Average k         
  =========   ========   ====================   
        1/1    1.03852                       
        2/1    0.99743                       
        3/1    1.02987                       
        4/1    1.04397                       
        5/1    1.06262                       
        6/1    1.06657                       
        7/1    0.98574                       
        8/1    1.04364                       
        9/1    1.01253                       
       10/1    1.02094                       
       11/1    0.99586                       
       12/1    1.00508    1.00047 +/- 0.00461
       13/1    1.05292    1.01795 +/- 0.01769
       14/1    1.04732    1.02530 +/- 0.01450
       15/1    1.04886    1.03001 +/- 0.01218
       16/1    1.00948    1.02659 +/- 0.01052
       17/1    1.02684    1.02662 +/- 0.00889
       18/1    0.97234    1.01984 +/- 0.01026
       19/1    0.99754    1.01736 +/- 0.00938
       20/1    0.98964    1.01459 +/- 0.00884
       21/1    1.04140    1.01703 +/- 0.00836
       22/1    1.03854    1.01882 +/- 0.00784
       23/1    1.05917    1.02192 +/- 0.00785
       24/1    1.02413    1.02208 +/- 0.00727
       25/1    1.03113    1.02268 +/- 0.00679
       26/1    1.05113    1.02446 +/- 0.00660
       27/1    1.03252    1.02494 +/- 0.00622
       28/1    1.05196    1.02644 +/- 0.00605
       29/1    0.99663    1.02487 +/- 0.00593
       30/1    1.01820    1.02454 +/- 0.00564
       31/1    1.02753    1.02468 +/- 0.00537
       32/1    1.02162    1.02454 +/- 0.00512
       33/1    1.04083    1.02525 +/- 0.00494
       34/1    1.03335    1.02558 +/- 0.00474
       35/1    1.01304    1.02508 +/- 0.00458
       36/1    0.99299    1.02385 +/- 0.00457
       37/1    1.04936    1.02479 +/- 0.00450
       38/1    1.02856    1.02493 +/- 0.00433
       39/1    1.03706    1.02535 +/- 0.00420
       40/1    1.08118    1.02721 +/- 0.00447
       41/1    1.00149    1.02638 +/- 0.00440
       42/1    1.00233    1.02563 +/- 0.00433
       43/1    1.03023    1.02577 +/- 0.00419
       44/1    1.03230    1.02596 +/- 0.00407
       45/1    0.98123    1.02468 +/- 0.00416
       46/1    1.02126    1.02458 +/- 0.00404
       47/1    0.99772    1.02386 +/- 0.00400
       48/1    1.02773    1.02396 +/- 0.00389
       49/1    1.01690    1.02378 +/- 0.00379
       50/1    1.02890    1.02391 +/- 0.00370
 Creating state point statepoint.50.h5...

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

 Total time for initialization     =  3.7616E-01 seconds
   Reading cross sections          =  3.2363E-01 seconds
 Total time in simulation          =  5.8159E+01 seconds
   Time in transport only          =  5.7959E+01 seconds
   Time in inactive batches        =  3.9349E+00 seconds
   Time in active batches          =  5.4224E+01 seconds
   Time synchronizing fission bank =  3.6903E-03 seconds
     Sampling source sites         =  2.4990E-03 seconds
     SEND/RECV source sites        =  1.1328E-03 seconds
   Time accumulating tallies       =  1.7598E-01 seconds
 Total time for finalization       =  3.5126E-03 seconds
 Total time elapsed                =  5.8551E+01 seconds
 Calculation Rate (inactive)       =  6353.42 neutrons/second
 Calculation Rate (active)         =  1844.20 neutrons/second

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

 k-effective (Collision)     =  1.02621 +/-  0.00393
 k-effective (Track-length)  =  1.02391 +/-  0.00370
 k-effective (Absorption)    =  1.02077 +/-  0.00423
 Combined k-effective        =  1.02331 +/-  0.00353
 Leakage Fraction            =  0.00000 +/-  0.00000

Out[15]:
0

Tally Data Processing

Our simulation ran successfully and created statepoint and summary output files. We begin our analysis by instantiating a StatePoint object.


In [16]:
# Load the last statepoint file
sp = openmc.StatePoint('statepoint.50.h5')

The statepoint is now ready to be analyzed by the Library. We simply have to load the tallies from the statepoint into the Library and our MGXS objects will compute the cross sections for us under-the-hood.


In [17]:
# Initialize MGXS Library with OpenMC statepoint data
mgxs_lib.load_from_statepoint(sp)

# Extrack the current tally separately
current_tally = sp.get_tally(name='current tally')

Using Tally Arithmetic to Compute the Delayed Neutron Precursor Concentrations

Finally, we illustrate how one can leverage OpenMC's tally arithmetic data processing feature with MGXS objects. The openmc.mgxs module uses tally arithmetic to compute multi-group cross sections with automated uncertainty propagation. Each MGXS object includes an xs_tally attribute which is a "derived" Tally based on the tallies needed to compute the cross section type of interest. These derived tallies can be used in subsequent tally arithmetic operations. For example, we can use tally artithmetic to compute the delayed neutron precursor concentrations using the Beta and DelayedNuFissionXS objects. The delayed neutron precursor concentrations are modeled using the following equations:

$$\frac{\partial}{\partial t} C_{k,d} (t) = \int_{0}^{\infty}\mathrm{d}E'\int_{\mathbf{r} \in V_{k}}\mathrm{d}\mathbf{r} \beta_{k,d} (t) \nu_d \sigma_{f,x}(\mathbf{r},E',t)\Phi(\mathbf{r},E',t) - \lambda_{d} C_{k,d} (t) $$$$C_{k,d} (t=0) = \frac{1}{\lambda_{d}} \int_{0}^{\infty}\mathrm{d}E'\int_{\mathbf{r} \in V_{k}}\mathrm{d}\mathbf{r} \beta_{k,d} (t=0) \nu_d \sigma_{f,x}(\mathbf{r},E',t=0)\Phi(\mathbf{r},E',t=0) $$

In [18]:
# Set the time constants for the delayed precursors (in seconds^-1)
precursor_halflife = np.array([55.6, 24.5, 16.3, 2.37, 0.424, 0.195])
precursor_lambda = math.log(2.0) / precursor_halflife

beta = mgxs_lib.get_mgxs(mesh, 'beta')

# Create a tally object with only the delayed group filter for the time constants
beta_filters = [f for f in beta.xs_tally.filters if type(f) is not openmc.DelayedGroupFilter]
lambda_tally = beta.xs_tally.summation(nuclides=beta.xs_tally.nuclides)
for f in beta_filters:
    lambda_tally = lambda_tally.summation(filter_type=type(f), remove_filter=True) * 0. + 1.

# Set the mean of the lambda tally and reshape to account for nuclides and scores
lambda_tally._mean = precursor_lambda
lambda_tally._mean.shape = lambda_tally.std_dev.shape

# Set a total nuclide and lambda score
lambda_tally.nuclides = [openmc.Nuclide(name='total')]
lambda_tally.scores = ['lambda']

delayed_nu_fission = mgxs_lib.get_mgxs(mesh, 'delayed-nu-fission')

# Use tally arithmetic to compute the precursor concentrations
precursor_conc = beta.xs_tally.summation(filter_type=openmc.EnergyFilter, remove_filter=True) * \
    delayed_nu_fission.xs_tally.summation(filter_type=openmc.EnergyFilter, remove_filter=True) / lambda_tally
    
# The difference is a derived tally which can generate Pandas DataFrames for inspection
precursor_conc.get_pandas_dataframe().head(10)


/home/romano/openmc/openmc/tallies.py:1875: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
/home/romano/openmc/openmc/tallies.py:1876: RuntimeWarning: invalid value encountered in true_divide
  other_rel_err = data['other']['std. dev.'] / data['other']['mean']
/home/romano/openmc/openmc/tallies.py:1877: RuntimeWarning: invalid value encountered in true_divide
  new_tally._mean = data['self']['mean'] / data['other']['mean']
/home/romano/openmc/openmc/tallies.py:1869: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
/home/romano/openmc/openmc/tallies.py:1870: RuntimeWarning: invalid value encountered in true_divide
  other_rel_err = data['other']['std. dev.'] / data['other']['mean']
Out[18]:
mesh 1 delayedgroup nuclide score mean std. dev.
x y z
0 1 1 1 1 total (((delayed-nu-fission / nu-fission) * (delayed... 0.011301 0.003184
1 1 1 1 2 total (((delayed-nu-fission / nu-fission) * (delayed... 0.000854 0.000078
2 1 1 1 3 total (((delayed-nu-fission / nu-fission) * (delayed... 0.007267 0.000713
3 1 1 1 4 total (((delayed-nu-fission / nu-fission) * (delayed... 0.088903 0.005501
4 1 1 1 5 total (((delayed-nu-fission / nu-fission) * (delayed... 0.044257 0.002890
5 1 1 1 6 total (((delayed-nu-fission / nu-fission) * (delayed... 0.001897 0.000282
6 1 2 1 1 total (((delayed-nu-fission / nu-fission) * (delayed... 0.007362 0.001968
7 1 2 1 2 total (((delayed-nu-fission / nu-fission) * (delayed... 0.000872 0.000064
8 1 2 1 3 total (((delayed-nu-fission / nu-fission) * (delayed... 0.007015 0.000702
9 1 2 1 4 total (((delayed-nu-fission / nu-fission) * (delayed... 0.101115 0.007008

Another useful feature of the Python API is the ability to extract the surface currents for the interfaces and surfaces of a mesh. We can inspect the currents for the mesh by getting the pandas dataframe.


In [19]:
current_tally.get_pandas_dataframe().head(10)


Out[19]:
mesh 1 surface nuclide score mean std. dev.
x y z
0 1 1 1 x-min out total current 0.00000 0.000000
1 1 1 1 x-max out total current 0.00000 0.000000
2 1 1 1 y-min out total current 0.03153 0.000566
3 1 1 1 y-max out total current 0.03153 0.000593
4 1 1 1 z-min out total current 0.00000 0.000000
5 1 1 1 z-max out total current 0.00000 0.000000
6 1 1 1 x-min in total current 0.03265 0.000753
7 1 1 1 x-max in total current 0.03216 0.000698
8 1 1 1 y-min in total current 0.00000 0.000000
9 1 1 1 y-max in total current 0.00000 0.000000

Cross Section Visualizations

In addition to inspecting the data in the tallies by getting the pandas dataframe, we can also plot the tally data on the domain mesh. Below is the delayed neutron fraction tallied in each mesh cell for each delayed group.


In [20]:
# Extract the energy-condensed delayed neutron fraction tally
beta_by_group = beta.get_condensed_xs(one_group).xs_tally.summation(filter_type='energy', remove_filter=True)
beta_by_group.mean.shape = (17, 17, 6)
beta_by_group.mean[beta_by_group.mean == 0] = np.nan

# Plot the betas
plt.figure(figsize=(18,9))
fig = plt.subplot(231)
plt.imshow(beta_by_group.mean[:,:,0], interpolation='none', cmap='jet')
plt.colorbar()
plt.title('Beta - delayed group 1')

fig = plt.subplot(232)
plt.imshow(beta_by_group.mean[:,:,1], interpolation='none', cmap='jet')
plt.colorbar()
plt.title('Beta - delayed group 2')

fig = plt.subplot(233)
plt.imshow(beta_by_group.mean[:,:,2], interpolation='none', cmap='jet')
plt.colorbar()
plt.title('Beta - delayed group 3')

fig = plt.subplot(234)
plt.imshow(beta_by_group.mean[:,:,3], interpolation='none', cmap='jet')
plt.colorbar()
plt.title('Beta - delayed group 4')

fig = plt.subplot(235)
plt.imshow(beta_by_group.mean[:,:,4], interpolation='none', cmap='jet')
plt.colorbar()
plt.title('Beta - delayed group 5')

fig = plt.subplot(236)
plt.imshow(beta_by_group.mean[:,:,5], interpolation='none', cmap='jet')
plt.colorbar()
plt.title('Beta - delayed group 6')


/home/romano/openmc/openmc/tallies.py:1875: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
/home/romano/openmc/openmc/tallies.py:1876: RuntimeWarning: invalid value encountered in true_divide
  other_rel_err = data['other']['std. dev.'] / data['other']['mean']
/home/romano/openmc/openmc/tallies.py:1877: RuntimeWarning: invalid value encountered in true_divide
  new_tally._mean = data['self']['mean'] / data['other']['mean']
Out[20]:
<matplotlib.text.Text at 0x7fecd32c1da0>