6. Defining Geometry

6.1. Surfaces and Regions

The geometry of a model in OpenMC is defined using constructive solid geometry (CSG), also sometimes referred to as combinatorial geometry. CSG allows a user to create complex regions using Boolean operators (intersection, union, and complement) on simpler regions. In order to define a region that we can assign to a cell, we must first define surfaces which bound the region. A surface is a locus of zeros of a function of Cartesian coordinates \(x,y,z\), e.g.

  • A plane perpendicular to the \(x\) axis: \(x - x_0 = 0\)
  • A cylinder parallel to the \(z\) axis: \((x - x_0)^2 + (y - y_0)^2 - R^2 = 0\)
  • A sphere: \((x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\)

Defining a surface alone is not sufficient to specify a volume – in order to define an actual volume, one must reference the half-space of a surface. A surface half-space is the region whose points satisfy a positive or negative inequality of the surface equation. For example, for a sphere of radius one centered at the origin, the surface equation is \(f(x,y,z) = x^2 + y^2 + z^2 - 1 = 0\). Thus, we say that the negative half-space of the sphere, is defined as the collection of points satisfying \(f(x,y,z) < 0\), which one can reason is the inside of the sphere. Conversely, the positive half-space of the sphere would correspond to all points outside of the sphere, satisfying \(f(x,y,z) > 0\).

In the Python API, surfaces are created via subclasses of openmc.Surface. The available surface types and their corresponding classes are listed in the following table.

Surface types available in OpenMC.
Surface Equation Class
Plane perpendicular to \(x\)-axis \(x - x_0 = 0\) openmc.XPlane
Plane perpendicular to \(y\)-axis \(y - y_0 = 0\) openmc.YPlane
Plane perpendicular to \(z\)-axis \(z - z_0 = 0\) openmc.ZPlane
Arbitrary plane \(Ax + By + Cz = D\) openmc.Plane
Infinite cylinder parallel to \(x\)-axis \((y-y_0)^2 + (z-z_0)^2 - R^2 = 0\) openmc.XCylinder
Infinite cylinder parallel to \(y\)-axis \((x-x_0)^2 + (z-z_0)^2 - R^2 = 0\) openmc.YCylinder
Infinite cylinder parallel to \(z\)-axis \((x-x_0)^2 + (y-y_0)^2 - R^2 = 0\) openmc.ZCylinder
Sphere \((x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 - R^2 = 0\) openmc.Sphere
Cone parallel to the \(x\)-axis \((y-y_0)^2 + (z-z_0)^2 - R^2(x-x_0)^2 = 0\) openmc.XCone
Cone parallel to the \(y\)-axis \((x-x_0)^2 + (z-z_0)^2 - R^2(y-y_0)^2 = 0\) openmc.YCone
Cone parallel to the \(z\)-axis \((x-x_0)^2 + (y-y_0)^2 - R^2(z-z_0)^2 = 0\) openmc.ZCone
General quadric surface \(Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0\) openmc.Quadric

Each surface is characterized by several parameters. As one example, the parameters for a sphere are the \(x,y,z\) coordinates of the center of the sphere and the radius of the sphere. All of these parameters can be set either as optional keyword arguments to the class constructor or via attributes:

sphere = openmc.Sphere(r=10.0)

# This is equivalent
sphere = openmc.Sphere()
sphere.r = 10.0

Once a surface has been created, half-spaces can be obtained by applying the unary - or + operators, corresponding to the negative and positive half-spaces, respectively. For example:

>>> sphere = openmc.Sphere(r=10.0)
>>> inside_sphere = -sphere
>>> outside_sphere = +sphere
>>> type(inside_sphere)
<class 'openmc.surface.Halfspace'>

Instances of openmc.Halfspace can be combined together using the Boolean operators & (intersection), | (union), and ~ (complement):

>>> inside_sphere = -openmc.Sphere()
>>> above_plane = +openmc.ZPlane()
>>> northern_hemisphere = inside_sphere & above_plane
>>> type(northern_hemisphere)
<class 'openmc.region.Intersection'>

For many regions, a bounding-box can be determined automatically:

>>> northern_hemisphere.bounding_box
(array([-1., -1., 0.]), array([1., 1., 1.]))

While a bounding box can be determined for regions involving half-spaces of spheres, cylinders, and axis-aligned planes, it generally cannot be determined if the region involves cones, non-axis-aligned planes, or other exotic second-order surfaces. For example, the openmc.model.hexagonal_prism() function returns the interior region of a hexagonal prism; because it is bounded by a openmc.Plane, trying to get its bounding box won’t work:

>>> hex = openmc.model.hexagonal_prism()
>>> hex.bounding_box
(array([-0.8660254,       -inf,       -inf]),
 array([ 0.8660254,        inf,        inf]))

6.1.1. Boundary Conditions

When a surface is created, by default particles that pass through the surface will consider it to be transmissive, i.e., they pass through the surface freely. If your model does not extend to infinity in all spatial dimensions, you may want to specify different behavior for particles passing through a surface. To specify a vacuum boundary condition, simply change the Surface.boundary_type attribute to ‘vacuum’:

outer_surface = openmc.Sphere(r=100.0, boundary_type='vacuum')

# This is equivalent
outer_surface = openmc.Sphere(r=100.0)
outer_surface.boundary_type = 'vacuum'

Reflective and periodic boundary conditions can be set with the strings ‘reflective’ and ‘periodic’. Vacuum and reflective boundary conditions can be applied to any type of surface. Periodic boundary conditions can be applied to pairs of planar surfaces. For axis-aligned planes, matching periodic surfaces can be determined automatically. For non-axis-aligned planes, it is necessary to specify pairs explicitly using the Surface.periodic_surface attribute as in the following example:

p1 = openmc.Plane(a=0.3, b=5.0, d=1.0, boundary_type='periodic')
p2 = openmc.Plane(a=0.3, b=5.0, d=-1.0, boundary_type='periodic')
p1.periodic_surface = p2

Rotationally-periodic boundary conditions can be specified for a pair of XPlane and YPlane; in that case, the Surface.periodic_surface attribute must be specified manually as well.

Caution

When using rotationally-periodic boundary conditions, your geometry must be defined in the first quadrant, i.e., above the y-plane and to the right of the x-plane.

6.2. Cells

Once you have a material created and a region of space defined, you need to define a cell that assigns the material to the region. Cells are created using the openmc.Cell class:

fuel = openmc.Cell(fill=uo2, region=pellet)

# This is equivalent
fuel = openmc.Cell()
fuel.fill = uo2
fuel.region = pellet

In this example, an instance of openmc.Material is assigned to the Cell.fill attribute. One can also fill a cell with a universe or lattice.

The classes Halfspace, Intersection, Union, and Complement and all instances of openmc.Region and can be assigned to the Cell.region attribute.

6.3. Universes

Similar to MCNP and Serpent, OpenMC is capable of using universes, collections of cells that can be used as repeatable units of geometry. At a minimum, there must be one “root” universe present in the model. To define a universe, an instance of openmc.Universe is created and then cells can be added using the Universe.add_cells() or Universe.add_cell() methods. Alternatively, a list of cells can be specified in the constructor:

universe = openmc.Universe(cells=[cell1, cell2, cell3])

# This is equivalent
universe = openmc.Universe()
universe.add_cells([cell1, cell2])
universe.add_cell(cell3)

Universes are generally used in three ways:

  1. To be assigned to a Geometry object (see Exporting a Geometry Model),
  2. To be assigned as the fill for a cell via the Cell.fill attribute, and
  3. To be used in a regular arrangement of universes in a lattice.

Once a universe is constructed, it can actually be used to determine what cell or material is found at a given location by using the Universe.find() method, which returns a list of universes, cells, and lattices which are traversed to find a given point. The last element of that list would contain the lowest-level cell at that location:

>>> universe.find((0., 0., 0.))[-1]
Cell
        ID             =    10000
        Name           =    cell 1
        Fill           =    Material 10000
        Region         =    -10000
        Rotation       =    None
        Temperature    =    None
        Translation    =    None

As you are building a geometry, it is also possible to display a plot of single universe using the Universe.plot() method. This method requires that you have matplotlib installed.

6.4. Lattices

Many particle transport models involve repeated structures that occur in a regular pattern such as a rectangular or hexagonal lattice. In such a case, it would be cumbersome to have to define the boundaries of each of the cells to be filled with a universe. OpenMC provides a means to define lattice structures through the openmc.RectLattice and openmc.HexLattice classes.

6.4.1. Rectangular Lattices

A rectangular lattice defines a two-dimensional or three-dimensional array of universes that are filled into rectangular prisms (lattice elements) each of which has the same width, length, and height. To completely define a rectangular lattice, one needs to specify

  • The coordinates of the lower-left corner of the lattice (RectLattice.lower_left),
  • The pitch of the lattice, i.e., the distance between the center of adjacent lattice elements (RectLattice.pitch),
  • What universes should fill each lattice element (RectLattice.universes), and
  • A universe that is used to fill any lattice position outside the well-defined portion of the lattice (RectLattice.outer).

For example, to create a 3x3 lattice centered at the origin in which each lattice element is 5cm by 5cm and is filled by a universe u, one could run:

lattice = openmc.RectLattice()
lattice.lower_left = (-7.5, -7.5)
lattice.pitch = (5.0, 5.0)
lattice.universes = [[u, u, u],
                     [u, u, u],
                     [u, u, u]]

Note that because this is a two-dimensional lattice, the lower-left coordinates and pitch only need to specify the \(x,y\) values. The order that the universes appear is such that the first row corresponds to lattice elements with the highest \(y\) -value. Note that the RectLattice.universes attribute expects a doubly-nested iterable of type openmc.Universe — this can be normal Python lists, as shown above, or a NumPy array can be used as well:

lattice.universes = np.tile(u, (3, 3))

For a three-dimensional lattice, the \(x,y,z\) coordinates of the lower-left coordinate need to be given and the pitch should also give dimensions for all three axes. For example, to make a 3x3x3 lattice where the bottom layer is universe u, the middle layer is universe q and the top layer is universe z would look like:

lat3d = openmc.RectLattice()
lat3d.lower_left = (-7.5, -7.5, -7.5)
lat3d.pitch = (5.0, 5.0, 5.0)
lat3d.universes = [
    [[u, u, u],
     [u, u, u],
     [u, u, u]],
    [[q, q, q],
     [q, q, q],
     [q, q, q]],
    [[z, z, z],
     [z, z, z]
     [z, z, z]]]

Again, using NumPy can make things easier:

lat3d.universes = np.empty((3, 3, 3), dtype=openmc.Universe)
lat3d.universes[0, ...] = u
lat3d.universes[1, ...] = q
lat3d.universes[2, ...] = z

Finally, it’s possible to specify that lattice positions that aren’t normally without the bounds of the lattice be filled with an “outer” universe. This allows one to create a truly infinite lattice if desired. An outer universe is set with the RectLattice.outer attribute.

6.4.2. Hexagonal Lattices

OpenMC also allows creation of 2D and 3D hexagonal lattices. Creating a hexagonal lattice is similar to creating a rectangular lattice with a few differences:

  • The center of the lattice must be specified (HexLattice.center).
  • For a 2D hexagonal lattice, a single value for the pitch should be specified, although it still needs to appear in a list. For a 3D hexagonal lattice, the pitch in the radial and axial directions should be given.
  • For a hexagonal lattice, the HexLattice.universes attribute cannot be given as a NumPy array for reasons explained below.
  • As with rectangular lattices, the HexLattice.outer attribute will specify an outer universe.

For a 2D hexagonal lattice, the HexLattice.universes attribute should be set to a two-dimensional list of universes filling each lattice element. Each sub-list corresponds to one ring of universes and is ordered from the outermost ring to the innermost ring. The universes within each sub-list are ordered from the “top” (position with greatest y value) and proceed in a clockwise fashion around the ring. The HexLattice.show_indices() static method can be used to help figure out how to place universes:

>>> print(openmc.HexLattice.show_indices(3))
            (0, 0)
      (0,11)      (0, 1)
(0,10)      (1, 0)      (0, 2)
      (1, 5)      (1, 1)
(0, 9)      (2, 0)      (0, 3)
      (1, 4)      (1, 2)
(0, 8)      (1, 3)      (0, 4)
      (0, 7)      (0, 5)
            (0, 6)

Note that by default, hexagonal lattices are positioned such that each lattice element has two faces that are parallel to the \(y\) axis. As one example, to create a three-ring lattice centered at the origin with a pitch of 10 cm where all the lattice elements centered along the \(y\) axis are filled with universe u and the remainder are filled with universe q, the following code would work:

hexlat = openmc.HexLattice()
hexlat.center = (0, 0)
hexlat.pitch = [10]

outer_ring = [u, q, q, q, q, q, u, q, q, q, q, q]
middle_ring = [u, q, q, u, q, q]
inner_ring = [u]
hexlat.universes = [outer_ring, middle_ring, inner_ring]

If you need to create a hexagonal boundary (composed of six planar surfaces) for a hexagonal lattice, openmc.model.hexagonal_prism() can be used.

6.5. Exporting a Geometry Model

Once you have finished building your geometry by creating surfaces, cell, and, if needed, lattices, the last step is to create an instance of openmc.Geometry and export it to an XML file that the openmc executable can read using the Geometry.export_to_xml() method. This can be done as follows:

geom = openmc.Geometry(root_univ)
geom.export_to_xml()

# This is equivalent
geom = openmc.Geometry()
geom.root_universe = root_univ
geom.export_to_xml()

Note that it’s not strictly required to manually create a root universe. You can also pass a list of cells to the openmc.Geometry constructor and it will handle creating the unverse:

geom = openmc.Geometry([cell1, cell2, cell3])
geom.export_to_xml()

6.6. Using CAD-based Geometry

OpenMC relies on the Direct Accelerated Geometry Monte Carlo toolkit (DAGMC) to represent CAD-based geometry in a surface mesh format. A DAGMC run can be enabled in OpenMC by setting the dagmc property to True in the model Settings either via the Python openmc.settings Python class:

settings = openmc.Settings()
settings.dagmc = True

or in the settings.xml file:

<dagmc>true</dagmc>

With dagmc set to true, OpenMC will load the DAGMC model (from a local file named dagmc.h5m) when initializing a simulation. If a geometry.xml is present as well, it will be ignored.

Note: DAGMC geometries used in OpenMC are currently required to be clean, meaning that all surfaces have been imprinted and merged successfully and that the model is watertight. Future implementations of DAGMC geometry will support small volume overlaps and un-merged surfaces.