MGXS Part III: LibrariesΒΆ

mgxs-part-iii

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-group cross sections for a fuel assembly
  • Automated creation, manipulation and storage of MGXS with openmc.mgxs.Library
  • Validation of multi-group cross sections with OpenMOC
  • Steady-state pin-by-pin fission rates comparison between OpenMC and OpenMOC

Note: This Notebook was created using OpenMOC to verify the multi-group cross-sections generated by OpenMC. You must install OpenMOC on your system to run this Notebook in its entirety. In addition, this Notebook illustrates the use of Pandas DataFrames to containerize multi-group cross section data.

Generate Input Files

In [1]:
import math
import pickle

from IPython.display import Image
import matplotlib.pyplot as plt
import numpy as np

import openmc
import openmc.mgxs
from openmc.openmoc_compatible import get_openmoc_geometry
import openmoc
import openmoc.process
from openmoc.materialize import load_openmc_mgxs_lib

%matplotlib inline
/home/wbinventor/miniconda3/lib/python3.5/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)

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 [2]:
# Instantiate some Nuclides
h1 = openmc.Nuclide('H1')
b10 = openmc.Nuclide('B10')
o16 = openmc.Nuclide('O16')
u235 = openmc.Nuclide('U235')
u238 = openmc.Nuclide('U238')
zr90 = openmc.Nuclide('Zr90')

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

In [3]:
# 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 [4]:
# Instantiate a Materials object
materials_file = openmc.Materials((fuel, water, zircaloy))

# Export to "materials.xml"
materials_file.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 [5]:
# 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=-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 [6]:
# 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 [7]:
# 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 [8]:
# 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 [9]:
# 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])

# Initialize an empty 17x17 array of the lattice universes
universes = np.empty((17, 17), dtype=openmc.Universe)

# Fill the array with the fuel pin and guide tube universes
universes[:,:] = fuel_pin_universe
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 [10]:
# Create root Cell
root_cell = openmc.Cell(name='root cell')
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 [11]:
# Create Geometry and set root Universe
geometry = openmc.Geometry()
geometry.root_universe = root_universe
In [12]:
# 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 40 active batches each with 2500 particles.

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

# Instantiate a Settings object
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.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_file.source = openmc.source.Source(space=uniform_dist)

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

Let us also create a Plots file that we can use to verify that our fuel assembly geometry was created successfully.

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

# Instantiate a Plots object, add Plot, and export to "plots.xml"
plot_file = openmc.Plots([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 [15]:
# Run openmc in plotting mode
openmc.plot_geometry(output=False)
Out[15]:
0
In [16]:
# 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[16]:

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 2-group structure using the built-in EnergyGroups class.

In [17]:
# Instantiate a 2-group EnergyGroups object
groups = openmc.mgxs.EnergyGroups()
groups.group_edges = np.array([0., 0.625, 20.0e6])

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

In [18]:
# Initialize an 2-group MGXS Library for OpenMOC
mgxs_lib = openmc.mgxs.Library(geometry)
mgxs_lib.energy_groups = groups

Now, we must specify to the Library which types of cross sections to compute. In particular, the following are the multi-group cross section MGXS subclasses that are mapped to string codes accepted by the Library class:

  • TotalXS ("total")
  • TransportXS ("transport" or "nu-transport with nu set to True)
  • AbsorptionXS ("absorption")
  • CaptureXS ("capture")
  • FissionXS ("fission" or "nu-fission" with nu set to True)
  • KappaFissionXS ("kappa-fission")
  • ScatterXS ("scatter" or "nu-scatter" with nu set to True)
  • ScatterMatrixXS ("scatter matrix" or "nu-scatter matrix" with nu set to True)
  • Chi ("chi")
  • ChiPrompt ("chi prompt")
  • InverseVelocity ("inverse-velocity")
  • PromptNuFissionXS ("prompt-nu-fission")
  • DelayedNuFissionXS ("delayed-nu-fission")
  • ChiDelayed ("chi-delayed")
  • Beta ("beta")

In this case, let's create the multi-group cross sections needed to run an OpenMOC simulation to verify the accuracy of our cross sections. In particular, we will define "transport", "nu-fission", '"fission", "nu-scatter matrix" and "chi" cross sections for our Library.

Note: A variety of different approximate transport-corrected total multi-group cross sections (and corresponding scattering matrices) can be found in the literature. At the present time, the openmc.mgxs module only supports the "P0" transport correction. This correction can be turned on and off through the boolean Library.correction property which may take values of "P0" (default) or None.

In [19]:
# Specify multi-group cross section types to compute
mgxs_lib.mgxs_types = ['transport', 'nu-fission', 'fission', 'nu-scatter matrix', 'chi']

Now we must specify the type of domain over which we would like the Library to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the Library supports "material", "cell", "universe", and "mesh" domain types. We will use a "cell" domain type here to compute cross sections in each of the cells in the fuel assembly geometry.

Note: By default, the Library class will instantiate MGXS objects for each and every domain (material, cell or universe) in the geometry of interest. However, one may specify a subset of these domains to the Library.domains property. In our case, we wish to compute multi-group cross sections in each and every cell since they will be needed in our downstream OpenMOC calculation on the identical combinatorial geometry mesh.

In [20]:
# Specify a "cell" domain type for the cross section tally filters
mgxs_lib.domain_type = 'cell'

# Specify the cell domains over which to compute multi-group cross sections
mgxs_lib.domains = geometry.get_all_material_cells().values()

We can easily instruct the Library to compute multi-group cross sections on a nuclide-by-nuclide basis with the boolean Library.by_nuclide property. By default, by_nuclide is set to False, but we will set it to True here.

In [21]:
# Compute cross sections on a nuclide-by-nuclide basis
mgxs_lib.by_nuclide = True

Lastly, we use the Library to construct the tallies needed to compute all of the requested multi-group cross sections in each domain and nuclide.

In [22]:
# Construct all tallies needed for the multi-group cross section library
mgxs_lib.build_library()

The tallies can now be export to a "tallies.xml" input file for OpenMC.

NOTE: At this point the Library has constructed nearly 100 distinct Tally objects. The overhead to tally in OpenMC scales as $O(N)$ for $N$ tallies, which can become a bottleneck for large tally datasets. To compensate for this, the Python API's Tally, Filter and Tallies classes allow for the smart merging of tallies when possible. The Library class supports this runtime optimization with the use of the optional merge paramter (False by default) for the Library.add_to_tallies_file(...) method, as shown below.

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

In addition, we instantiate a fission rate mesh tally to compare with OpenMOC.

In [24]:
# Instantiate a tally Mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [17, 17]
mesh.lower_left = [-10.71, -10.71]
mesh.upper_right = [+10.71, +10.71]

# Instantiate tally Filter
mesh_filter = openmc.MeshFilter(mesh)

# Instantiate the Tally
tally = openmc.Tally(name='mesh tally')
tally.filters = [mesh_filter]
tally.scores = ['fission', 'nu-fission']

# Add tally to collection
tallies_file.append(tally)
In [25]:
# Export all tallies to a "tallies.xml" file
tallies_file.export_to_xml()
In [26]:
# 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 | 647bf77a57a3cc5cce24b39cb192e1b99f52e499
         Date/Time | 2017-02-27 14:21:38
    OpenMP Threads | 4

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

 Reading settings XML file...
 Reading geometry XML file...
 Reading materials XML file...
 Reading cross sections XML file...
 Reading U235 from
 /home/wbinventor/Documents/NSE-CRPG-Codes/openmc/data/nndc_hdf5/U235.h5
 Reading U238 from
 /home/wbinventor/Documents/NSE-CRPG-Codes/openmc/data/nndc_hdf5/U238.h5
 Reading O16 from
 /home/wbinventor/Documents/NSE-CRPG-Codes/openmc/data/nndc_hdf5/O16.h5
 Reading H1 from
 /home/wbinventor/Documents/NSE-CRPG-Codes/openmc/data/nndc_hdf5/H1.h5
 Reading B10 from
 /home/wbinventor/Documents/NSE-CRPG-Codes/openmc/data/nndc_hdf5/B10.h5
 Reading Zr90 from
 /home/wbinventor/Documents/NSE-CRPG-Codes/openmc/data/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...

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


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

 Total time for initialization     =  3.4887E-01 seconds
   Reading cross sections          =  2.1990E-01 seconds
 Total time in simulation          =  3.2195E+01 seconds
   Time in transport only          =  3.1778E+01 seconds
   Time in inactive batches        =  1.9903E+00 seconds
   Time in active batches          =  3.0205E+01 seconds
   Time synchronizing fission bank =  5.9614E-03 seconds
     Sampling source sites         =  4.8344E-03 seconds
     SEND/RECV source sites        =  1.0392E-03 seconds
   Time accumulating tallies       =  1.5849E-03 seconds
 Total time for finalization       =  3.9664E-05 seconds
 Total time elapsed                =  3.2560E+01 seconds
 Calculation Rate (inactive)       =  12561.1 neutrons/second
 Calculation Rate (active)         =  3310.69 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[26]:
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 [27]:
# 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 [28]:
# Initialize MGXS Library with OpenMC statepoint data
mgxs_lib.load_from_statepoint(sp)

Voila! Our multi-group cross sections are now ready to rock 'n roll!

Extracting and Storing MGXS Data

The Library supports a rich API to automate a variety of tasks, including multi-group cross section data retrieval and storage. We will highlight a few of these features here. First, the Library.get_mgxs(...) method allows one to extract an MGXS object from the Library for a particular domain and cross section type. The following cell illustrates how one may extract the NuFissionXS object for the fuel cell.

Note: The MGXS.get_mgxs(...) method will accept either the domain or the integer domain ID of interest.

In [29]:
# Retrieve the NuFissionXS object for the fuel cell from the library
fuel_mgxs = mgxs_lib.get_mgxs(fuel_cell, 'nu-fission')

The NuFissionXS object supports all of the methods described previously in the openmc.mgxs tutorials, such as Pandas DataFrames: Note that since so few histories were simulated, we should expect a few division-by-error errors as some tallies have not yet scored any results.

In [30]:
df = fuel_mgxs.get_pandas_dataframe()
df
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1835: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
Out[30]:
cell group in nuclide mean std. dev.
3 10000 1 U235 8.096764e-03 3.130177e-05
4 10000 1 U238 7.364515e-03 4.510564e-05
5 10000 1 O16 0.000000e+00 0.000000e+00
0 10000 2 U235 3.611153e-01 2.048312e-03
1 10000 2 U238 6.735070e-07 3.780177e-09
2 10000 2 O16 0.000000e+00 0.000000e+00

Similarly, we can use the MGXS.print_xs(...) method to view a string representation of the multi-group cross section data.

In [31]:
fuel_mgxs.print_xs()
Multi-Group XS
	Reaction Type  =	nu-fission
	Domain Type    =	cell
	Domain ID      =	10000
	Nuclide        =	U235
	Cross Sections [cm^-1]:
            Group 1 [0.625      - 20000000.0eV]:	8.10e-03 +/- 3.87e-01%
            Group 2 [0.0        - 0.625     eV]:	3.61e-01 +/- 5.67e-01%

	Nuclide        =	U238
	Cross Sections [cm^-1]:
            Group 1 [0.625      - 20000000.0eV]:	7.36e-03 +/- 6.12e-01%
            Group 2 [0.0        - 0.625     eV]:	6.74e-07 +/- 5.61e-01%

	Nuclide        =	O16
	Cross Sections [cm^-1]:
            Group 1 [0.625      - 20000000.0eV]:	0.00e+00 +/- 0.00e+00%
            Group 2 [0.0        - 0.625     eV]:	0.00e+00 +/- 0.00e+00%



/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1506: RuntimeWarning: invalid value encountered in true_divide
  data = self.std_dev[indices] / self.mean[indices]

One can export the entire Library to HDF5 with the Library.build_hdf5_store(...) method as follows:

In [32]:
# Store the cross section data in an "mgxs/mgxs.h5" HDF5 binary file
mgxs_lib.build_hdf5_store(filename='mgxs.h5', directory='mgxs')
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1835: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1836: RuntimeWarning: invalid value encountered in true_divide
  other_rel_err = data['other']['std. dev.'] / data['other']['mean']
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1837: RuntimeWarning: invalid value encountered in true_divide
  new_tally._mean = data['self']['mean'] / data['other']['mean']

The HDF5 store will contain the numerical multi-group cross section data indexed by domain, nuclide and cross section type. Some data workflows may be optimized by storing and retrieving binary representations of the MGXS objects in the Library. This feature is supported through the Library.dump_to_file(...) and Library.load_from_file(...) routines which use Python's pickle module. This is illustrated as follows.

In [33]:
# Store a Library and its MGXS objects in a pickled binary file "mgxs/mgxs.pkl"
mgxs_lib.dump_to_file(filename='mgxs', directory='mgxs')
In [34]:
# Instantiate a new MGXS Library from the pickled binary file "mgxs/mgxs.pkl"
mgxs_lib = openmc.mgxs.Library.load_from_file(filename='mgxs', directory='mgxs')

The Library class may be used to leverage the energy condensation features supported by the MGXS class. In particular, one can use the Library.get_condensed_library(...) with a coarse group structure which is a subset of the original "fine" group structure as shown below.

In [35]:
# Create a 1-group structure
coarse_groups = openmc.mgxs.EnergyGroups(group_edges=[0., 20.0e6])

# Create a new MGXS Library on the coarse 1-group structure
coarse_mgxs_lib = mgxs_lib.get_condensed_library(coarse_groups)
In [36]:
# Retrieve the NuFissionXS object for the fuel cell from the 1-group library
coarse_fuel_mgxs = coarse_mgxs_lib.get_mgxs(fuel_cell, 'nu-fission')

# Show the Pandas DataFrame for the 1-group MGXS
coarse_fuel_mgxs.get_pandas_dataframe()
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1835: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
Out[36]:
cell group in nuclide mean std. dev.
0 10000 1 U235 0.074393 0.000308
1 10000 1 U238 0.005982 0.000036
2 10000 1 O16 0.000000 0.000000

Verification with OpenMOC

Of course it is always a good idea to verify that one's cross sections are accurate. We can easily do so here with the deterministic transport code OpenMOC. We first construct an equivalent OpenMOC geometry.

In [37]:
# Create an OpenMOC Geometry from the OpenMC Geometry
openmoc_geometry = get_openmoc_geometry(mgxs_lib.geometry)

Now, we can inject the multi-group cross sections into the equivalent fuel assembly OpenMOC geometry. The openmoc.materialize module supports the loading of Library objects from OpenMC as illustrated below.

In [38]:
# Load the library into the OpenMOC geometry
materials = load_openmc_mgxs_lib(mgxs_lib, openmoc_geometry)
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1835: RuntimeWarning: invalid value encountered in true_divide
  self_rel_err = data['self']['std. dev.'] / data['self']['mean']
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1836: RuntimeWarning: invalid value encountered in true_divide
  other_rel_err = data['other']['std. dev.'] / data['other']['mean']
/home/wbinventor/Documents/NSE-CRPG-Codes/openmc/openmc/tallies.py:1837: RuntimeWarning: invalid value encountered in true_divide
  new_tally._mean = data['self']['mean'] / data['other']['mean']

We are now ready to run OpenMOC to verify our cross-sections from OpenMC.

In [39]:
# Generate tracks for OpenMOC
track_generator = openmoc.TrackGenerator(openmoc_geometry, num_azim=32, azim_spacing=0.1)
track_generator.generateTracks()

# Run OpenMOC
solver = openmoc.CPUSolver(track_generator)
solver.computeEigenvalue()
[  NORMAL ]  Importing ray tracing data from file...
[  NORMAL ]  Computing the eigenvalue...
[  NORMAL ]  Iteration 0:	k_eff = 0.823793	res = 0.000E+00
[  NORMAL ]  Iteration 1:	k_eff = 0.780554	res = 1.938E-01
[  NORMAL ]  Iteration 2:	k_eff = 0.739678	res = 6.539E-02
[  NORMAL ]  Iteration 3:	k_eff = 0.711003	res = 5.285E-02
[  NORMAL ]  Iteration 4:	k_eff = 0.689738	res = 3.931E-02
[  NORMAL ]  Iteration 5:	k_eff = 0.675038	res = 3.014E-02
[  NORMAL ]  Iteration 6:	k_eff = 0.665753	res = 2.147E-02
[  NORMAL ]  Iteration 7:	k_eff = 0.661013	res = 1.389E-02
[  NORMAL ]  Iteration 8:	k_eff = 0.660052	res = 7.293E-03
[  NORMAL ]  Iteration 9:	k_eff = 0.662216	res = 2.057E-03
[  NORMAL ]  Iteration 10:	k_eff = 0.666942	res = 3.576E-03
[  NORMAL ]  Iteration 11:	k_eff = 0.673743	res = 7.278E-03
[  NORMAL ]  Iteration 12:	k_eff = 0.682202	res = 1.030E-02
[  NORMAL ]  Iteration 13:	k_eff = 0.691961	res = 1.264E-02
[  NORMAL ]  Iteration 14:	k_eff = 0.702715	res = 1.438E-02
[  NORMAL ]  Iteration 15:	k_eff = 0.714203	res = 1.561E-02
[  NORMAL ]  Iteration 16:	k_eff = 0.726205	res = 1.641E-02
[  NORMAL ]  Iteration 17:	k_eff = 0.738532	res = 1.686E-02
[  NORMAL ]  Iteration 18:	k_eff = 0.751030	res = 1.703E-02
[  NORMAL ]  Iteration 19:	k_eff = 0.763567	res = 1.697E-02
[  NORMAL ]  Iteration 20:	k_eff = 0.776034	res = 1.674E-02
[  NORMAL ]  Iteration 21:	k_eff = 0.788344	res = 1.637E-02
[  NORMAL ]  Iteration 22:	k_eff = 0.800423	res = 1.591E-02
[  NORMAL ]  Iteration 23:	k_eff = 0.812215	res = 1.536E-02
[  NORMAL ]  Iteration 24:	k_eff = 0.823673	res = 1.477E-02
[  NORMAL ]  Iteration 25:	k_eff = 0.834764	res = 1.415E-02
[  NORMAL ]  Iteration 26:	k_eff = 0.845462	res = 1.350E-02
[  NORMAL ]  Iteration 27:	k_eff = 0.855748	res = 1.285E-02
[  NORMAL ]  Iteration 28:	k_eff = 0.865611	res = 1.220E-02
[  NORMAL ]  Iteration 29:	k_eff = 0.875045	res = 1.156E-02
[  NORMAL ]  Iteration 30:	k_eff = 0.884048	res = 1.093E-02
[  NORMAL ]  Iteration 31:	k_eff = 0.892623	res = 1.032E-02
[  NORMAL ]  Iteration 32:	k_eff = 0.900775	res = 9.727E-03
[  NORMAL ]  Iteration 33:	k_eff = 0.908511	res = 9.158E-03
[  NORMAL ]  Iteration 34:	k_eff = 0.915841	res = 8.613E-03
[  NORMAL ]  Iteration 35:	k_eff = 0.922778	res = 8.092E-03
[  NORMAL ]  Iteration 36:	k_eff = 0.929332	res = 7.596E-03
[  NORMAL ]  Iteration 37:	k_eff = 0.935519	res = 7.124E-03
[  NORMAL ]  Iteration 38:	k_eff = 0.941351	res = 6.676E-03
[  NORMAL ]  Iteration 39:	k_eff = 0.946843	res = 6.253E-03
[  NORMAL ]  Iteration 40:	k_eff = 0.952011	res = 5.852E-03
[  NORMAL ]  Iteration 41:	k_eff = 0.956869	res = 5.474E-03
[  NORMAL ]  Iteration 42:	k_eff = 0.961431	res = 5.118E-03
[  NORMAL ]  Iteration 43:	k_eff = 0.965712	res = 4.783E-03
[  NORMAL ]  Iteration 44:	k_eff = 0.969727	res = 4.467E-03
[  NORMAL ]  Iteration 45:	k_eff = 0.973489	res = 4.170E-03
[  NORMAL ]  Iteration 46:	k_eff = 0.977013	res = 3.892E-03
[  NORMAL ]  Iteration 47:	k_eff = 0.980310	res = 3.631E-03
[  NORMAL ]  Iteration 48:	k_eff = 0.983394	res = 3.386E-03
[  NORMAL ]  Iteration 49:	k_eff = 0.986277	res = 3.156E-03
[  NORMAL ]  Iteration 50:	k_eff = 0.988971	res = 2.942E-03
[  NORMAL ]  Iteration 51:	k_eff = 0.991487	res = 2.740E-03
[  NORMAL ]  Iteration 52:	k_eff = 0.993835	res = 2.552E-03
[  NORMAL ]  Iteration 53:	k_eff = 0.996026	res = 2.376E-03
[  NORMAL ]  Iteration 54:	k_eff = 0.998069	res = 2.212E-03
[  NORMAL ]  Iteration 55:	k_eff = 0.999974	res = 2.059E-03
[  NORMAL ]  Iteration 56:	k_eff = 1.001749	res = 1.915E-03
[  NORMAL ]  Iteration 57:	k_eff = 1.003403	res = 1.782E-03
[  NORMAL ]  Iteration 58:	k_eff = 1.004943	res = 1.657E-03
[  NORMAL ]  Iteration 59:	k_eff = 1.006377	res = 1.540E-03
[  NORMAL ]  Iteration 60:	k_eff = 1.007711	res = 1.432E-03
[  NORMAL ]  Iteration 61:	k_eff = 1.008952	res = 1.331E-03
[  NORMAL ]  Iteration 62:	k_eff = 1.010107	res = 1.237E-03
[  NORMAL ]  Iteration 63:	k_eff = 1.011181	res = 1.149E-03
[  NORMAL ]  Iteration 64:	k_eff = 1.012179	res = 1.067E-03
[  NORMAL ]  Iteration 65:	k_eff = 1.013107	res = 9.909E-04
[  NORMAL ]  Iteration 66:	k_eff = 1.013969	res = 9.201E-04
[  NORMAL ]  Iteration 67:	k_eff = 1.014769	res = 8.542E-04
[  NORMAL ]  Iteration 68:	k_eff = 1.015513	res = 7.929E-04
[  NORMAL ]  Iteration 69:	k_eff = 1.016204	res = 7.358E-04
[  NORMAL ]  Iteration 70:	k_eff = 1.016845	res = 6.828E-04
[  NORMAL ]  Iteration 71:	k_eff = 1.017440	res = 6.335E-04
[  NORMAL ]  Iteration 72:	k_eff = 1.017993	res = 5.877E-04
[  NORMAL ]  Iteration 73:	k_eff = 1.018505	res = 5.451E-04
[  NORMAL ]  Iteration 74:	k_eff = 1.018980	res = 5.056E-04
[  NORMAL ]  Iteration 75:	k_eff = 1.019422	res = 4.688E-04
[  NORMAL ]  Iteration 76:	k_eff = 1.019831	res = 4.347E-04
[  NORMAL ]  Iteration 77:	k_eff = 1.020210	res = 4.030E-04
[  NORMAL ]  Iteration 78:	k_eff = 1.020562	res = 3.736E-04
[  NORMAL ]  Iteration 79:	k_eff = 1.020888	res = 3.463E-04
[  NORMAL ]  Iteration 80:	k_eff = 1.021190	res = 3.209E-04
[  NORMAL ]  Iteration 81:	k_eff = 1.021471	res = 2.974E-04
[  NORMAL ]  Iteration 82:	k_eff = 1.021730	res = 2.756E-04
[  NORMAL ]  Iteration 83:	k_eff = 1.021971	res = 2.553E-04
[  NORMAL ]  Iteration 84:	k_eff = 1.022194	res = 2.365E-04
[  NORMAL ]  Iteration 85:	k_eff = 1.022400	res = 2.191E-04
[  NORMAL ]  Iteration 86:	k_eff = 1.022592	res = 2.030E-04
[  NORMAL ]  Iteration 87:	k_eff = 1.022769	res = 1.880E-04
[  NORMAL ]  Iteration 88:	k_eff = 1.022933	res = 1.741E-04
[  NORMAL ]  Iteration 89:	k_eff = 1.023085	res = 1.612E-04
[  NORMAL ]  Iteration 90:	k_eff = 1.023226	res = 1.493E-04
[  NORMAL ]  Iteration 91:	k_eff = 1.023356	res = 1.382E-04
[  NORMAL ]  Iteration 92:	k_eff = 1.023477	res = 1.279E-04
[  NORMAL ]  Iteration 93:	k_eff = 1.023588	res = 1.184E-04
[  NORMAL ]  Iteration 94:	k_eff = 1.023692	res = 1.097E-04
[  NORMAL ]  Iteration 95:	k_eff = 1.023788	res = 1.015E-04
[  NORMAL ]  Iteration 96:	k_eff = 1.023876	res = 9.392E-05
[  NORMAL ]  Iteration 97:	k_eff = 1.023958	res = 8.694E-05
[  NORMAL ]  Iteration 98:	k_eff = 1.024034	res = 8.045E-05
[  NORMAL ]  Iteration 99:	k_eff = 1.024104	res = 7.446E-05
[  NORMAL ]  Iteration 100:	k_eff = 1.024169	res = 6.890E-05
[  NORMAL ]  Iteration 101:	k_eff = 1.024229	res = 6.374E-05
[  NORMAL ]  Iteration 102:	k_eff = 1.024285	res = 5.897E-05
[  NORMAL ]  Iteration 103:	k_eff = 1.024336	res = 5.457E-05
[  NORMAL ]  Iteration 104:	k_eff = 1.024384	res = 5.049E-05
[  NORMAL ]  Iteration 105:	k_eff = 1.024428	res = 4.672E-05
[  NORMAL ]  Iteration 106:	k_eff = 1.024469	res = 4.321E-05
[  NORMAL ]  Iteration 107:	k_eff = 1.024506	res = 3.994E-05
[  NORMAL ]  Iteration 108:	k_eff = 1.024541	res = 3.696E-05
[  NORMAL ]  Iteration 109:	k_eff = 1.024574	res = 3.418E-05
[  NORMAL ]  Iteration 110:	k_eff = 1.024603	res = 3.164E-05
[  NORMAL ]  Iteration 111:	k_eff = 1.024631	res = 2.925E-05
[  NORMAL ]  Iteration 112:	k_eff = 1.024657	res = 2.705E-05
[  NORMAL ]  Iteration 113:	k_eff = 1.024680	res = 2.501E-05
[  NORMAL ]  Iteration 114:	k_eff = 1.024702	res = 2.313E-05
[  NORMAL ]  Iteration 115:	k_eff = 1.024722	res = 2.140E-05
[  NORMAL ]  Iteration 116:	k_eff = 1.024741	res = 1.981E-05
[  NORMAL ]  Iteration 117:	k_eff = 1.024758	res = 1.827E-05
[  NORMAL ]  Iteration 118:	k_eff = 1.024774	res = 1.690E-05
[  NORMAL ]  Iteration 119:	k_eff = 1.024789	res = 1.564E-05
[  NORMAL ]  Iteration 120:	k_eff = 1.024802	res = 1.445E-05
[  NORMAL ]  Iteration 121:	k_eff = 1.024815	res = 1.338E-05
[  NORMAL ]  Iteration 122:	k_eff = 1.024827	res = 1.236E-05
[  NORMAL ]  Iteration 123:	k_eff = 1.024837	res = 1.142E-05
[  NORMAL ]  Iteration 124:	k_eff = 1.024847	res = 1.057E-05

We report the eigenvalues computed by OpenMC and OpenMOC here together to summarize our results.

In [40]:
# Print report of keff and bias with OpenMC
openmoc_keff = solver.getKeff()
openmc_keff = sp.k_combined[0]
bias = (openmoc_keff - openmc_keff) * 1e5

print('openmc keff = {0:1.6f}'.format(openmc_keff))
print('openmoc keff = {0:1.6f}'.format(openmoc_keff))
print('bias [pcm]: {0:1.1f}'.format(bias))
openmc keff = 1.023307
openmoc keff = 1.024847
bias [pcm]: 154.0

There is a non-trivial bias between the eigenvalues computed by OpenMC and OpenMOC. One can show that these biases do not converge to <100 pcm with more particle histories. For heterogeneous geometries, additional measures must be taken to address the following three sources of bias:

  • Appropriate transport-corrected cross sections
  • Spatial discretization of OpenMOC's mesh
  • Constant-in-angle multi-group cross sections

Flux and Pin Power Visualizations

We will conclude this tutorial by illustrating how to visualize the fission rates computed by OpenMOC and OpenMC. First, we extract volume-integrated fission rates from OpenMC's mesh fission rate tally for each pin cell in the fuel assembly.

In [41]:
# Get the OpenMC fission rate mesh tally data
mesh_tally = sp.get_tally(name='mesh tally')
openmc_fission_rates = mesh_tally.get_values(scores=['nu-fission'])

# Reshape array to 2D for plotting
openmc_fission_rates.shape = (17,17)

# Normalize to the average pin power
openmc_fission_rates /= np.mean(openmc_fission_rates)

Next, we extract OpenMOC's volume-averaged fission rates into a 2D 17x17 NumPy array.

In [42]:
# Create OpenMOC Mesh on which to tally fission rates
openmoc_mesh = openmoc.process.Mesh()
openmoc_mesh.dimension = np.array(mesh.dimension)
openmoc_mesh.lower_left = np.array(mesh.lower_left)
openmoc_mesh.upper_right = np.array(mesh.upper_right)
openmoc_mesh.width = openmoc_mesh.upper_right - openmoc_mesh.lower_left
openmoc_mesh.width /= openmoc_mesh.dimension

# Tally OpenMOC fission rates on the Mesh
openmoc_fission_rates = openmoc_mesh.tally_fission_rates(solver)
openmoc_fission_rates = np.squeeze(openmoc_fission_rates)
openmoc_fission_rates = np.fliplr(openmoc_fission_rates)

# Normalize to the average pin fission rate
openmoc_fission_rates /= np.mean(openmoc_fission_rates)

Now we can easily use Matplotlib to visualize the fission rates from OpenMC and OpenMOC side-by-side.

In [43]:
# Ignore zero fission rates in guide tubes with Matplotlib color scheme
openmc_fission_rates[openmc_fission_rates == 0] = np.nan
openmoc_fission_rates[openmoc_fission_rates == 0] = np.nan

# Plot OpenMC's fission rates in the left subplot
fig = plt.subplot(121)
plt.imshow(openmc_fission_rates, interpolation='none', cmap='jet')
plt.title('OpenMC Fission Rates')

# Plot OpenMOC's fission rates in the right subplot
fig2 = plt.subplot(122)
plt.imshow(openmoc_fission_rates, interpolation='none', cmap='jet')
plt.title('OpenMOC Fission Rates')
Out[43]:
<matplotlib.text.Text at 0x7f3a95649710>