Unstructured Mesh Tallies in OpenMC

In this example we'll look at how to setup and use unstructured mesh tallies in OpenMC. Unstructured meshes are able to provide results over spatial regions of a problem while conforming to a specific geometric features -- something that is often difficult to do using the regular and rectilinear meshes in OpenMC.

Here, we'll apply an unstructured mesh tally to the PWR assembly model from the OpenMC examples.

NOTE: This notebook will not run successfully if OpenMC has not been built with DAGMC support enabled.


In [1]:
from IPython.display import Image
import openmc
import openmc.lib

assert(openmc.lib._dagmc_enabled())

We'll need to download the unstructured mesh file used in this notebook. We'll be retrieving those using the function and URLs below.


In [2]:
%matplotlib inline
from matplotlib import pyplot as plt
plt.rcParams["figure.figsize"] = (30,10)

import urllib.request

pin_mesh_url = 'https://tinyurl.com/u9ce9d7' # 1.2 MB

def download(url, filename='dagmc.h5m'):
    """
    Helper function for retrieving dagmc models
    """
    u = urllib.request.urlopen(url)
    
    if u.status != 200:
        raise RuntimeError("Failed to download file.")
    
    # save file as dagmc.h5m
    with open(filename, 'wb') as f:
        f.write(u.read())

First we'll import that model from the set of OpenMC examples.


In [3]:
model = openmc.examples.pwr_assembly()

We'll make a couple of adjustments to this 2D model as it won't play very well with the 3D mesh we'll be looking at. First, we'll bound the pincell between +/- 10 cm in the Z dimension.


In [4]:
min_z = openmc.ZPlane(z0=-10.0)
max_z = openmc.ZPlane(z0=10.0)

z_region = +min_z & -max_z

cells = model.geometry.get_all_cells()
for cell in cells.values():
    cell.region &= z_region

The other adjustment we'll make is to remove the reflective boundary conditions on the X and Y boundaries. (This is purely to generate a more interesting flux profile.)


In [5]:
surfaces = model.geometry.get_all_surfaces()
# modify the boundary condition of the
# planar surfaces bounding the assembly
for surface in surfaces.values():
    if isinstance(surface, openmc.Plane):
        surface.boundary_type = 'vacuum'

Let's take a quick look at the model to ensure our changs have been added properly.


In [6]:
root_univ = model.geometry.root_universe

# axial image
root_univ.plot(width=(22.0, 22.0),
               pixels=(200, 300),
               basis='xz',
               color_by='material',
               seed=0)


Out[6]:
<matplotlib.image.AxesImage at 0x7fe9650cd320>

In [7]:
# radial image
root_univ.plot(width=(22.0, 22.0),
               pixels=(400, 400),
               basis='xy',
               color_by='material',
               seed=0)


Out[7]:
<matplotlib.image.AxesImage at 0x7fe967252e48>

Looks good! Let's run some particles through the problem.


In [8]:
model.run()


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

                   | The OpenMC Monte Carlo Code
         Copyright | 2011-2020 MIT and OpenMC contributors
           License | http://openmc.readthedocs.io/en/latest/license.html
           Version | 0.12.0-dev
          Git SHA1 | e7eceff2aa3a9bbfd4dd553f585d1a31b051ddb3
         Date/Time | 2020-03-27 11:18:32
    OpenMP Threads | 2

 Reading settings XML file...
 Reading cross sections XML file...
 Reading materials XML file...
 Reading geometry XML file...
 Reading U234 from /home/shriwise/opt/openmc/xs/nndc_hdf5/U234.h5
 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 Zr90 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr90.h5
 Reading Zr91 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr91.h5
 Reading Zr92 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr92.h5
 Reading Zr94 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr94.h5
 Reading Zr96 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr96.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 B11 from /home/shriwise/opt/openmc/xs/nndc_hdf5/B11.h5
 Reading c_H_in_H2O from /home/shriwise/opt/openmc/xs/nndc_hdf5/c_H_in_H2O.h5
 Maximum neutron transport energy: 20000000.000000 eV for U235
 Minimum neutron data temperature: 294.000000 K
 Maximum neutron data temperature: 294.000000 K
 Preparing distributed cell instances...
 Writing summary.h5 file...
 Initializing source particles...

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

  Bat./Gen.      k            Average k
  =========   ========   ====================
        1/1    0.20444
        2/1    0.15502
        3/1    0.19804
        4/1    0.22159
        5/1    0.19776
        6/1    0.20086
        7/1    0.21896    0.20991 +/- 0.00905
        8/1    0.23134    0.21706 +/- 0.00885
        9/1    0.29029    0.23536 +/- 0.01935
       10/1    0.20094    0.22848 +/- 0.01649
 Creating state point statepoint.10.h5...

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

 Total time for initialization     = 6.8458e-01 seconds
   Reading cross sections          = 6.7039e-01 seconds
 Total time in simulation          = 1.9462e-02 seconds
   Time in transport only          = 1.7538e-02 seconds
   Time in inactive batches        = 9.2816e-03 seconds
   Time in active batches          = 1.0181e-02 seconds
   Time synchronizing fission bank = 5.2594e-05 seconds
     Sampling source sites         = 4.6704e-05 seconds
     SEND/RECV source sites        = 3.1580e-06 seconds
   Time accumulating tallies       = 2.0290e-06 seconds
 Total time for finalization       = 1.0320e-06 seconds
 Total time elapsed                = 7.0439e-01 seconds
 Calculation Rate (inactive)       = 53869.9 particles/second
 Calculation Rate (active)         = 49113.3 particles/second

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

 k-effective (Collision)     = 0.25094 +/- 0.02010
 k-effective (Track-length)  = 0.22848 +/- 0.01649
 k-effective (Absorption)    = 0.21556 +/- 0.04156
 Combined k-effective        = 0.20707 +/- 0.01965
 Leakage Fraction            = 0.79200 +/- 0.03382

Out[8]:
0.20706788967181938+/-0.019653037754288005

Now it's time to apply our mesh tally to the problem. We'll be using the tetrahedral mesh "pins1-4.h5m" shown below:


In [9]:
Image("./images/pin_mesh.png", width=600)


Out[9]:

This mesh was generated using Trelis with radii that match the fuel/coolant channels of the PWR model. These four channels correspond to the highlighted channels of the assembly below.

Two of the channels are coolant and the other two are fuel.


In [10]:
from matplotlib.patches import Rectangle
from matplotlib import pyplot as plt

pitch = 1.26 # cm

img = root_univ.plot(width=(22.0, 22.0),
                           pixels=(600, 600),
                           basis='xy',
                           color_by='material',
                           seed=0)

# highlight channels
for i in range(0, 4):
    corner = (i * pitch - pitch / 2.0, -i * pitch - pitch / 2.0)
    rect = Rectangle(corner,
                     pitch,
                     pitch,
                     edgecolor='blue',
                     fill=False)
    img.axes.add_artist(rect)


Applying an unstructured mesh tally

To use this mesh, we'll create an unstructured mesh instance and apply it to a mesh filter.


In [11]:
download(pin_mesh_url, "pins1-4.h5m")
umesh = openmc.UnstructuredMesh("pins1-4.h5m")
mesh_filter = openmc.MeshFilter(umesh)

We can now apply this filter like any other. For this demonstration we'll score both the flux and heating in these pins.


In [12]:
tally = openmc.Tally()
tally.filters = [mesh_filter]
tally.scores = ['heating', 'flux']
tally.estimator = 'tracklength'
model.tallies = [tally]

Now we'll run this model with the unstructured mesh tally applied. Notice that the simulation takes some time to start due to some additional data structures used by the unstructured mesh tally. Additionally, the particle rate drops dramatically during the active cycles of this simulation.

Unstructured meshes are useful, but they can be computationally expensive!


In [13]:
model.settings.particles = 100_000
model.settings.inactive = 20
model.settings.batches = 100
model.run(output=False)


Out[13]:
0.2311872124292096+/-0.00017018932313754055

At the end of the simulation, we see the statepoint file along with a file named "tally_1.100.vtk". This file contains the results of the unstructured mesh tally with convenient labels for the scores applied. In our case the following scores will be present in the VTK:

  • flux_total_value
  • flux_total_std_dev
  • heating_total_value
  • heading_total_std_dev

    Where "total" represents

Currently, an unstructured VTK file will only be generated for tallies if the unstructured mesh is is the only filter applied to that tally. All results for the unstructured mesh tally are present in the statepoint file regardless of the number of filters applied, however.

These files can be viewed using free tools like Paraview and VisIt to examine the results.


In [14]:
!ls *.vtk


tally_1.100.vtk

Flux


In [15]:
Image("./images/umesh_flux.png", width=600)


Out[15]:

Heating

Here is an image of the heating score as viewed in VisIt. Note that no heating is scored in the water-filled channels as expected.


In [16]:
Image("./images/umesh_heating.png", width=600)


Out[16]:

Statepoint Data


In [17]:
with openmc.StatePoint("statepoint.100.h5") as sp:
    tally = sp.tallies[1]
    umesh = sp.meshes[1]

Enough information for visualization of results on the unstructured mesh is also provided in the statepoint file. Namely, the mesh element volumes and centroids are available.


In [18]:
print(umesh.volumes)
print(umesh.centroids)


[1.43381086e-04 1.48043747e-04 1.60408339e-04 ... 7.04197023e-05
 7.04197023e-05 7.04197023e-05]
[[ 2.88485691 -2.55429784  9.97768184]
 [ 2.87565092 -2.60469781  9.8884092 ]
 [ 2.85832254 -2.65291228  9.97768184]
 ...
 [ 1.46082175 -1.15569203 -3.62914358]
 [ 1.4443143  -1.1321793  -3.65475081]
 [ 1.46884412 -1.15657736 -3.68206543]]

The combination of these values can provide for an appoxmiate visualization of the unstructured mesh without its explicit representation or use of an additional mesh library.

We hope you've found this example notebook useful!


In [19]:
Image("./images/umesh_w_assembly.png", width=600)


Out[19]: