diff --git a/openmc/model/triso.py b/openmc/model/triso.py index f344b8edd42..896e44a8558 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -14,7 +14,7 @@ import scipy.spatial import openmc -from ..checkvalue import check_type +from ..checkvalue import check_greater_than, check_type MAX_PF_RSP = 0.38 # maximum packing fraction for random sequential packing @@ -802,6 +802,295 @@ def repel_spheres(self, p, q, d, d_new): q[:] = (q - c)*ll[0]/r + c +class _RegionContainer(_Container): + """Generic region container in which to pack spheres. + + Parameters + ---------- + region : openmc.Region + Region defining the container boundary. + bounding_box : openmc.BoundingBox + Finite bounding box enclosing the region. + volume : float + Volume of the region. + sphere_radius : float + Radius of spheres to be packed in container. + + Attributes + ---------- + region : openmc.Region + Region defining the container boundary. + bounding_box : openmc.BoundingBox + Finite bounding box enclosing the region. + volume : float + Volume of the region. + sphere_radius : float + Radius of spheres to be packed in container. + center : list of float + Cartesian coordinates of the center of the container bounding box. + 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. + + """ + + _DIRECTION_COUNT = 50 + _SAFETY_FACTOR = 1.001 + _MAX_POINT_ATTEMPTS = 10000 + _MAX_BACKTRACKS = 20 + + def __init__(self, region, bounding_box, volume, sphere_radius): + super().__init__(sphere_radius, center=bounding_box.center) + self.region = region + self.bounding_box = bounding_box + self._volume = float(volume) + self._boundary_radius = None + self._analytic_surfaces = self._build_analytic_surfaces() + self._directions = None + + @property + def region(self): + return self._region + + @region.setter + def region(self, region): + check_type('region', region, openmc.Region) + self._region = region + + @property + def bounding_box(self): + return self._bounding_box + + @bounding_box.setter + def bounding_box(self, bounding_box): + check_type('bounding_box', bounding_box, openmc.BoundingBox) + self._bounding_box = bounding_box + self._limits = None + self._cell_length = None + + @property + def boundary_radius(self): + if self._boundary_radius is None: + return self.sphere_radius + return self._boundary_radius + + @boundary_radius.setter + def boundary_radius(self, boundary_radius): + self._boundary_radius = float(boundary_radius) + + @property + def limits(self): + if self._limits is None: + r = self.sphere_radius + ll = self.bounding_box.lower_left + r + ur = self.bounding_box.upper_right - r + self._limits = [ll, ur] + return self._limits + + @limits.setter + def limits(self, limits): + self._limits = limits + + @property + def cell_length(self): + if self._cell_length is None: + mesh_length = self.bounding_box.width + 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._volume + + @classmethod + def from_region(self, region, sphere_radius, bounding_box=None, volume=None): + check_type('region', region, openmc.Region) + + if volume is None: + raise ValueError('`volume` must be provided for generic regions.') + check_type('volume', volume, Real) + check_greater_than('volume', volume, 0.0) + + if bounding_box is None: + bounding_box = region.bounding_box + check_type('bounding_box', bounding_box, openmc.BoundingBox) + if (np.isinf(bounding_box.lower_left).any() or + np.isinf(bounding_box.upper_right).any()): + raise ValueError('`bounding_box` must be finite for generic regions.') + + return _RegionContainer(region, bounding_box, volume, sphere_radius) + + def random_point(self): + ll, ul = self.limits + for _ in range(self._MAX_POINT_ATTEMPTS): + p = np.array([uniform(ll[0], ul[0]), + uniform(ll[1], ul[1]), + uniform(ll[2], ul[2])]) + if self._is_valid_center(p): + return p.tolist() + raise ValueError('Failed to sample a valid point in the region. ' + 'Consider providing a tighter bounding_box or ' + 'verifying the region is closed.') + + def repel_spheres(self, p, q, d, d_new): + s = (d_new - d)/2 + + v = (p - q)/d + p_old = p.copy() + q_old = q.copy() + p_target = p_old + s*v + q_target = q_old - s*v + if (self._is_valid_center(p_target) and + self._is_valid_center(q_target)): + p[:] = p_target + q[:] = q_target + return + + self._backtrack_pair(p, q, p_old, q_old, p_target, q_target) + + def _backtrack_pair(self, p, q, p_old, q_old, p_target, q_target): + for i in range(self._MAX_BACKTRACKS): + alpha = 0.5**(i + 1) + p_candidate = p_old + alpha*(p_target - p_old) + q_candidate = q_old + alpha*(q_target - q_old) + if (self._is_valid_center(p_candidate) and + self._is_valid_center(q_candidate)): + p[:] = p_candidate + q[:] = q_candidate + return + p[:] = p_old + q[:] = q_old + + def _is_valid_center(self, point): + point = np.asarray(point, dtype=float) + boundary_radius = self.boundary_radius + if self._analytic_surfaces is not None: + for entry in self._analytic_surfaces: + if self._analytic_clearance(entry, point) < boundary_radius: + return False + return True + + for direction in self._direction_vectors(): + if (point + self._SAFETY_FACTOR*boundary_radius*direction + not in self.region): + return False + return True + + def _direction_vectors(self): + if self._directions is None: + self._directions = self._fibonacci_directions(self._DIRECTION_COUNT) + return self._directions + + @staticmethod + def _fibonacci_directions(count): + directions = [] + offset = 2.0/count + increment = pi*(3.0 - sqrt(5.0)) + for i in range(count): + y = ((i*offset) - 1) + (offset/2) + r = sqrt(max(0.0, 1 - y*y)) + phi = i*increment + x = cos(phi)*r + z = sin(phi)*r + directions.append(np.array([x, y, z])) + return directions + + def _build_analytic_surfaces(self): + if isinstance(self.region, openmc.Halfspace): + nodes = [self.region] + elif (isinstance(self.region, openmc.Intersection) and + all(isinstance(node, openmc.Halfspace) for node in self.region)): + nodes = list(self.region) + else: + return None + + entries = [] + for node in nodes: + surface = node.surface + side = node.side + + if isinstance(surface, (openmc.Plane, openmc.XPlane, openmc.YPlane, + openmc.ZPlane)): + a, b, c = surface.a, surface.b, surface.c + norm = sqrt(a*a + b*b + c*c) + if np.isclose(norm, 0.0): + return None + entries.append(('plane', side, surface, norm)) + elif isinstance(surface, openmc.Sphere): + center = np.array([surface.x0, surface.y0, surface.z0]) + entries.append(('sphere', side, center, surface.r)) + elif isinstance(surface, openmc.XCylinder): + entries.append(('xcylinder', side, surface.y0, surface.z0, + surface.r)) + elif isinstance(surface, openmc.YCylinder): + entries.append(('ycylinder', side, surface.x0, surface.z0, + surface.r)) + elif isinstance(surface, openmc.ZCylinder): + entries.append(('zcylinder', side, surface.x0, surface.y0, + surface.r)) + elif isinstance(surface, (openmc.Cone, openmc.XCone, openmc.YCone, + openmc.ZCone)): + axis = np.array([surface.dx, surface.dy, surface.dz], + dtype=float) + axis_norm = np.linalg.norm(axis) + if np.isclose(axis_norm, 0.0): + return None + axis = axis/axis_norm + apex = np.array([surface.x0, surface.y0, surface.z0]) + k = sqrt(surface.r2) + entries.append(('cone', side, apex, axis, k)) + else: + return None + + return entries + + def _analytic_clearance(self, entry, point): + kind = entry[0] + side = entry[1] + + if kind == 'plane': + surface = entry[2] + norm = entry[3] + value = surface.evaluate(point) + return (value if side == '+' else -value)/norm + + if kind == 'sphere': + center = entry[2] + radius = entry[3] + r = np.linalg.norm(point - center) + return radius - r if side == '-' else r - radius + + if kind == 'xcylinder': + y0, z0, radius = entry[2], entry[3], entry[4] + r = sqrt((point[1] - y0)**2 + (point[2] - z0)**2) + return radius - r if side == '-' else r - radius + + if kind == 'ycylinder': + x0, z0, radius = entry[2], entry[3], entry[4] + r = sqrt((point[0] - x0)**2 + (point[2] - z0)**2) + return radius - r if side == '-' else r - radius + + if kind == 'zcylinder': + x0, y0, radius = entry[2], entry[3], entry[4] + r = sqrt((point[0] - x0)**2 + (point[1] - y0)**2) + return radius - r if side == '-' else r - radius + + apex = entry[2] + axis = entry[3] + k = entry[4] + w = point - apex + a = np.dot(w, axis) + b = np.linalg.norm(w - a*axis) + denom = sqrt(1 + k*k) + if side == '-': + return (k*abs(a) - b)/denom + return (b - k*abs(a))/denom + + def create_triso_lattice(trisos, lower_left, pitch, shape, background): """Create a lattice containing TRISO particles for optimized tracking. @@ -1208,7 +1497,8 @@ def update_rod_list(i): def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, - contraction_rate=1.e-3, seed=None): + contraction_rate=1.e-3, seed=None, bounding_box=None, + volume=None): """Generate a random, non-overlapping configuration of spheres within a container. @@ -1218,7 +1508,9 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, Outer radius of spheres. region : openmc.Region Container in which the spheres are packed. Supported shapes are - rectangular prism, cylinder, sphere, and spherical shell. + rectangular prism, cylinder, sphere, and spherical shell. Other closed + regions are supported when `bounding_box` is finite and `volume` is + provided. 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' @@ -1238,6 +1530,12 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, longer to converge. seed : int, optional Pseudorandom number generator seed passed to :func:`random.seed` + bounding_box : openmc.BoundingBox, optional + Bounding box used to sample sphere centers for generic regions. This + is required when the region's bounding box is unbounded. + volume : float, optional + Volume of the region. This is required for generic regions and is + ignored for analytic container shapes. Returns ------ @@ -1285,27 +1583,37 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, # Create container with the correct shape based on the supplied region domain = None - for cls in _Container.__subclasses__(): + for cls in (_RectangularPrism, _Cylinder, _SphericalShell): try: domain = cls.from_region(region, radius) + break 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)) + try: + domain = _RegionContainer.from_region( + region, radius, bounding_box=bounding_box, volume=volume) + except ValueError as exc: + raise ValueError( + f'Could not create container from region {region}: {exc}' + ) from exc # Determine the packing fraction/number of spheres - volume = _volume_sphere(radius) + sphere_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 + pf = sphere_volume*num_spheres/domain.volume else: pf = float(pf) - num_spheres = int(pf*domain.volume//volume) + num_spheres = int(pf*domain.volume//sphere_volume) + + pf_actual = sphere_volume*num_spheres/domain.volume + if isinstance(domain, _RegionContainer) and pf <= MAX_PF_RSP: + pf = pf_actual + initial_pf = pf # Make sure initial packing fraction is less than packing fraction if initial_pf > pf: @@ -1332,6 +1640,9 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, domain.limits = [[x - initial_radius + radius for x in domain.limits[0]], [x + initial_radius - radius for x in domain.limits[1]]] + if isinstance(domain, _RegionContainer): + domain.boundary_radius = radius + # Generate non-overlapping spheres for an initial inner radius using # random sequential packing algorithm spheres = _random_sequential_pack(domain, num_spheres) diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index c1f8c039e57..243ca02d932 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -from math import pi +from math import pi, sqrt import numpy as np from numpy.linalg import norm @@ -12,13 +12,19 @@ _RADIUS = 0.1 _PACKING_FRACTION = 0.35 +_TETRA_SIZE = 1 +_TETRA_VOLUME = _TETRA_SIZE**3/6 +_TETRA_NUM_SPHERES = 8 +_TETRA_PACKING_FRACTION = 0.2 _PARAMS = [ {'shape': 'rectangular_prism', 'volume': 1**3}, {'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': 'sphere', 'volume': 4/3*pi*1**3}, - {'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)} + {'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)}, + {'shape': 'tetrahedron', 'volume': _TETRA_VOLUME, + 'pf': _TETRA_PACKING_FRACTION} ] @@ -92,6 +98,20 @@ def centers_spherical_shell(): pf=_PACKING_FRACTION, initial_pf=0.2) +@pytest.fixture(scope='module') +def centers_tetrahedron(): + min_x = openmc.XPlane(0) + min_y = openmc.YPlane(0) + min_z = openmc.ZPlane(0) + max_plane = openmc.Plane(a=1, b=1, c=1, d=_TETRA_SIZE) + region = +min_x & +min_y & +min_z & -max_plane + bounding_box = openmc.BoundingBox((0, 0, 0), + (_TETRA_SIZE, _TETRA_SIZE, _TETRA_SIZE)) + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + num_spheres=_TETRA_NUM_SPHERES, volume=_TETRA_VOLUME, + bounding_box=bounding_box) + + @pytest.fixture(scope='module') def triso_universe(): sphere = openmc.Sphere(r=_RADIUS) @@ -170,10 +190,24 @@ def test_contained_spherical_shell(centers_spherical_shell): assert r_min > 0.5 or r_min == pytest.approx(0.5) +def test_contained_tetrahedron(centers_tetrahedron): + """Make sure all spheres are entirely contained within the domain.""" + x_min = min(centers_tetrahedron[:, 0]) - _RADIUS + y_min = min(centers_tetrahedron[:, 1]) - _RADIUS + z_min = min(centers_tetrahedron[:, 2]) - _RADIUS + limit = _TETRA_SIZE - _RADIUS*sqrt(3) + sum_max = max(centers_tetrahedron.sum(axis=1)) + assert x_min > 0 or x_min == pytest.approx(0) + assert y_min > 0 or y_min == pytest.approx(0) + assert z_min > 0 or z_min == pytest.approx(0) + assert sum_max < limit or sum_max == pytest.approx(limit) + + def test_packing_fraction(container, centers): """Check that the actual PF is close to the requested PF.""" pf = len(centers) * 4/3 * pi *_RADIUS**3 / container['volume'] - assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2) + target_pf = container.get('pf', _PACKING_FRACTION) + assert pf == pytest.approx(target_pf, rel=1e-2) def test_num_spheres():