diff --git a/docs/source/_images/jt-truncated-cone.svg b/docs/source/_images/jt-truncated-cone.svg new file mode 100644 index 00000000000..9bd5b756192 --- /dev/null +++ b/docs/source/_images/jt-truncated-cone.svg @@ -0,0 +1,247 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rmajor + rminor + z0 + rlim + rcone + rs + z + zs + + diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 05cda4b6423..a03dda50ba5 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -178,6 +178,153 @@ 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`_. +Random Sphere Packing +--------------------- + +Some models will have regions consisting of a randomly arranged packing of +monosized spheres. This is particularly common in models simulating TRISO fuels. +OpenMC has two methods of random sphere packing: Random Sequential Packing (RSP) +and Close Random Packing (CRP). + +In RSP, sphere center coordinates are randomly generated one at a time +within the container's boundaries. New spheres which overlap with previously +generated spheres are rejected. This random coordinate selection continues +until the target number of random sphere centers is generated. + +In OpenMC, the CRP method implemented uses the Jodrey-Tory algorithm [Jodrey-Tory]_. +In Jodrey-Tory, each sphere is assigned two diameters: the inner diameter, +:math:`d_{inner}`, and the outer diameter, :math:`d_{outer}`, +which approach each other during the simulation. The inner diameter is defined +as the minimum center-to-center distance, or rod, between all spheres in the +container. The outer diameter is an arbitrary value that is larger than the +target sphere diameter. It is usually defined such that packing the container +with a target number of spheres of diameter :math:`d_{outer}` would yield a +packing fraction, :math:`pf`, of 1.0. + +The Jodrey-Tory algorithm prioritizes fixing the worst overlaps in each iteration. +Spheres are considered overlapping if the rod between them is less than :math:`d_{outer}`. +However, a sphere is only moved once each iteration, i.e., if sphere A is +overlapping with spheres B and C, but rod AC is shorter than rod AB, then +the overlap between A and C is addressed, but not A and B. Overlaps are eliminated +moving the spheres apart along the line joining their centers until they are +:math:`d_{outer}` apart. At the end of an iteration, :math:`d_{outer}` is +reduced according to: + +.. math:: + :label: jt-dout + + d_{out}^{i+1} = d_{out}^{i} - \frac{1}{2}^{j} * d_{out0} * \frac{k}{n} + +where k is the contraction rate, n is the number of spheres, and + +.. math:: + :label: jt-dout-j + + j = floor(-log_{10}(pf_{out} - pf_{in})) + +Iterations continue until the :math:`d_{outer}` and :math:`d_{inner}` converge +or the target :math:`pf` is reached, which is determined by the current value +of :math:`d_{inner}`. + +The OpenMC implementations of both RSP and CRP use a mesh over the domain +to speed up the nearest neighbor search by only searching for a sphere's +neighbors within that mesh cell. + +RSP performs better than CRP for lower :math:`pf`, but it becomes +prohibitively slow as it approaches its packing limit (~0.38). CRP can +achieve higher :math:`pf` of up to ~0.64 and scales better with increasing +:math:`pf`. If the desired :math:`pf` is below the threshold for which RSP +will be faster than CRP, only RSP is used. If a higher :math:`pf` is required, +RSP is used to pack the target number of spheres using a smaller than desired +final radius. This initial configuration of spheres is then used as a starting +point for CRP using Jodrey-Tory. + +After each sphere is moved, its position must be checked to confirm that it is +still within the boundaries of the container. In OpenMC, this check is performed +by calculating the limits of the container, which is the minimum and maximum +values the sphere center coordinate can have in each direction. In general, +limits can be defined as a surface that is parallel to the container's surface, +and is located a perpendicular distance inside the container equal to the +radius of the packed spheres. + +For example, if the container is a cube centered at (0, 0, 0) and defined with +six planes at :math:`x0, -x0, y0, -y0, z0`, and :math:`-z0`, then the limits are defined by a +cube centered on the origin with sides at :math:`(x0-r_{s}), -(x0-r_{s}), +y0-r_{s}, -(y0-r_{s}), z0-r_{s}`, and :math:`-(z0-r_{s})`, where :math:`r_{s}` +is the radius of the packed spheres. + +For most shapes this is straightforward, as the value of the limit is constant +for a given direction. This is not true for the limit in r of a conical container, +however. OpenMC allows for the packing of a truncated conical container in which +the top and bottom surfaces are capped by planes (similar to a cylinder). In +the case of such a container, the limit in r is dependent upon the sphere's +current z-coordinate relative to the top and bottom surfaces of the container. + +.. _fig-cone-container: + +.. figure:: ../_images/jt-truncated-cone.svg + :align: center + :figclass: align-center + :width: 400px + +Consider a conical container with dimensions labeled with the above convention. +In this convention, :math:`r_{major}` refers to the base at the +postive-most end of the cone, and :math:`r_{minor}` refers to the base at the +negative-most end of the cone. This does not refer to their magnitudes relative +to each other. + +The limit in r can be defined as a linear function of z using the equation for +the perpendicular distance between two parallel lines. With two parallel lines +defined by: + +.. math:: + :label: parallel-sys-eqn + + ar_{cone} + bz_{s} + c_{cone} = 0 + + ar_{lim} + bz_{s} + c_{lim} = 0 + +The parallel distance between them, :math:`d = r_{s}` is: + +.. math:: + :label: perpendicular-distance + + r_{s} = \frac{| c_{lim} - c_{cone} |}{\sqrt{a^2 + b^2}} + +The equation for :math:`r_{cone}` can be found using :math:`(r_minor, z0-z)` and +:math:`(r_major, z0+z)` as two points along the line, giving: + +.. math:: + :label: rcone-line + + z_{s} - (z_{0} - z) = (r_{cone} - r_{minor})\frac{2z}{r_{major} - r_{minor}} + + \frac{2z}{r_{major} - r_{minor}}r_{cone} - z_{s} - \frac{2z}{r_{major} - r_{minor}}r_{minor} + (z_{0} - z) = 0 + + mr_{cone} - z_{s} - mr_{minor} + (z_{0} - z) = 0 + +With :math:`m = \frac{2z}{r_{major} - r_{minor}}`, the constant :math:`c_{lim}` +is found to be: + +.. math:: + :label: rlim-clim + + r_{s} = \frac{|c_{lim} - c_{cone}|}{\sqrt{m^2 + 1^2}} + + c_{lim} = r_{s}\sqrt{m^2 + 1} \pm ((z_{0} - z) - r_{minor}) + +Which can be used to define the radial limit as a function of the sphere's current +z-position relative to the bases of the cone. The :math:`\pm` symbol is resolved +by choosing the sign that results in adding :math:`r_{minor}`, which gives an +:math:`r_{lim}` line the appropriate distance inside the conical container. + +.. math:: + :label: rlim-eqn + + r_{lim}(z_{s}) = \frac{1}{m}z_{s} - \frac{r_{s}}{m}\sqrt{m^2 + 1} + r_{minor} - \frac{(z_{0} - z)}{m} + + + ------------------------------------------ Computing the Distance to Nearest Boundary ------------------------------------------ @@ -1068,3 +1215,5 @@ surface is known as in :ref:`reflection`. .. _MCNP: https://mcnp.lanl.gov .. _Serpent: https://serpent.vtt.fi .. _Monte Carlo Performance benchmark: https://github.com/mit-crpg/benchmarks/tree/master/mc-performance/openmc +.. [Jodrey-Tory] W. S. Jodrey and E. M. Tory, "Computer simulation of close random + packing of equal spheres", Phys. Rev. A 32 (1985) 2347-2351. diff --git a/openmc/model/triso.py b/openmc/model/triso.py index f344b8edd42..0fc12ac4118 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -633,6 +633,270 @@ def repel_spheres(self, p, q, d, d_new): q[j] = (q[j] - c[j])*ul[1]/r + c[j] q[k] = np.clip(q[k], ll[0], ul[0]) +class _Cone(_Container): + """Conical container in which to pack spheres. + + Parameters + ---------- + length : float + Length of the conical container. + radius_major_ : float + Radius of the conical container at the positive end of the axis along + which the cone is aligned. + radius_minor_ : float + Radius of the conical container at the negative end of the axis along + which the cone is aligned. + axis : string + Axis along which the length of the cone is aligned. + sphere_radius : float + Radius of spheres to be packed in container. + center : Iterable of float + Cartesian coordinates of the center of the container. Default is + (0., 0., 0.) + + Attributes + ---------- + length : float + Length of the conical container. + radius_major_ : float + Radius of the conical container at the positive end of the axis along + which the cone is aligned. + radius_minor_ : float + Radius of the conical container at the negative end of the axis along + which the cone is aligned. + axis : string + Axis along which the length of the cylinder is aligned. + sphere_radius : float + Radius of spheres to be packed in container. + center : list of float + Cartesian coordinates of the center of the container. Default is + (0., 0., 0.) + shift : list of int + Rolled indices of the x-, y-, and z- coordinates of a sphere so the + configuration is aligned with the correct axis. No shift corresponds to + a cone along the z-axis. + cell_length : list of float + Length in x-, y-, and z- directions of each cell in mesh overlaid on + domain. + limits : list of float + Minimum and maximum distance in the direction parallel to the axis + where sphere center can be placed and maximum radial distance as a + function of position relative to the center of the conical container. + volume : float + Volume of the container. + + """ + + def __init__(self, length, radius_major, radius_minor, axis, + sphere_radius, center=(0., 0., 0.)): + super().__init__(sphere_radius, center) + self._shift = None + self.length = length + self.radius_major = radius_major + self.radius_minor = radius_minor + self.axis = axis + + @property + def length(self): + return self._length + + @length.setter + def length(self, length): + self._length = float(length) + self._limits = None + self._cell_length = None + + @property + def radius_major(self): + return self._radius_major + + @radius_major.setter + def radius_major(self, radius_major): + self._radius_major = float(radius_major) + self._limits = None + self._cell_length = None + + @property + def radius_minor(self): + return self._radius_minor + + @radius_minor.setter + def radius_minor(self, radius_minor): + self._radius_minor = float(radius_minor) + self._limits = None + self._cell_length = None + + @property + def axis(self): + return self._axis + + @axis.setter + def axis(self, axis): + self._axis = axis + self._shift = None + + @property + def shift(self): + if self._shift is None: + if self.axis == 'x': + self._shift = [1, 2, 0] + elif self.axis == 'y': + self._shift = [2, 0, 1] + else: + self._shift = [0, 1, 2] + return self._shift + + @property + def limits(self): + if self._limits is None: + z0 = self.center[self.shift[2]] + z = self.length/2 + r = self.sphere_radius + r_major = self.radius_major + r_minor = self.radius_minor + m = (2*z)/(r_major-r_minor) + # r limit depends on z-coord of packed sphere. + # the equation is broken into [factor, constant] + # so it can be stored without symbols + self._limits = [[z0 - z + r], + [z0 + z - r, + [1/m, + -(r/m)*(m**2 + 1)**0.5 - (z0-z)/m + r_minor]]] + return self._limits + + @limits.setter + def limits(self, limits): + self._limits = limits + + @property + def cell_length(self): + if self._cell_length is None: + h = 4*self.sphere_radius + #r_avg = 0.5*(self.radius_major + self.radius_minor) + i, j, k = self.shift + self._cell_length = [None]*3 + self._cell_length[i] = 2*self.radius_major/int(2*self.radius_major/h) + self._cell_length[j] = 2*self.radius_major/int(2*self.radius_major/h) + self._cell_length[k] = self.length/int(self.length/h) + return self._cell_length + + @property + def volume(self): + return (self.length*pi/3)*(self.radius_major**2 + + self.radius_major*self.radius_minor + + self.radius_minor**2) + + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) + + # Assume the simplest case where the conical volume is the + # intersection of the half-spaces of a cone and two planes + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 3: + raise ValueError + + # Identify the axis that the cone lies along + axis = region[0].surface.type[0] + + # Make sure the region is composed of a cone and two planes on the + # same axis + count = Counter(node.surface.type for node in region) + if count[axis + '-cone'] != 1 or count[axis + '-plane'] != 2: + raise ValueError + + # Sort the half-spaces by surface type + cone, p1, p2 = sorted(region, key=lambda x: x.surface.type) + + # Calculate the parameters for a cone along the x-axis + if axis == 'x': + if p1.surface.x0 > p2.surface.x0: + p1, p2 = p2, p1 + length = p2.surface.x0 - p1.surface.x0 + center = (p1.surface.x0 + length/2, cone.surface.y0, cone.surface.z0) + radius_major = (cone.surface.r2* + (p2.surface.x0 - cone.surface.x0)**2)**0.5 + radius_minor = (cone.surface.r2* + (p1.surface.x0 - cone.surface.x0)**2)**0.5 + + # Calculate the parameters for a cone along the y-axis + elif axis == 'y': + if p1.surface.y0 > p2.surface.y0: + p1, p2 = p2, p1 + length = p2.surface.y0 - p1.surface.y0 + center = (cone.surface.x0, p1.surface.y0 + length/2, cone.surface.z0) + radius_major = (cone.surface.r2* + (p2.surface.y0 - cone.surface.y0)**2)**0.5 + radius_minor = (cone.surface.r2* + (p1.surface.y0 - cone.surface.y0)**2)**0.5 + + # Calculate the parameters for a cone along the z-axis + else: + if p1.surface.z0 > p2.surface.z0: + p1, p2 = p2, p1 + length = p2.surface.z0 - p1.surface.z0 + center = (cone.surface.x0, cone.surface.y0, p1.surface.z0 + length/2) + radius_major = (cone.surface.r2* + (p2.surface.z0 - cone.surface.z0)**2)**0.5 + radius_minor = (cone.surface.r2* + (p1.surface.z0 - cone.surface.z0)**2)**0.5 + + # Make sure the half-spaces are on the correct side of the surfaces + if cone.side != '-' or p1.side != '+' or p2.side != '-': + raise ValueError + + # The region is the volume of a truncated cone, so create container + return _Cone(length, radius_major, radius_minor, axis, + sphere_radius, center) + + def random_point(self): + ll, ul = self.limits + i, j, k = self.shift + p = [None]*3 + p[k] = uniform(ll[0], ul[0]) + rl = p[k]*ul[1][0] + ul[1][1] + r = sqrt(uniform(0, rl**2)) + t = uniform(0, 2*pi) + p[i] = r*cos(t) + self.center[i] + p[j] = r*sin(t) + self.center[j] + return p + + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is + # equal to the outer diameter + s = (d_new - d)/2 + + v = (p - q)/d + p += s*v + q -= s*v + + # Enforce the rigid boundary by moving each sphere along the plane + # surface normal if they overlap, then checking and moving the + # sphere center to the new r limit if it overlaps with the side + # of the cone. + ll, ul = self.limits + c = self.center + i, j, k = self.shift + + p[k] = np.clip(p[k], ll[0], ul[0]) + r = sqrt((p[i] - c[i])**2 + (p[j] - c[j])**2) + rl = p[k]*ul[1][0] + ul[1][1] + if r > rl: + p[i] = (p[i] - c[i])*rl/r + c[i] + p[j] = (p[j] - c[j])*rl/r + c[j] + + q[k] = np.clip(q[k], ll[0], ul[0]) + r = sqrt((q[i] - c[i])**2 + (q[j] - c[j])**2) + rl = q[k]*ul[1][0] + ul[1][1] + if r > rl: + q[i] = (q[i] - c[i])*rl/r + c[i] + q[j] = (q[j] - c[j])*rl/r + c[j] class _SphericalShell(_Container): """Spherical shell container in which to pack spheres. @@ -1294,7 +1558,7 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, if not domain: raise ValueError('Could not map region {} to a container: supported ' 'container shapes are rectangular prism, cylinder, ' - 'sphere, and spherical shell.'.format(region)) + 'cone, sphere, and spherical shell.'.format(region)) # Determine the packing fraction/number of spheres volume = _volume_sphere(radius) @@ -1329,8 +1593,14 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, # Recalculate the limits for the initial random sequential packing using # the desired final sphere radius to ensure spheres are fully contained # within the domain during the close random pack - domain.limits = [[x - initial_radius + radius for x in domain.limits[0]], - [x + initial_radius - radius for x in domain.limits[1]]] + try: + domain.limits = [[x - initial_radius + radius for x in domain.limits[0]], + [x + initial_radius - radius for x in domain.limits[1]]] + except TypeError: + ll, ul = domain.limits + domain.limits = [[ll[0] - initial_radius + radius], + [ul[0] + initial_radius - radius, + [ul[1][0], ul[1][1] + initial_radius - radius]]] # Generate non-overlapping spheres for an initial inner radius using # random sequential packing algorithm diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index c1f8c039e57..4a16e264bab 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -17,6 +17,9 @@ {'shape': 'x_cylinder', 'volume': 1*pi*1**2}, {'shape': 'y_cylinder', 'volume': 1*pi*1**2}, {'shape': 'z_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'x_cone', 'volume': 2/3*pi*(1**2 + 1*3 + 3**2)}, + {'shape': 'y_cone', 'volume': 2/3*pi*(1**2 + 1*3 + 3**2)}, + {'shape': 'z_cone', 'volume': 2/3*pi*(1**2 + 1*3 + 3**2)}, {'shape': 'sphere', 'volume': 4/3*pi*1**3}, {'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)} ] @@ -74,6 +77,35 @@ def centers_z_cylinder(): return openmc.model.pack_spheres(radius=_RADIUS, region=region, pf=_PACKING_FRACTION, initial_pf=0.2) +@pytest.fixture(scope='module') +def centers_x_cone(): + cone = openmc.XCone(x0=-2, y0=1, z0=2, r2=1) + min_x = openmc.XPlane(-1) + max_x = openmc.XPlane(1) + region = +min_x & -max_x & -cone + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_y_cone(): + cone = openmc.YCone(x0=1, y0=-2, z0=2, r2=1) + min_y = openmc.YPlane(-1) + max_y = openmc.YPlane(1) + region = +min_y & -max_y & -cone + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_z_cone(): + cone = openmc.ZCone(x0=1, y0=2, z0=-2, r2=1) + min_z = openmc.ZPlane(-1) + max_z = openmc.ZPlane(1) + region = +min_z & -max_z & -cone + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + @pytest.fixture(scope='module') def centers_sphere(): @@ -153,6 +185,41 @@ def test_contained_z_cylinder(centers_z_cylinder): assert z_max < 1 or z_max == pytest.approx(1) assert z_min > 0 or z_min == pytest.approx(0) +def test_contained_x_cone(centers_x_cone): + """Make sure all spheres are entirely contained within the domain.""" + for p in centers_x_cone: + r = ((p[1] - 1)**2 + (p[2] - 2)**2)**0.5 + _RADIUS + rmax = p[0] + 1 - (0 - 1) + assert r < rmax or r == pytest.approx(rmax) + x_max = max(centers_x_cone[:,0]) + _RADIUS + x_min = min(centers_x_cone[:,0]) - _RADIUS + assert x_max < 1 or x_max == pytest.approx(1) + assert x_min > -1 or x_min == pytest.approx(-1) + + +def test_contained_y_cone(centers_y_cone): + """Make sure all spheres are entirely contained within the domain.""" + for p in centers_y_cone: + r = ((p[0] - 1)**2 + (p[2] - 2)**2)**0.5 + _RADIUS + rmax = p[1] + 1 - (0 - 1) + assert r < rmax or r == pytest.approx(rmax) + y_max = max(centers_y_cone[:,1]) + _RADIUS + y_min = min(centers_y_cone[:,1]) - _RADIUS + assert y_max < 1 or y_max == pytest.approx(1) + assert y_min > 0 or y_min == pytest.approx(-1) + + +def test_contained_z_cone(centers_z_cone): + """Make sure all spheres are entirely contained within the domain.""" + for p in centers_z_cone: + r = ((p[0] - 1)**2 + (p[1] - 2)**2)**0.5 + _RADIUS + rmax = p[2] + 1 - (0 - 1) + assert r < rmax or r == pytest.approx(rmax) + z_max = max(centers_z_cone[:,2]) + _RADIUS + z_min = min(centers_z_cone[:,2]) - _RADIUS + assert z_max < 1 or z_max == pytest.approx(1) + assert z_min > 0 or z_min == pytest.approx(-1) + def test_contained_sphere(centers_sphere): """Make sure all spheres are entirely contained within the domain."""