from abc import ABC, abstractproperty, abstractmethod
from collections import Counter, defaultdict
from collections.abc import Iterable
import copy
from heapq import heappush, heappop
import itertools
from math import pi, sin, cos, floor, log10, sqrt
from numbers import Real
import random
from random import uniform, gauss
import warnings
import numpy as np
import scipy.spatial
import openmc
from ..checkvalue import check_type
MAX_PF_RSP = 0.38 # maximum packing fraction for random sequential packing
MAX_PF_CRP = 0.64 # maximum packing fraction for close random packing
def _volume_sphere(r):
"""Return volume of a sphere of radius r"""
return 4/3 * pi * r**3
[docs]class TRISO(openmc.Cell):
"""Tristructural-isotopic (TRISO) micro fuel particle
Parameters
----------
outer_radius : float
Outer radius of TRISO particle
fill : openmc.Universe
Universe which contains all layers of the TRISO particle
center : Iterable of float
Cartesian coordinates of the center of the TRISO particle in cm
Attributes
----------
id : int
Unique identifier for the TRISO cell
name : str
Name of the TRISO cell
center : numpy.ndarray
Cartesian coordinates of the center of the TRISO particle in cm
fill : openmc.Universe
Universe that contains the TRISO layers
region : openmc.Region
Region of space within the TRISO particle
"""
def __init__(self, outer_radius, fill, center=(0., 0., 0.)):
self._surface = openmc.Sphere(r=outer_radius)
super().__init__(fill=fill, region=-self._surface)
self.center = np.asarray(center)
@property
def center(self):
return self._center
@center.setter
def center(self, center):
check_type('TRISO center', center, Iterable, Real)
self._surface.x0 = center[0]
self._surface.y0 = center[1]
self._surface.z0 = center[2]
self.translation = center
self._center = center
[docs] 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.region.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]))
class _Container(ABC):
"""Container in which to pack spheres.
Parameters
----------
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
----------
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.)
cell_length : list of float
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Constraint on where sphere center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, sphere_radius, center=(0., 0., 0.)):
self._cell_length = None
self._limits = None
self.sphere_radius = sphere_radius
self.center = center
@property
def sphere_radius(self):
return self._sphere_radius
@property
def center(self):
return self._center
@abstractproperty
def limits(self):
pass
@abstractproperty
def cell_length(self):
pass
@abstractproperty
def volume(self):
pass
@sphere_radius.setter
def sphere_radius(self, sphere_radius):
self._sphere_radius = float(sphere_radius)
self._limits = None
self._cell_length = None
@center.setter
def center(self, center):
self._center = center
def mesh_cell(self, p):
"""Calculate the index of the cell in a mesh overlaid on the domain in
which the given sphere center falls.
Parameters
----------
p : Iterable of float
Cartesian coordinates of sphere center.
Returns
-------
tuple of int
Indices of mesh cell.
"""
return tuple(int(p[i]/self.cell_length[i]) for i in range(3))
def nearby_mesh_cells(self, p):
"""Calculates the indices of all cells in a mesh overlaid on the domain
within one diameter of the given sphere.
Parameters
----------
p : Iterable of float
Cartesian coordinates of sphere center.
Returns
-------
list of tuple of int
Indices of mesh cells.
"""
d = 2*self.sphere_radius
r = [[a/self.cell_length[i] for a in [p[i]-d, p[i], p[i]+d]]
for i in range(3)]
return list(itertools.product(*({int(x) for x in y} for y in r)))
@abstractmethod
def from_region(self, region, sphere_radius):
"""Create a container to pack spheres in based on a region.
Parameters
----------
region : openmc.Region
Region to create container from.
sphere_radius : float
Outer radius of spheres.
"""
pass
@abstractmethod
def random_point(self):
"""Generate Cartesian coordinates of center of a sphere that is
contained entirely within the domain with uniform probability.
Returns
-------
list of float
Cartesian coordinates of sphere center.
"""
pass
@abstractmethod
def repel_spheres(self, p, q, d, d_new):
"""Move spheres p and q apart according to the following
transformation (accounting for boundary conditions on domain):
r_i^(n+1) = r_i^(n) + 1/2(d_out^(n+1) - d^(n))
r_j^(n+1) = r_j^(n) - 1/2(d_out^(n+1) - d^(n))
Parameters
----------
p, q : numpy.ndarray
Cartesian coordinates of sphere center.
d : float
distance between centers of spheres i and j.
d_new : float
final distance between centers of spheres i and j.
"""
pass
class _RectangularPrism(_Container):
"""Rectangular prism container in which to pack spheres.
Parameters
----------
width : float
Prism length along the x-axis
depth : float
Prism length along the y-axis
height : float
Prism length along the z-axis
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
----------
width : float
Prism length along the x-axis
depth : float
Prism length along the y-axis
height : float
Prism length along the z-axis
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.)
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 x-, y-, and z-direction where sphere
center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, width, depth, height, sphere_radius, center=(0., 0., 0.)):
super().__init__(sphere_radius, center)
self.width = width
self.depth = depth
self.height = height
@property
def width(self):
return self._width
@property
def depth(self):
return self._depth
@property
def height(self):
return self._height
@property
def limits(self):
if self._limits is None:
c = self.center
r = self.sphere_radius
x, y, z = self.width/2, self.depth/2, self.height/2
self._limits = [[c[0] - x + r, c[1] - y + r, c[2] - z + r],
[c[0] + x - r, c[1] + y - r, c[2] + z - r]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [self.width, self.depth, self.height]
self._cell_length = [x/int(x/(4*self.sphere_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return self.width*self.depth*self.height
@width.setter
def width(self, width):
self._width = float(width)
self._limits = None
self._cell_length = None
@depth.setter
def depth(self, depth):
self._depth = float(depth)
self._limits = None
self._cell_length = None
@height.setter
def height(self, height):
self._height = float(height)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
self._limits = limits
@classmethod
def from_region(self, region, sphere_radius):
check_type('region', region, openmc.Region)
# Assume the simplest case where the prism volume is the intersection
# of the half-spaces of six 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) != 6:
raise ValueError
# Sort half-spaces by surface type
px1, px2, py1, py2, pz1, pz2 = sorted(region, key=lambda x: x.surface.type)
# Make sure the region consists of the correct surfaces
if (not isinstance(px1.surface, openmc.XPlane) or
not isinstance(px2.surface, openmc.XPlane) or
not isinstance(py1.surface, openmc.YPlane) or
not isinstance(py2.surface, openmc.YPlane) or
not isinstance(pz1.surface, openmc.ZPlane) or
not isinstance(pz2.surface, openmc.ZPlane)):
raise ValueError
# Make sure the half-spaces are on the correct side of the surfaces
ll, ur = region.bounding_box
if any(x in ll or x in ur for x in (-np.inf, np.inf)):
raise ValueError
# Calculate the parameters for the container
width, depth, height = ur - ll
center = ll + [width/2, depth/2, height/2]
# The region is the volume of a rectangular prism, so create container
return _RectangularPrism(width, depth, height, sphere_radius, center)
def random_point(self):
ll, ul = self.limits
return [uniform(ll[0], ul[0]),
uniform(ll[1], ul[1]),
uniform(ll[2], ul[2])]
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 back along the
# surface normal until it is completely within the container if it
# overlaps the surface
p[:] = np.clip(p, self.limits[0], self.limits[1])
q[:] = np.clip(q, self.limits[0], self.limits[1])
class _Cylinder(_Container):
"""Cylindrical container in which to pack spheres.
Parameters
----------
length : float
Length of the cylindrical container.
radius : float
Radius of the cylindrical container.
axis : string
Axis along which the length of the cylinder 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 cylindrical container.
radius : float
Radius of the cylindrical container.
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 cylinder 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
Maximum radial distance and minimum and maximum distance in the
direction parallel to the axis where sphere center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, length, radius, axis, sphere_radius, center=(0., 0., 0.)):
super().__init__(sphere_radius, center)
self._shift = None
self.length = length
self.radius = radius
self.axis = axis
@property
def length(self):
return self._length
@property
def radius(self):
return self._radius
@property
def axis(self):
return self._axis
@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
self._limits = [[z0 - z + r], [z0 + z - r, self.radius - r]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
h = 4*self.sphere_radius
i, j, k = self.shift
self._cell_length = [None]*3
self._cell_length[i] = 2*self.radius/int(2*self.radius/h)
self._cell_length[j] = 2*self.radius/int(2*self.radius/h)
self._cell_length[k] = self.length/int(self.length/h)
return self._cell_length
@property
def volume(self):
return self.length*pi*self.radius**2
@length.setter
def length(self, length):
self._length = float(length)
self._limits = None
self._cell_length = None
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@axis.setter
def axis(self, axis):
self._axis = axis
self._shift = None
@limits.setter
def limits(self, limits):
self._limits = limits
@classmethod
def from_region(self, region, sphere_radius):
check_type('region', region, openmc.Region)
# Assume the simplest case where the cylinder volume is the
# intersection of the half-spaces of a cylinder 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 cylinder lies along
axis = region[0].surface.type[0]
# Make sure the region is composed of a cylinder and two planes on the
# same axis
count = Counter(node.surface.type for node in region)
if count[axis + '-cylinder'] != 1 or count[axis + '-plane'] != 2:
raise ValueError
# Sort the half-spaces by surface type
cyl, p1, p2 = sorted(region, key=lambda x: x.surface.type)
# Calculate the parameters for a cylinder 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, cyl.surface.y0, cyl.surface.z0)
# Calculate the parameters for a cylinder 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 = (cyl.surface.x0, p1.surface.y0 + length/2, cyl.surface.z0)
# Calculate the parameters for a cylinder along the z-axis
else:
if p1.surface.z0 > p2.surface.z0:
p1, p2 = p2, p1
length = p2.surface.z0 - p1.surface.z0
center = (cyl.surface.x0, cyl.surface.y0, p1.surface.z0 + length/2)
# Make sure the half-spaces are on the correct side of the surfaces
if cyl.side != '-' or p1.side != '+' or p2.side != '-':
raise ValueError
radius = cyl.surface.r
# The region is the volume of a cylinder, so create container
return _Cylinder(length, radius, axis, sphere_radius, center)
def random_point(self):
ll, ul = self.limits
r = sqrt(uniform(0, ul[1]**2))
t = uniform(0, 2*pi)
i, j, k = self.shift
p = [None]*3
p[i] = r*cos(t) + self.center[i]
p[j] = r*sin(t) + self.center[j]
p[k] = uniform(ll[0], ul[0])
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 back along the
# surface normal until it is completely within the container if it
# overlaps the surface
ll, ul = self.limits
c = self.center
i, j, k = self.shift
r = sqrt((p[i] - c[i])**2 + (p[j] - c[j])**2)
if r > ul[1]:
p[i] = (p[i] - c[i])*ul[1]/r + c[i]
p[j] = (p[j] - c[j])*ul[1]/r + c[j]
p[k] = np.clip(p[k], ll[0], ul[0])
r = sqrt((q[i] - c[i])**2 + (q[j] - c[j])**2)
if r > ul[1]:
q[i] = (q[i] - c[i])*ul[1]/r + c[i]
q[j] = (q[j] - c[j])*ul[1]/r + c[j]
q[k] = np.clip(q[k], ll[0], ul[0])
class _SphericalShell(_Container):
"""Spherical shell container in which to pack spheres.
Parameters
----------
radius : float
Outer radius of the spherical shell container.
inner_radius : float
Inner radius of the spherical shell container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
(0., 0., 0.)
Attributes
----------
radius : float
Outer radius of the spherical shell container.
inner_radius : float
Inner radius of the spherical shell container.
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.)
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 radial distance where sphere center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, radius, inner_radius, sphere_radius,
center=(0., 0., 0.)):
super().__init__(sphere_radius, center)
self.radius = radius
self.inner_radius = inner_radius
@property
def radius(self):
return self._radius
@property
def inner_radius(self):
return self._inner_radius
@property
def limits(self):
if self._limits is None:
r_max = self.radius - self.sphere_radius
if self.inner_radius == 0:
r_min = 0
else:
r_min = self.inner_radius + self.sphere_radius
self._limits = [[r_min], [r_max]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = 3*[2*self.radius]
self._cell_length = [x/int(x/(4*self.sphere_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return _volume_sphere(self.radius) - _volume_sphere(self.inner_radius)
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@inner_radius.setter
def inner_radius(self, inner_radius):
self._inner_radius = float(inner_radius)
self._limits = None
@limits.setter
def limits(self, limits):
self._limits = limits
@classmethod
def from_region(self, region, sphere_radius):
check_type('region', region, openmc.Region)
# First check if the region is the volume inside a sphere. Assume the
# simplest case where the sphere volume is the negative half-space of a
# sphere.
if (isinstance(region, openmc.Halfspace)
and isinstance(region.surface, openmc.Sphere)
and region.side == '-'):
# The region is the volume of a sphere, so create container
radius = region.surface.r
center = (region.surface.x0, region.surface.y0, region.surface.z0)
return _SphericalShell(radius, 0., sphere_radius, center)
# Next check for a spherical shell volume. Assume the simplest case
# where the spherical shell volume is the intersection of the
# half-spaces of two spheres.
if not isinstance(region, openmc.Intersection):
raise ValueError
if any(not isinstance(node, openmc.Halfspace) for node in region):
raise ValueError
if len(region) != 2:
raise ValueError
if any(not isinstance(node.surface, openmc.Sphere) for node in region):
raise ValueError
s1, s2 = sorted(region, key=lambda x: x.surface.r)
radius = s2.surface.r
inner_radius = s1.surface.r
center = (s1.surface.x0, s1.surface.y0, s1.surface.z0)
if center != (s2.surface.x0, s2.surface.y0, s2.surface.z0):
raise ValueError
if s1.side != '+' or s2.side != '-':
raise ValueError
# The region is the volume of a spherical shell, so create container
return _SphericalShell(radius, inner_radius, sphere_radius, center)
def random_point(self):
c = self.center
ll, ul = self.limits
x, y, z = (gauss(0, 1), gauss(0, 1), gauss(0, 1))
r = (uniform(ll[0]**3, ul[0]**3)**(1/3)/sqrt(x**2 + y**2 + z**2))
return [r*x + c[0], r*y + c[1], r*z + c[2]]
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 back along the
# surface normal until it is completely within the container if it
# overlaps the surface
c = self.center
ll, ul = self.limits
r = sqrt((p[0] - c[0])**2 + (p[1] - c[1])**2 + (p[2] - c[2])**2)
if r > ul[0]:
p[:] = (p - c)*ul[0]/r + c
elif r < ll[0]:
p[:] = (p - c)*ll[0]/r + c
r = sqrt((q[0] - c[0])**2 + (q[1] - c[1])**2 + (q[2] - c[2])**2)
if r > ul[0]:
q[:] = (q - c)*ul[0]/r + c
elif r < ll[0]:
q[:] = (q - c)*ll[0]/r + c
[docs]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 triso_locations:
# Create copy of TRISO particle with materials preserved and
# different cell/surface IDs
t_copy = copy.copy(t)
t_copy.id = None
t_copy.fill = t.fill
t_copy._surface = openmc.Sphere(r=t._surface.r,
x0=t._surface.x0,
y0=t._surface.y0,
z0=t._surface.z0)
t_copy.region = -t_copy._surface
triso_locations[idx].append(t_copy)
else:
warnings.warn('TRISO particle is partially or completely '
'outside of the lattice.')
# 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.region 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_cell(t)
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
def _random_sequential_pack(domain, num_spheres):
"""Random sequential packing of spheres within a container.
Parameters
----------
domain : openmc.model._Container
Container in which to pack spheres.
num_spheres : int
Number of spheres to pack.
Returns
------
numpy.ndarray
Cartesian coordinates of centers of spheres.
"""
sqd = (2*domain.sphere_radius)**2
spheres = []
mesh = defaultdict(list)
for i in range(num_spheres):
# Randomly sample new center coordinates while there are any overlaps
while True:
p = domain.random_point()
idx = domain.mesh_cell(p)
if any((p[0]-q[0])**2 + (p[1]-q[1])**2 + (p[2]-q[2])**2 < sqd
for q in mesh[idx]):
continue
else:
break
spheres.append(p)
for idx in domain.nearby_mesh_cells(p):
mesh[idx].append(p)
return np.array(spheres)
def _close_random_pack(domain, spheres, contraction_rate):
"""Close random packing of spheres using the Jodrey-Tory algorithm.
Parameters
----------
domain : openmc.model._Container
Container in which to pack spheres.
spheres : numpy.ndarray
Initial Cartesian coordinates of centers of spheres.
contraction_rate : float
Contraction rate of outer diameter.
"""
def add_rod(d, i, j):
"""Add a new rod to the priority queue.
Parameters
----------
d : float
distance between centers of spheres i and j.
i, j : int
Index of spheres in spheres array.
"""
rod = [d, i, j]
rods_map[i] = (j, rod)
rods_map[j] = (i, rod)
heappush(rods, rod)
def remove_rod(i):
"""Mark the rod containing sphere i as removed.
Parameters
----------
i : int
Index of sphere in spheres array.
"""
if i in rods_map:
j, rod = rods_map.pop(i)
del rods_map[j]
rod[1] = removed
rod[2] = removed
def pop_rod():
"""Remove and return the shortest rod.
Returns
-------
d : float
distance between centers of spheres i and j.
i, j : int
Index of spheres in spheres array.
"""
while rods:
d, i, j = heappop(rods)
if i != removed and j != removed:
del rods_map[i]
del rods_map[j]
return d, i, j
return None, None, None
def create_rod_list():
"""Generate sorted list of rods (distances between sphere centers).
Rods are arranged in a heap where each element contains the rod length
and the sphere indices. A rod between spheres p and q is only
included if the distance between p and q could not be changed by the
elimination of a greater overlap, i.e. q has no nearer neighbors than p.
A mapping of sphere ids to rods is maintained in 'rods_map'. Each key
in the dict is the id of a sphere that is in the rod list, and the
value is the id of its nearest neighbor and the rod that contains them.
The dict is used to find rods in the priority queue and to mark removed
rods so rods can be "removed" without breaking the heap structure
invariant.
"""
# Create KD tree for quick nearest neighbor search
tree = scipy.spatial.cKDTree(spheres)
# Find distance to nearest neighbor and index of nearest neighbor for
# all spheres
d, n = tree.query(spheres, k=2)
d = d[:, 1]
n = n[:, 1]
# Array of sphere indices, indices of nearest neighbors, and
# distances to nearest neighbors
a = np.vstack((list(range(n.size)), n, d)).T
# Sort along second column and swap first and second columns to create
# array of nearest neighbor indices, indices of spheres they are
# nearest neighbors of, and distances between them
b = a[a[:, 1].argsort()]
b[:, [0, 1]] = b[:, [1, 0]]
# Find the intersection between 'a' and 'b': a list of spheres who
# are each other's nearest neighbors and the distance between them
r = list({tuple(x) for x in a} & {tuple(x) for x in b})
# Remove duplicate rods and sort by distance
r = map(list, set([(x[2], int(min(x[0:2])), int(max(x[0:2])))
for x in r]))
# Clear priority queue and add rods
del rods[:]
rods_map.clear()
for d, i, j in r:
if d < outer_diameter and not np.isclose(d, outer_diameter, atol=1.0e-14):
add_rod(d, i, j)
def update_mesh(i):
"""Update which mesh cells the sphere is in based on new sphere
center coordinates.
'mesh'/'mesh_map' is a two way dictionary used to look up which
spheres are located within one diameter of a given mesh cell and
which mesh cells a given sphere center is within one diameter of.
This is used to speed up the nearest neighbor search.
Parameters
----------
i : int
Index of sphere in spheres array.
"""
# Determine which mesh cells the sphere is in and remove the
# sphere id from those cells
for idx in mesh_map[i]:
mesh[idx].remove(i)
del mesh_map[i]
# Determine which mesh cells are within one diameter of sphere's
# center and add this sphere to the list of spheres in those cells
for idx in domain.nearby_mesh_cells(spheres[i]):
mesh[idx].add(i)
mesh_map[i].add(idx)
def reduce_outer_diameter():
"""Reduce the outer diameter so that at the (i+1)-st iteration it is:
d_out^(i+1) = d_out^(i) - (1/2)^(j) * d_out0 * k / n,
where k is the contraction rate, n is the number of spheres, and
j = floor(-log10(pf_out - pf_in)).
Returns
-------
float
New outer diameter
"""
inner_pf = _volume_sphere(inner_diameter/2)*num_spheres / domain.volume
outer_pf = _volume_sphere(outer_diameter/2)*num_spheres / domain.volume
j = floor(-log10(outer_pf - inner_pf))
return (outer_diameter - 0.5**j * contraction_rate *
initial_outer_diameter / num_spheres)
def nearest(i):
"""Find index of nearest neighbor of sphere i.
Parameters
----------
i : int
Index in spheres array of sphere for which to find nearest
neighbor.
Returns
-------
int
Index in spheres array of nearest neighbor of i
float
distance between i and nearest neighbor.
"""
# Need the second nearest neighbor of i since the nearest neighbor
# will be itself. Using argpartition, the k-th nearest neighbor is
# placed at index k.
idx = list(mesh[domain.mesh_cell(spheres[i])])
dists = scipy.spatial.distance.cdist([spheres[i]], spheres[idx])[0]
if dists.size > 1:
j = dists.argpartition(1)[1]
return idx[j], dists[j]
else:
return None, None
def update_rod_list(i):
"""Update the rod list with the new nearest neighbors of sphere since
its overlap was eliminated.
Parameters
----------
i : int
Index of sphere in spheres array.
"""
# If the nearest neighbor k of sphere i has no nearer neighbors,
# remove the rod currently containing k from the rod list and add rod
# k-i, keeping the rod list sorted
k, d_ik = nearest(i)
if (k and nearest(k)[0] == i and d_ik < outer_diameter
and not np.isclose(d, outer_diameter, atol=1.0e-14)):
remove_rod(k)
add_rod(d_ik, i, k)
num_spheres = len(spheres)
diameter = 2*domain.sphere_radius
# Flag for marking rods that have been removed from priority queue
removed = -1
# Outer diameter initially set to arbitrary value that yields pf of 1
initial_outer_diameter = 2*(domain.volume/(num_spheres*4/3*pi))**(1/3)
# Inner and outer diameter of spheres will change during packing
outer_diameter = initial_outer_diameter
inner_diameter = 0.
# List of rods arranged in a heap and mapping of sphere ids to rods
rods = []
rods_map = {}
# Initialize two-way dictionary that identifies which spheres are near a
# given mesh cell and which mesh cells a sphere is near
mesh = defaultdict(set)
mesh_map = defaultdict(set)
for i in range(num_spheres):
for idx in domain.nearby_mesh_cells(spheres[i]):
mesh[idx].add(i)
mesh_map[i].add(idx)
while True:
# Rebuild the sorted list of rods according to the current sphere
# configuration
create_rod_list()
# Set the inner diameter to the shortest center-to-center distance
# between any two spheres
if rods:
inner_diameter = rods[0][0]
# Reached the desired sphere radius
if inner_diameter >= diameter:
break
# The algorithm converged before reaching the desired sphere radius.
# This can happen when the desired packing fraction is close to the
# packing fraction limit. The packing fraction is a random variable
# that is determined by the sphere locations and the contraction
# rate. A higher packing fraction can be achieved with a smaller
# contraction rate, though at the cost of a longer simulation time --
# the number of iterations needed to remove all overlaps is inversely
# proportional to the contraction rate.
if inner_diameter >= outer_diameter or not rods:
warnings.warn('Close random pack converged before reaching true '
'sphere radius; some spheres may overlap. Try '
'reducing contraction rate or packing fraction.')
break
while True:
d, i, j = pop_rod()
if not d:
break
outer_diameter = reduce_outer_diameter()
domain.repel_spheres(spheres[i], spheres[j], d, outer_diameter)
update_mesh(i)
update_mesh(j)
update_rod_list(i)
update_rod_list(j)
if not rods:
break
inner_diameter = rods[0][0]
if inner_diameter >= diameter or inner_diameter >= outer_diameter:
break
[docs]def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3,
contraction_rate=1.e-3, seed=1):
"""Generate a random, non-overlapping configuration of spheres within a
container.
Parameters
----------
radius : float
Outer radius of spheres.
region : openmc.Region
Container in which the spheres are packed. Supported shapes are
rectangular prism, cylinder, sphere, and spherical shell.
pf : float
Packing fraction of the spheres. One of 'pf' and 'num_spheres' must
be specified; the other will be calculated. If both are specified, 'pf'
takes precedence over 'num_spheres'.
num_spheres : int
Number of spheres to pack in the domain. One of 'num_spheres' and 'pf'
must be specified; the other will be calculated.
initial_pf : float, optional
Packing fraction used to initialize the configuration of spheres in
the domain. Default value is 0.3. It is not recommended to set the
initial packing fraction much higher than 0.3 as the random sequential
packing algorithm becomes prohibitively slow as it approaches its limit
(~0.38).
contraction_rate : float, optional
Contraction rate of the outer diameter. Higher packing fractions can be
reached using a smaller contraction rate, but the algorithm will take
longer to converge.
seed : int, optional
RNG seed.
Returns
------
numpy.ndarray
Cartesian coordinates of sphere centers.
Notes
-----
The sphere configuration is generated using a combination of random
sequential packing (RSP) and close random packing (CRP). RSP performs
better than CRP for lower packing fractions (pf), but it becomes
prohibitively slow as it approaches its packing limit (~0.38). CRP can
achieve higher pf of up to ~0.64 and scales better with increasing pf.
If the desired pf is below some threshold for which RSP will be faster than
CRP ('initial_packing_fraction'), only RSP is used. If a higher pf is
required, spheres with a radius smaller than the desired final radius
(and therefore with a smaller pf) are initialized within the domain using
RSP. This initial configuration of spheres is then used as a starting
point for CRP using Jodrey and Tory's algorithm [1]_.
In RSP, sphere centers are placed one by one at random, and placement
attempts for a sphere are made until the sphere is not overlapping any
others. This implementation of the algorithm uses a mesh over the domain
to speed up the nearest neighbor search by only searching for a sphere's
neighbors within that mesh cell.
In CRP, each sphere is assigned two diameters, an inner and an outer,
which approach each other during the simulation. The inner diameter,
defined as the minimum center-to-center distance, is the true diameter of
the spheres and defines the pf. At each iteration the worst overlap
between spheres based on outer diameter is eliminated by moving the
spheres apart along the line joining their centers. Iterations continue
until the two diameters converge or until the desired pf is reached.
References
----------
.. [1] W. S. Jodrey and E. M. Tory, "Computer simulation of close random
packing of equal spheres", Phys. Rev. A 32 (1985) 2347-2351.
"""
# Seed RNG
random.seed(seed)
# Create container with the correct shape based on the supplied region
domain = None
for cls in _Container.__subclasses__():
try:
domain = cls.from_region(region, radius)
except ValueError:
pass
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))
# Determine the packing fraction/number of spheres
volume = _volume_sphere(radius)
if pf is None and num_spheres is None:
raise ValueError('`pf` or `num_spheres` must be specified.')
elif pf is None:
num_spheres = int(num_spheres)
pf = volume*num_spheres/domain.volume
else:
pf = float(pf)
num_spheres = int(pf*domain.volume//volume)
# Make sure initial packing fraction is less than packing fraction
if initial_pf > pf:
initial_pf = pf
# Check packing fraction for close random packing
if pf > MAX_PF_CRP:
raise ValueError('Packing fraction {0} is greater than the limit for '
'close random packing, {1}'.format(pf, MAX_PF_CRP))
# Check packing fraction for random sequential packing
if initial_pf > MAX_PF_RSP:
raise ValueError('Initial packing fraction {0} is greater than the '
'limit for random sequential packing, '
'{1}'.format(initial_pf, MAX_PF_RSP))
# Calculate the sphere radius used in the initial random sequential
# packing from the initial packing fraction
initial_radius = (3/4*initial_pf*domain.volume/(pi*num_spheres))**(1/3)
domain.sphere_radius = initial_radius
# 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]]]
# Generate non-overlapping spheres for an initial inner radius using
# random sequential packing algorithm
spheres = _random_sequential_pack(domain, num_spheres)
# Use the sphere configuration produced in random sequential packing as a
# starting point for close random pack with the desired final sphere
# radius
if initial_pf != pf:
domain.sphere_radius = radius
_close_random_pack(domain, spheres, contraction_rate)
return spheres