2. Geometry

2.1. Constructive Solid Geometry

OpenMC uses a technique known as constructive solid geometry (CSG) to build arbitrarily complex three-dimensional models in Euclidean space. In a CSG model, every unique object is described as the union and/or intersection of half-spaces created by bounding surfaces. Every surface divides all of space into exactly two half-spaces. We can mathematically define a surface as a collection of points that satisfy an equation of the form \(f(x,y,z) = 0\) where \(f(x,y,z)\) is a given function. All coordinates for which \(f(x,y,z) < 0\) are referred to as the negative half-space (or simply the negative side) and coordinates for which \(f(x,y,z) > 0\) are referred to as the positive half-space.

Let us take the example of a sphere centered at the point \((x_0,y_0,z_0)\) with radius \(R\). One would normally write the equation of the sphere as

(1)\[(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2\]

By subtracting the right-hand term from both sides of equation (1), we can then write the surface equation for the sphere:

(2)\[f(x,y,z) = (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\]

One can confirm that any point inside this sphere will correspond to \(f(x,y,z) < 0\) and any point outside the sphere will correspond to \(f(x,y,z) > 0\).

In OpenMC, every surface defined by the user is assigned an integer to uniquely identify it. We can then refer to either of the two half-spaces created by a surface by a combination of the unique ID of the surface and a positive/negative sign. Figure 5 shows an example of an ellipse with unique ID 1 dividing space into two half-spaces.

../_images/halfspace.svg

Figure 5: Example of an ellipse and its associated half-spaces.

References to half-spaces created by surfaces are used to define regions of space of uniform composition, which are then assigned to cells. OpenMC allows regions to be defined using union, intersection, and complement operators. As in MCNP, the intersection operator is implicit as doesn’t need to be written in a region specification. A defined region is then associated with a material composition in a cell. Figure 6 shows an example of a cell region defined as the intersection of an ellipse and two planes.

../_images/union.svg

Figure 6: The shaded region represents a cell bounded by three surfaces.

The ability to form regions based on bounding quadratic surfaces enables OpenMC to model arbitrarily complex three-dimensional objects. In practice, one is limited only by the different surface types available in OpenMC. The following table lists the available surface types, the identifier used to specify them in input files, the corresponding surface equation, and the input parameters needed to fully define the surface.

Surface types available in OpenMC.

Surface

Identifier

Equation

Parameters

Plane perpendicular to \(x\)-axis

x-plane

\(x - x_0 = 0\)

\(x_0\)

Plane perpendicular to \(y\)-axis

y-plane

\(y - y_0 = 0\)

\(y_0\)

Plane perpendicular to \(z\)-axis

z-plane

\(z - z_0 = 0\)

\(z_0\)

Arbitrary plane

plane

\(Ax + By + Cz = D\)

\(A\;B\;C\;D\)

Infinite cylinder parallel to \(x\)-axis

x-cylinder

\((y-y_0)^2 + (z-z_0)^2 = R^2\)

\(y_0\;z_0\;R\)

Infinite cylinder parallel to \(y\)-axis

y-cylinder

\((x-x_0)^2 + (z-z_0)^2 = R^2\)

\(x_0\;z_0\;R\)

Infinite cylinder parallel to \(z\)-axis

z-cylinder

\((x-x_0)^2 + (y-y_0)^2 = R^2\)

\(x_0\;y_0\;R\)

Sphere

sphere

\((x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 = R^2\)

\(x_0 \; y_0 \; z_0 \; R\)

Cone parallel to the \(x\)-axis

x-cone

\((y-y_0)^2 + (z-z_0)^2 = R^2(x-x_0)^2\)

\(x_0 \; y_0 \; z_0 \; R^2\)

Cone parallel to the \(y\)-axis

y-cone

\((x-x_0)^2 + (z-z_0)^2 = R^2(y-y_0)^2\)

\(x_0 \; y_0 \; z_0 \; R^2\)

Cone parallel to the \(z\)-axis

z-cone

\((x-x_0)^2 + (y-y_0)^2 = R^2(z-z_0)^2\)

\(x_0 \; y_0 \; z_0 \; R^2\)

General quadric surface

quadric

\(Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0\)

\(A \; B \; C \; D \; E \; F \; G \; H \; J \; K\)

Torus parallel to the \(x\)-axis

x-torus

\((x-x_0)^2/B^2+\frac{( \sqrt{(y-y_0)^2+(z-z_0)^2} - A)^2}{C^2} - 1 = 0\)

\(x_0 \; y_0 \; z_0 \; A \; B \; C\)

Torus parallel to the \(y\)-axis

y-torus

\((y-y_0)^2/B^2+\frac{( \sqrt{(x-x_0)^2+(z-z_0)^2} - A)^2}{C^2} - 1 = 0\)

\(x_0 \; y_0 \; z_0 \; A \; B \; C\)

Torus parallel to the \(z\)-axis

z-torus

\((z-z_0)^2/B^2+\frac{( \sqrt{(x-x_0)^2+(y-y_0)^2} - A)^2}{C^2} - 1 = 0\)

\(x_0 \; y_0 \; z_0 \; A \; B \; C\)

2.1.1. Universes

OpenMC supports universe-based geometry similar to the likes of MCNP and Serpent. This capability enables user to model any identical repeated structures once and then fill them in various spots in the geometry. A prototypical example of a repeated structure would be a fuel pin within a fuel assembly or a fuel assembly within a core.

Each cell in OpenMC can either be filled with a normal material or with a universe. If the cell is filled with a universe, only the region of the universe that is within the defined boundaries of the parent cell will be present in the geometry. That is to say, even though a collection of cells in a universe may extend to infinity, not all of the universe will be “visible” in the geometry since it will be truncated by the boundaries of the cell that contains it.

When a cell is filled with a universe, it is possible to specify that the universe filling the cell should be rotated and translated. This is done through the rotation and translation attributes on a cell (note though that these can only be specified on a cell that is filled with another universe, not a material).

It is not necessary to use or assign universes in a geometry if there are no repeated structures. Any cell in the geometry that is not assigned to a specified universe is automatically part of the base universe whose coordinates are just the normal coordinates in Euclidean space.

2.1.2. Lattices

Often times, repeated structures in a geometry occur in a regular pattern such as a rectangular or hexagonal lattice. In such a case, it would be cumbersome for a user to have to define the boundaries of each of the cells to be filled with a universe. Thus, OpenMC provides a lattice capability similar to that used in MCNP and Serpent.

The implementation of lattices is similar in principle to universes — instead of a cell being filled with a universe, the user can specify that it is filled with a finite lattice. The lattice is then defined by a two-dimensional array of universes that are to fill each position in the lattice. A good example of the use of lattices and universes can be seen in the OpenMC model for the Monte Carlo Performance benchmark.

2.2. Computing the Distance to Nearest Boundary

One of the most basic algorithms in any Monte Carlo code is determining the distance to the nearest surface within a cell. Since each cell is defined by the surfaces that bound it, if we compute the distance to all surfaces bounding a cell, we can determine the nearest one.

With the possibility of a particle having coordinates on multiple levels (universes) in a geometry, we must exercise care when calculating the distance to the nearest surface. Each different level of geometry has a set of boundaries with which the particle’s direction of travel may intersect. Thus, it is necessary to check the distance to the surfaces bounding the cell in each level. This should be done starting the highest (most global) level going down to the lowest (most local) level. That ensures that if two surfaces on different levels are coincident, by default the one on the higher level will be selected as the nearest surface. Although they are not explicitly defined, it is also necessary to check the distance to surfaces representing lattice boundaries if a lattice exists on a given level.

The following procedure is used to calculate the distance to each bounding surface. Suppose we have a particle at \((x_0,y_0,z_0)\) traveling in the direction \(u_0,v_0,w_0\). To find the distance \(d\) to a surface \(f(x,y,z) = 0\), we need to solve the equation:

(3)\[f(x_0 + du_0, y_0 + dv_0, z_0 + dw_0) = 0\]

If no solutions to equation (3) exist or the only solutions are complex, then the particle’s direction of travel will not intersect the surface. If the solution to equation (3) is negative, this means that the surface is “behind” the particle, i.e. if the particle continues traveling in its current direction, it will not hit the surface. The complete derivation for different types of surfaces used in OpenMC will be presented in the following sections.

Since \(f(x,y,z)\) in general is quadratic in \(x\), \(y\), and \(z\), this implies that \(f(x_0 + du_0, y + dv_0, z + dw_0)\) is quadratic in \(d\). Thus we expect at most two real solutions to (3). If no solutions to (3) exist or the only solutions are complex, then the particle’s direction of travel will not intersect the surface. If the solution to (3) is negative, this means that the surface is “behind” the particle, i.e. if the particle continues traveling in its current direction, it will not hit the surface.

Once a distance has been computed to a surface, we need to check if it is closer than previously-computed distances to surfaces. Unfortunately, we cannot just use the minimum function because some of the calculated distances, which should be the same in theory (e.g. coincident surfaces), may be slightly different due to the use of floating-point arithmetic. Consequently, we should first check for floating-point equality of the current distance calculated and the minimum found thus far. This is done by checking if

(4)\[\frac{| d - d_{min} |}{d_{min}} < \epsilon\]

where \(d\) is the distance to a surface just calculated, \(d_{min}\) is the minimum distance found thus far, and \(\epsilon\) is a small number. In OpenMC, this parameter is set to \(\epsilon = 10^{-14}\) since all floating calculations are done on 8-byte floating point numbers.

2.2.1. Plane Perpendicular to an Axis

The equation for a plane perpendicular to, for example, the x-axis is simply \(x - x_0 = 0\). As such, we need to solve \(x + du - x_0 = 0\). The solution for the distance is

(5)\[d = \frac{x_0 - x}{u}\]

Note that if the particle’s direction of flight is parallel to the x-axis, i.e. \(u = 0\), the distance to the surface will be infinity. While the example here was for a plane perpendicular to the x-axis, the same formula can be applied for the surfaces \(y = y_0\) and \(z = z_0\).

2.2.2. Generic Plane

The equation for a generic plane is \(Ax + By + Cz = D\). Thus, we need to solve the equation \(A(x + du) + B(y + dv) + C(z + dw) = D\). The solution to this equation for the distance is

(6)\[d = \frac{D - Ax - By - Cz}{Au + Bv + Cw}\]

Again, we need to check whether the denominator is zero. If so, this means that the particle’s direction of flight is parallel to the plane and it will therefore never hit the plane.

2.2.3. Cylinder Parallel to an Axis

The equation for a cylinder parallel to, for example, the x-axis is \((y - y_0)^2 + (z - z_0)^2 = R^2\). Thus, we need to solve \((y + dv - y_0)^2 + (z + dw - z_0)^2 = R^2\). Let us define \(\bar{y} = y - y_0\) and \(\bar{z} = z - z_0\). We then have

(7)\[(\bar{y} + dv)^2 + (\bar{z} + dw)^2 = R^2\]

Expanding equation (7) and rearranging terms, we obtain

(8)\[(v^2 + w^2) d^2 + 2 (\bar{y}v + \bar{z}w) d + (\bar{y}^2 + \bar{z}^2 - R^2) = 0\]

This is a quadratic equation for \(d\). To simplify notation, let us define \(a = v^2 + w^2\), \(k = \bar{y}v + \bar{z}w\), and \(c = \bar{y}^2 + \bar{z}^2 - R^2\). Thus, the distance is just the solution to \(ad^2 + 2kd + c = 0\):

(9)\[d = \frac{-k \pm \sqrt{k^2 - ac}}{a}\]

A few conditions must be checked for. If \(a = 0\), this means the particle is parallel to the cylinder and will thus never intersect it. Also, if \(k^2 - ac < 0\), this means that both solutions to the quadratic are complex. In physical terms, this means that the ray along which the particle is traveling does not make any intersections with the cylinder.

If we do have intersections and \(c < 0\), this means that the particle is inside the cylinder. Thus, one solution should be positive and one should be negative. Clearly, the positive distance will occur when the sign on the square root of the discriminant is positive since \(a > 0\).

If we have intersections and \(c > 0\) this means that the particle is outside the cylinder. Thus, the solutions to the quadratic are either both positive or both negative. If they are both positive, the smaller (closer) one will be the solution with a negative sign on the square root of the discriminant.

The same equations and logic here can be used for cylinders that are parallel to the y- or z-axis with appropriate substitution of constants.

2.2.4. Sphere

The equation for a sphere is \((x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2\). Thus, we need to solve the equation

(10)\[(x + du - x_0)^2 + (y + dv - y_0)^2 + (z + dw - z_0)^2 = R^2\]

Let us define \(\bar{x} = x - x_0\), \(\bar{y} = y - y_0\), and \(\bar{z} = z - z_0\). We then have

(11)\[(\bar{x} + du)^2 + (\bar{y} + dv)^2 + (\bar{z} - dw)^2 = R^2\]

Expanding equation (11) and rearranging terms, we obtain

(12)\[d^2 + 2 (\bar{x}u + \bar{y}v + \bar{z}w) d + (\bar{x}^2 + \bar{y}^2 + \bar{z}^2 - R^2) = 0\]

This is a quadratic equation for \(d\). To simplify notation, let us define \(k = \bar{x}u + \bar{y}v + \bar{z}w\) and \(c = \bar{x}^2 + \bar{y}^2 + \bar{z}^2 - R^2\). Thus, the distance is just the solution to \(d^2 + 2kd + c = 0\):

(13)\[d = -k \pm \sqrt{k^2 - c}\]

If the discriminant \(k^2 - c < 0\), this means that both solutions to the quadratic are complex. In physical terms, this means that the ray along which the particle is traveling does not make any intersections with the sphere.

If we do have intersections and \(c < 0\), this means that the particle is inside the sphere. Thus, one solution should be positive and one should be negative. The positive distance will occur when the sign on the square root of the discriminant is positive. If we have intersections but \(c > 0\) this means that the particle is outside the sphere. The solutions to the quadratic will then be either both positive or both negative. If they are both positive, the smaller (closer) one will be the solution with a negative sign on the square root of the discriminant.

2.2.5. Cone Parallel to an Axis

The equation for a cone parallel to, for example, the x-axis is \((y - y_0)^2 + (z - z_0)^2 = R^2(x - x_0)^2\). Thus, we need to solve \((y + dv - y_0)^2 + (z + dw - z_0)^2 = R^2(x + du - x_0)^2\). Let us define \(\bar{x} = x - x_0\), \(\bar{y} = y - y_0\), and \(\bar{z} = z - z_0\). We then have

(14)\[(\bar{y} + dv)^2 + (\bar{z} + dw)^2 = R^2(\bar{x} + du)^2\]

Expanding equation (14) and rearranging terms, we obtain

(15)\[(v^2 + w^2 - R^2u^2) d^2 + 2 (\bar{y}v + \bar{z}w - R^2\bar{x}u) d + (\bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2) = 0\]

Defining the terms

(16)\[a = v^2 + w^2 - R^2u^2 k = \bar{y}v + \bar{z}w - R^2\bar{x}u c = \bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2\]

we then have the simple quadratic equation \(ad^2 + 2kd + c = 0\) which can be solved as described in Cylinder Parallel to an Axis.

2.2.6. General Quadric

The equation for a general quadric surface is \(Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0\). Thus, we need to solve the equation

(17)\[A(x+du)^2 + B(y+dv)^2 + C(z+dw)^2 + D(x+du)(y+dv) + E(y+dv)(z+dw) + \\ F(x+du)(z+dw) + G(x+du) + H(y+dv) + J(z+dw) + K = 0\]

Expanding equation (17) and rearranging terms, we obtain

(18)\[d^2(uv + vw + uw) + 2d(Aux + Bvy + Cwx + (D(uv + vx) + E(vz + wy) + \\ F(wx + uz))/2) + (x(Ax + Dy) + y(By + Ez) + z(Cz + Fx)) = 0\]

Defining the terms

(19)\[a = uv + vw + uw k = Aux + Bvy + Cwx + (D(uv + vx) + E(vz + wy) + F(wx + uz))/2 c = x(Ax + Dy) + y(By + Ez) + z(Cz + Fx)\]

we then have the simple quadratic equation \(ad^2 + 2kd + c = 0\) which can be solved as described in Cylinder Parallel to an Axis.

2.2.7. Torus Parallel to an Axis

The equation for a torus parallel to, for example, the x-axis is

(20)\[\frac{(x-x_0)^2}{B^2} + \frac{(\sqrt{(y-y_0)^2 + (z-z_0)^2} - A)^2}{C^2} - 1 = 0.\]

First, it needs to be cast into a polynomial form. Rearranging terms,

(21)\[ (D\bar{x}^2 + \bar{y}^2 + \bar{z}^2 + A^2 - C^2)^2 = 4A^2(\bar{y}^2 + \bar{z}^2)\]

where \(D = (C/B)^2\), \(\bar{x} = x - x_0\), \(\bar{y} = y - y_0\), and \(\bar{z} = z - z_0\). To find the distance to the surface, we thus need to solve

(22)\[(D(\bar{x} + du)^2 + (\bar{y} + dv)^2 + (\bar{z} + dw)^2 + A^2 - C^2)^2 = 4A^2((\bar{y} + dv)^2 + (\bar{z} + dw)^2).\]

Expanding and collecting like powers of \(d\) yields

(23)\[(c_2d^2 + c_1d + c_0)^2 = c_2'd^2 + c_1'd + c_0'\]

where

(24)\[\begin{align*} c_2 &= Du^2 + v^2 + w^2 \\ c_1 &= 2(Du\bar{x} + v\bar{y} + w\bar{z}) \\ c_0 &= D\bar{x}^2 + \bar{y}^2 + \bar{z}^2 + A^2 - C^2 \\ c_2' &= 4A^2 (v^2 + w^2) \\ c_1' &= 8A^2 (v\bar{y} + w\bar{z}) \\ c_0' &= 4A^2(\bar{y}^2 + \bar{z}^2). \end{align*}\]

Expanding the left-hand side and collecting like powers of \(d\) on one side, we obtain

(25)\[(c_2^2)d^4 + (2c_1c_2)d^3 + (c_1^2 + 2c_0c_2 - c_2')d^2 + (2c_0c_1 - c_1')d + (c_0^2 - c_0') = 0.\]

The above equation is a fourth-order (quartic) polynomial equation. Although there is an analytical solution to the general quartic equation, it can be subject to roundoff errors when evaluated numerically. OpenMC uses an external quartic equation solver developed by Orellana and De Michele that is based on the decomposition of the quartic polynomial into two quadratics.

2.3. Finding a Cell Given a Point

Another basic algorithm is to determine which cell contains a given point in the global coordinate system, i.e. if the particle’s position is \((x,y,z)\), what cell is it currently in. This is done in the following manner in OpenMC. With the possibility of multiple levels of coordinates, we must perform a recursive search for the cell. First, we start in the highest (most global) universe, which we call the base universe, and loop over each cell within that universe. For each cell, we check whether the specified point is inside the cell using the algorithm described in Finding a Lattice Tile. If the cell is filled with a normal material, the search is done and we have identified the cell containing the point. If the cell is filled with another universe, we then search all cells within that universe to see if any of them contain the specified point. If the cell is filled with a lattice, the position within the lattice is determined, and then whatever universe fills that lattice position is recursively searched. The search ends once a cell containing a normal material is found that contains the specified point.

2.4. Finding a Lattice Tile

If a particle is inside a lattice, its position inside the lattice must be determined before assigning it to a cell. Throughout this section, the volumetric units of the lattice will be referred to as “tiles”. Tiles are identified by thier indices, and the process of discovering which tile contains the particle is referred to as “indexing”.

2.4.1. Rectilinear Lattice Indexing

Indices are assigned to tiles in a rectilinear lattice based on the tile’s position along the \(x\), \(y\), and \(z\) axes. Figure 7 maps the indices for a 2D lattice. The indices, (1, 1), map to the lower-left tile. (5, 1) and (5, 5) map to the lower-right and upper-right tiles, respectively.

../_images/rect_lat.svg

Figure 7: Rectilinear lattice tile indices.

In general, a lattice tile is specified by the three indices, \((i_x, i_y, i_z)\). If a particle’s current coordinates are \((x, y, z)\) then the indices can be determined from these formulas:

(26)\[i_x = \left \lceil \frac{x - x_0}{p_0} \right \rceil i_y = \left \lceil \frac{y - y_0}{p_1} \right \rceil i_z = \left \lceil \frac{z - z_0}{p_2} \right \rceil\]

where \((x_0, y_0, z_0)\) are the coordinates to the lower-left-bottom corner of the lattice, and \(p_0, p_1, p_2\) are the pitches along the \(x\), \(y\), and \(z\) axes, respectively.

2.4.2. Hexagonal Lattice Indexing

A skewed coordinate system is used for indexing hexagonal lattice tiles. Rather than a \(y\)-axis, another axis is used that is rotated 30 degrees counter-clockwise from the \(y\)-axis. This axis is referred to as the \(\alpha\)-axis. Figure 8 shows how 2D hexagonal tiles are mapped with the \((x, \alpha)\) basis. In this system, (0, 0) maps to the center tile, (0, 2) to the top tile, and (2, -1) to the middle tile on the right side.

../_images/hex_lat.svg

Figure 8: Hexagonal lattice tile indices.

Unfortunately, the indices cannot be determined with one simple formula as before. Indexing requires a two-step process, a coarse step which determines a set of four tiles that contains the particle and a fine step that determines which of those four tiles actually contains the particle.

In the first step, indices are found using these formulas:

(27)\[\alpha = -\frac{x}{\sqrt{3}} + y i_x^* = \left \lfloor \frac{x}{p_0 \sqrt{3} / 2} \right \rfloor i_\alpha^* = \left \lfloor \frac{\alpha}{p_0} \right \rfloor\]

where \(p_0\) is the lattice pitch (in the \(x\)-\(y\) plane). The true index of the particle could be \((i_x^*, i_\alpha^*)\), \((i_x^* + 1, i_\alpha^*)\), \((i_x^*, i_\alpha^* + 1)\), or \((i_x^* + 1, i_\alpha^* + 1)\).

The second step selects the correct tile from that neighborhood of 4. OpenMC does this by calculating the distance between the particle and the centers of each of the 4 tiles, and then picking the closest tile. This works because regular hexagonal tiles form a Voronoi tessellation which means that all of the points within a tile are closest to the center of that same tile.

Indexing along the \(z\)-axis uses the same method from rectilinear lattices, i.e.

(28)\[i_z = \left \lceil \frac{z - z_0}{p_2} \right \rceil\]

2.5. Determining if a Coordinate is in a Cell

To determine which cell a particle is in given its coordinates, we need to be able to check whether a given cell contains a point. The algorithm for determining if a cell contains a point is as follows. For each surface that bounds a cell, we determine the particle’s sense with respect to the surface. As explained earlier, if we have a point \((x_0,y_0,z_0)\) and a surface \(f(x,y,z) = 0\), the point is said to have negative sense if \(f(x_0,y_0,z_0) < 0\) and positive sense if \(f(x_0,y_0,z_0) > 0\). If for all surfaces, the sense of the particle with respect to the surface matches the specified sense that defines the half-space within the cell, then the point is inside the cell. Note that this algorithm works only for simple cells defined as intersections of half-spaces.

It may help to illustrate this algorithm using a simple example. Let’s say we have a cell defined as

<surface id="1" type="sphere"  coeffs="0 0 0 10" />
<surface id="2" type="x-plane" coeffs="-3" />
<surface id="3" type="y-plane" coeffs="2" />
<cell id="1" surfaces="-1 2 -3" />

This means that the cell is defined as the intersection of the negative half space of a sphere, the positive half-space of an x-plane, and the negative half-space of a y-plane. Said another way, any point inside this cell must satisfy the following equations

(29)\[x^2 + y^2 + z^2 - 10^2 < 0 \\ x - (-3) > 0 \\ y - 2 < 0\]

In order to determine if a point is inside the cell, we would substitute its coordinates into equation (29). If the inequalities are satisfied, than the point is indeed inside the cell.

2.6. Handling Surface Crossings

A particle will cross a surface if the distance to the nearest surface is closer than the distance sampled to the next collision. A number of things happen when a particle hits a surface. First, we need to check if a non-transmissive boundary condition has been applied to the surface. If a vacuum boundary condition has been applied, the particle is killed and any surface current tallies are scored to as needed. If a reflective boundary condition has been applied to the surface, surface current tallies are scored to and then the particle’s direction is changed according to the procedure in Reflective Boundary Conditions. Note that the white boundary condition can be considered as the special case of reflective boundary condition, where the same processing method will be applied to deal with the surface current tallies scoring, except for determining the changes of particle’s direction according to the procedures in White Boundary Conditions.

Next, we need to determine what cell is beyond the surface in the direction of travel of the particle so that we can evaluate cross sections based on its material properties. At initialization, a list of neighboring cells is created for each surface in the problem as described in Building Neighbor Lists. The algorithm outlined in Finding a Cell Given a Point is used to find a cell containing the particle with one minor modification; rather than searching all cells in the base universe, only the list of neighboring cells is searched. If this search is unsuccessful, then a search is done over every cell in the base universe.

2.7. Building Neighbor Lists

Neighbor lists are data structures that are used to accelerate geometry searches when a particle crosses a boundary. Namely, they are used to constrain the number of cells that must be searched in order to determine which cell a particle is crossing into. Earlier versions of OpenMC relied on “surface-based” neighbor lists, where the cells that are adjacent to each surface are stored in lists, one for each side of a surface. As of version 0.11, OpenMC switched to using “cell-based” neighbor lists. For each cell, a list of the adjacent cells is stored and then used to limit future searches. Unlike surface-based neighbor lists, cell-based neighbor lists cannot be computed prior to transport. Thus, cell-based neighbor lists in OpenMC grow dynamically as particles are transported through the geometry and cross surfaces. Special care must be taken to ensure that these dynamic neighbor lists are populated in a threadsafe manner. Full details of the implementation in OpenMC can be found in a paper by Harper et al.

2.8. Reflective Boundary Conditions

If the velocity of a particle is \(\mathbf{v}\) and it crosses a surface of the form \(f(x,y,z) = 0\) with a reflective boundary condition, it can be shown based on geometric arguments that the velocity vector will then become

(30)\[\mathbf{v'} = \mathbf{v} - 2 (\mathbf{v} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}\]

where \(\hat{\mathbf{n}}\) is a unit vector normal to the surface at the point of the surface crossing. The rationale for this can be understood by noting that \((\mathbf{v} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}\) is the projection of the velocity vector onto the normal vector. By subtracting two times this projection, the velocity is reflected with respect to the surface normal. Since the magnitude of the velocity of the particle will not change as it undergoes reflection, we can work with the direction of the particle instead, simplifying equation (30) to

(31)\[\mathbf{\Omega'} = \mathbf{\Omega} - 2 (\mathbf{\Omega} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}\]

where \(\mathbf{v} = || \mathbf{v} || \mathbf{\Omega}\). The direction of the surface normal will be the gradient of the surface at the point of crossing, i.e. \(\mathbf{n} = \nabla f(x,y,z)\). Substituting this into equation (31), we get

(32)\[\mathbf{\Omega'} = \mathbf{\Omega} - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} \nabla f\]

If we write the initial and final directions in terms of their vector components, \(\mathbf{\Omega} = (u,v,w)\) and \(\mathbf{\Omega'} = (u', v', w')\), this allows us to represent equation (31) as a series of equations:

(33)\[u' = u - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} \frac{\partial f}{\partial x} \\ v' = v - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} \frac{\partial f}{\partial y} \\ w' = w - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} \frac{\partial f}{\partial z}\]

One can then use equation (33) to develop equations for transforming a particle’s direction given the equation of the surface.

2.8.1. Plane Perpendicular to an Axis

For a plane that is perpendicular to an axis, the rule for reflection is almost so simple that no derivation is needed at all. Nevertheless, we will proceed with the derivation to confirm that the rules of geometry agree with our intuition. The gradient of the surface \(f(x,y,z) = x - x_0 = 0\) is simply \(\nabla f = (1, 0, 0)\). Note that this vector is already normalized, i.e. \(|| \nabla f || = 1\). The second two equations in (33) tell us that \(v\) and \(w\) do not change and the first tell us that

(34)\[u' = u - 2u = -u\]

We see that reflection for a plane perpendicular to an axis only entails negating the directional cosine for that axis.

2.8.2. Generic Plane

A generic plane has the form \(f(x,y,z) = Ax + By + Cz - D = 0\). Thus, the gradient to the surface is simply \(\nabla f = (A,B,C)\) whose norm squared is \(A^2 + B^2 + C^2\). This implies that

(35)\[\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = \frac{2(Au + Bv + Cw)}{A^2 + B^2 + C^2}\]

Substituting equation (35) into equation (33) gives us the form of the solution. For example, the x-component of the reflected direction will be

(36)\[u' = u - \frac{2A(Au + Bv + Cw)}{A^2 + B^2 + C^2}\]

2.8.3. Cylinder Parallel to an Axis

A cylinder parallel to, for example, the x-axis has the form \(f(x,y,z) = (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\). Thus, the gradient to the surface is

(37)\[\nabla f = 2 \left ( \begin{array}{c} 0 \\ y - y_0 \\ z - z_0 \end{array} \right ) = 2 \left ( \begin{array}{c} 0 \\ \bar{y} \\ \bar{z} \end{array} \right )\]

where we have introduced the constants \(\bar{y}\) and \(\bar{z}\). Taking the square of the norm of the gradient, we find that

(38)\[|| \nabla f ||^2 = 4 \bar{y}^2 + 4 \bar{z}^2 = 4 R^2\]

This implies that

(39)\[\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = \frac{\bar{y}v + \bar{z}w}{R^2}\]

Substituting equations (39) and (37) into equation (33) gives us the form of the solution. In this case, the x-component will not change. The y- and z-components of the reflected direction will be

(40)\[v' = v - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{y}}{R^2} \\ w' = w - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{z}}{R^2}\]

2.8.4. Sphere

The surface equation for a sphere has the form \(f(x,y,z) = (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\). Thus, the gradient to the surface is

(41)\[\nabla f = 2 \left ( \begin{array}{c} x - x_0 \\ y - y_0 \\ z - z_0 \end{array} \right ) = 2 \left ( \begin{array}{c} \bar{x} \\ \bar{y} \\ \bar{z} \end{array} \right )\]

where we have introduced the constants \(\bar{x}, \bar{y}, \bar{z}\). Taking the square of the norm of the gradient, we find that

(42)\[|| \nabla f ||^2 = 4 \bar{x}^2 + 4 \bar{y}^2 + 4 \bar{z}^2 = 4 R^2\]

This implies that

(43)\[\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = \frac{\bar{x}u + \bar{y}v + \bar{z}w}{R^2}\]

Substituting equations (43) and (41) into equation (33) gives us the form of the solution:

(44)\[u' = u - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{x} }{R^2} \\ v' = v - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{y} }{R^2} \\ w' = w - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{z} }{R^2}\]

2.8.5. Cone Parallel to an Axis

A cone parallel to, for example, the z-axis has the form \(f(x,y,z) = (x - x_0)^2 + (y - y_0)^2 - R^2(z - z_0)^2 = 0\). Thus, the gradient to the surface is

(45)\[\nabla f = 2 \left ( \begin{array}{c} x - x_0 \\ y - y_0 \\ -R^2(z - z_0) \end{array} \right ) = 2 \left ( \begin{array}{c} \bar{x} \\ \bar{y} \\ -R^2\bar{z} \end{array} \right )\]

where we have introduced the constants \(\bar{x}\), \(\bar{y}\), and \(\bar{z}\). Taking the square of the norm of the gradient, we find that

(46)\[|| \nabla f ||^2 = 4 \bar{x}^2 + \bar{y}^2 + 4 R^4 \bar{z}^2 \\ = 4 R^2 \bar{z}^2 + 4 R^4 \bar{z}^2 \\ = 4 R^2 (1 + R^2) \bar{z}^2\]

This implies that

(47)\[\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = \frac{\bar{x}u + \bar{y}v - R^2\bar{z}w}{R^2 (1 + R^2) \bar{z}^2}\]

Substituting equations (47) and (45) into equation (33) gives us the form of the solution:

(48)\[u' = u - \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w) \bar{x}}{R^2 (1 + R^2) \bar{z}^2} v' = v - \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w) \bar{y}}{R^2 (1 + R^2) \bar{z}^2} w' = w + \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w)}{R^2 (1 + R^2) \bar{z}}\]

2.8.6. General Quadric

A general quadric surface has the form \(f(x,y,z) = Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0\). Thus, the gradient to the surface is

(49)\[\nabla f = \left ( \begin{array}{c} 2Ax + Dy + Fz + G \\ 2By + Dx + Ez + H \\ 2Cz + Ey + Fx + J \end{array} \right ).\]

2.8.7. Torus Parallel to an Axis

A torus parallel to, for example, the x-axis has the form

(50)\[f(x,y,z) = \frac{(x-x_0)^2}{B^2} + \frac{(\sqrt{(y-y_0)^2 + (z-z_0)^2} - A)^2}{C^2} - 1.\]

The gradient to the surface is therefore

(51)\[\nabla f = \left ( \begin{array}{c} 2\bar{x}/B^2 \\ 2\bar{y}(g - A)/(C^2g) \\ 2\bar{z}(g - A)/(C^2g) \end{array} \right )\]

where \(g = \sqrt{\bar{y}^2 + \bar{z}^2}\) and, as always, \(\bar{x} = x - x_0\), \(\bar{y} = y - y_0\), and \(\bar{z} = z - z_0\).

2.9. White Boundary Conditions

The white boundary condition is usually applied in deterministic codes, where the particle will hit the surface and travel back with isotropic angular distribution. The change in particle’s direction is sampled from a cosine distribution instead of uniform. Figure 9 shows an example of cosine-distribution reflection on the arbitrary surface relative to the surface normal.

../_images/cosine-dist.png

Figure 9: Cosine-distribution reflection on an arbitrary surface.

The probability density function (pdf) for the reflected direction can be expressed as follows,

(52)\[f(\mu, \phi) d\mu d\phi = \frac{\mu}{\pi} d\mu d\phi = 2\mu d\mu \frac{d\phi}{2\pi}\]

where \(\mu = \cos \theta\) is the cosine of the polar angle between reflected direction and the normal to the surface; and \(\theta\) is the azimuthal angle in \([0,2\pi]\). We can separate the multivariate probability density into two separate univariate density functions, one for the cosine of the polar angle,

(53)\[f(\mu) = 2\mu\]

and one for the azimuthal angle,

(54)\[f(\phi) = \frac{1}{2\pi}.\]

Each of these density functions can be sampled by analytical inversion of the cumulative distribution distribution, resulting in the following sampling scheme:

(55)\[\mu = \sqrt{\xi_1} \\ \phi = 2\pi\xi_2\]

where \(\xi_1\) and \(\xi_2\) are uniform random numbers on \([0,1)\). With the sampled values of \(\mu\) and \(\phi\), the final reflected direction vector can be computed via rotation of the surface normal using the equations from Transforming a Particle’s Coordinates. The white boundary condition can be applied to any kind of surface, as long as the normal to the surface is known as in Reflective Boundary Conditions.