Initial TRISO modeling capabilities

docs/source/pythonapi/index.rst

@@ -296,6 +296,33 @@ Multi-group Cross Section Libraries
 
     openmc.mgxs.Library
 
+-------------------------------------
+:mod:`openmc.model` -- Model Building
+-------------------------------------
+
+TRISO Fuel Modeling
+-------------------
+
+Classes
++++++++
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myclass.rst
+
+   openmc.model.TRISO
+
+Functions
++++++++++
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+
+   openmc.model.create_triso_lattice
+
+
 .. _Jupyter: https://jupyter.org/
 .. _NumPy: http://www.numpy.org/
 .. _Codecademy: https://www.codecademy.com/tracks/python

openmc/model/__init__.py

@@ -0,0 +1 @@
+from .triso import *

openmc/model/triso.py

@@ -0,0 +1,172 @@
+import copy
+from collections import Iterable
+from numbers import Real
+
+import numpy as np
+
+import openmc
+import openmc.checkvalue as cv
+
+class TRISO(object):
+    """Tristructural-isotopic (TRISO) micro fuel particle
+
+    Parameters
+    ----------
+    materials : Iterable of openmc.Material
+        Material to be assigned to each layer of the TRISO particle starting
+        with the innermost and proceeding outwards
+    radii : Iterable of float
+        Outer radii in cm of each layer of the TRISO particle in ascending order
+    center : Iterable of float
+        Cartesian coordinates of the center of the TRISO particle in cm
+
+    Attributes
+    ----------
+    cells : list of opemc.Cell
+        Each layer of the TRISO particle
+    center : numpy.ndarray
+        Cartesian coordinates of the center of the TRISO particle in cm
+    outside : openmc.Region
+        Region of space outside of the TRISO particle
+    bounding_box : tuple of numpy.ndarray
+        Lower-left and upper-right coordinates of an axis-aligned bounding box
+        for the TRISO particle
+
+    """
+
+    def __init__(self, materials, radii, center=(0., 0., 0.)):
+        surfaces = [openmc.Sphere(R=r) for r in radii]
+        cells = []
+        for i, m in enumerate(materials):
+            c = openmc.Cell(fill=m)
+            if i == 0:
+                c.region = -surfaces[i]
+            else:
+                c.region = +surfaces[i-1] & -surfaces[i]
+            cells.append(c)
+        self._cells = cells
+        self._surfaces = surfaces
+        self.center = np.asarray(center)
+
+    @property
+    def bounding_box(self):
+        return self.cells[-1].region.bounding_box
+
+    @property
+    def cells(self):
+        return self._cells
+
+    @property
+    def center(self):
+        return self._center
+
+    @property
+    def outside(self):
+        return ~self.cells[-1].region.nodes[-1]
+
+    @center.setter
+    def center(self, center):
+        cv.check_type('TRISO center', center, Iterable, Real)
+        for s in self._surfaces:
+            s.x0, s.y0, s.z0 = center
+        self._center = center
+
+    def classify(self, lattice):
+        """Determine lattice element indices which might contain the TRISO particle.
+
+        Parameters
+        ----------
+        lattice : openmc.RectLattice
+            Lattice to check
+
+        Returns
+        -------
+        list of tuple
+            (z,y,x) lattice element indices which might contain the TRISO
+            particle.
+
+        """
+
+        ll, ur = self.bounding_box
+        if lattice.ndim == 2:
+            (i_min, j_min), p = lattice.find_element(ll)
+            (i_max, j_max), p = lattice.find_element(ur)
+            return list(np.broadcast(*np.ogrid[
+                j_min:j_max+1, i_min:i_max+1]))
+        else:
+            (i_min, j_min, k_min), p = lattice.find_element(ll)
+            (i_max, j_max, k_max), p = lattice.find_element(ur)
+            return list(np.broadcast(*np.ogrid[
+                k_min:k_max+1, j_min:j_max+1, i_min:i_max+1]))
+
+
+def create_triso_lattice(trisos, lower_left, pitch, shape, background):
+    """Create a lattice containing TRISO particles for optimized tracking.
+
+    Parameters
+    ----------
+    trisos : list of openmc.model.TRISO
+        List of TRISO particles to put in lattice
+    lower_left : Iterable of float
+        Lower-left Cartesian coordinates of the lattice
+    pitch : Iterable of float
+        Pitch of the lattice elements in the x-, y-, and z-directions
+    shape : Iterable of float
+        Number of lattice elements in the x-, y-, and z-directions
+    background : openmc.Material
+        A background material that is used anywhere within the lattice but
+        outside a TRISO particle
+
+    Returns
+    -------
+    lattice : openmc.RectLattice
+        A lattice containing the TRISO particles
+
+    """
+
+    lattice = openmc.RectLattice()
+    lattice.lower_left = lower_left
+    lattice.pitch = pitch
+
+    indices = list(np.broadcast(*np.ogrid[:shape[2], :shape[1], :shape[0]]))
+    triso_locations = {idx: [] for idx in indices}
+    for t in trisos:
+        for idx in t.classify(lattice):
+            if idx in sorted(triso_locations):
+                # Create copy of TRISO particle with materials preserved and
+                # different cell/surface IDs
+                t_copy = copy.deepcopy(t)
+                for c, c_copy in zip(t.cells, t_copy.cells):
+                    c_copy.id = None
+                    c_copy.fill = c.fill
+                for s in t_copy._surfaces:
+                    s.id = None
+                triso_locations[idx].append(t_copy)
+
+    # Create universes
+    universes = np.empty(shape[::-1], dtype=openmc.Universe)
+    for idx, triso_list in sorted(triso_locations.items()):
+        if len(triso_list) > 0:
+            outside_trisos = openmc.Intersection(*[t.outside for t in triso_list])
+            background_cell = openmc.Cell(fill=background, region=outside_trisos)
+        else:
+            background_cell = openmc.Cell(fill=background)
+
+        u = openmc.Universe()
+        u.add_cell(background_cell)
+        for t in triso_list:
+            u.add_cells(t.cells)
+            iz, iy, ix = idx
+            t.center = lattice.get_local_coordinates(t.center, (ix, iy, iz))
+
+        if len(shape) == 2:
+            universes[-1 - idx[0], idx[1]] = u
+        else:
+            universes[idx[0], -1 - idx[1], idx[2]] = u
+    lattice.universes = universes
+
+    # Set outer universe
+    background_cell = openmc.Cell(fill=background)
+    lattice.outer = openmc.Universe(cells=[background_cell])
+
+    return lattice

setup.py

@@ -11,7 +11,8 @@ except ImportError:
 
 kwargs = {'name': 'openmc',
           'version': '0.7.1',
-          'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.stats'],
+          'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.model',
+                       'openmc.stats'],
           'scripts': glob.glob('scripts/openmc-*'),
 
           # Metadata

tests/test_triso/inputs_true.dat

@@ -0,0 +1 @@
+6c6fecedf3db0b91b69b7b7b57b3a52c42947253bea4ee3889d6bbd6e74935cc4e599c1c743eb589a1045868412f3f88c5fef7cf10e180bdc4bd182970c9ee15

tests/test_triso/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.685303E+00 1.121936E-01

tests/test_triso/test_triso.py

@@ -0,0 +1,107 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import glob
+import random
+from math import sqrt
+
+import numpy as np
+
+sys.path.insert(0, os.pardir)
+from testing_harness import PyAPITestHarness
+import openmc
+import openmc.model
+
+
+class TRISOTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+        # Define TRISO matrials
+        fuel = openmc.Material()
+        fuel.set_density('g/cm3', 10.5)
+        fuel.add_nuclide('U-235', 0.14154)
+        fuel.add_nuclide('U-238', 0.85846)
+        fuel.add_nuclide('C-Nat', 0.5)
+        fuel.add_nuclide('O-16', 1.5)
+
+        porous_carbon = openmc.Material()
+        porous_carbon.set_density('g/cm3', 1.0)
+        porous_carbon.add_nuclide('C-Nat', 1.0)
+        porous_carbon.add_s_alpha_beta('Graph', '71t')
+
+        ipyc = openmc.Material()
+        ipyc.set_density('g/cm3', 1.90)
+        ipyc.add_nuclide('C-Nat', 1.0)
+        ipyc.add_s_alpha_beta('Graph', '71t')
+
+        sic = openmc.Material()
+        sic.set_density('g/cm3', 3.20)
+        sic.add_element('Si', 1.0)
+        sic.add_nuclide('C-Nat', 1.0)
+
+        opyc = openmc.Material()
+        opyc.set_density('g/cm3', 1.87)
+        opyc.add_nuclide('C-Nat', 1.0)
+        opyc.add_s_alpha_beta('Graph', '71t')
+
+        graphite = openmc.Material()
+        graphite.set_density('g/cm3', 1.1995)
+        graphite.add_nuclide('C-Nat', 1.0)
+        graphite.add_s_alpha_beta('Graph', '71t')
+
+        # Create TRISO particles
+        materials = [fuel, porous_carbon, ipyc, sic, opyc]
+        radii = np.array([212.5, 312.5, 347.5, 382.5, 422.5])*1e-4
+        trisos = []
+        random.seed(1)
+        for i in range(100):
+            # Randomly sample location
+            x = random.uniform(-0.5, 0.5)
+            y = random.uniform(-0.5, 0.5)
+            z = random.uniform(-0.5, 0.5)
+            t = openmc.model.TRISO(materials, radii, (x, y, z))
+
+            # Make sure TRISO doesn't overlap with another
+            for tp in trisos:
+                xp, yp, zp = tp.center
+                distance = sqrt((x - xp)**2 + (y - yp)**2 + (z - zp)**2)
+                if distance <= 2*radii[-1]:
+                    break
+            else:
+                trisos.append(t)
+
+        # Define box to contain lattice
+        min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')
+        max_x = openmc.XPlane(x0=0.5, boundary_type='reflective')
+        min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective')
+        max_y = openmc.YPlane(y0=0.5, boundary_type='reflective')
+        min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective')
+        max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective')
+        box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z)
+
+        # Create lattice
+        ll, ur = box.region.bounding_box
+        shape = (3, 3, 3)
+        lattice = openmc.model.create_triso_lattice(
+            trisos, ll, (ur - ll)/shape, shape, graphite)
+        box.fill = lattice
+
+        root = openmc.Universe(0, cells=[box])
+        geom = openmc.Geometry(root)
+        geom.export_to_xml()
+
+        settings = openmc.Settings()
+        settings.batches = 5
+        settings.inactive = 0
+        settings.particles = 50
+        settings.source = openmc.Source(space=openmc.stats.Point())
+        settings.export_to_xml()
+
+        mats = openmc.Materials([fuel, porous_carbon, ipyc, sic, opyc, graphite])
+        mats.default_xs = '71c'
+        mats.export_to_xml()
+
+
+if __name__ == '__main__':
+    harness = TRISOTestHarness('statepoint.5.h5')
+    harness.main()