From d43e95def735f4d26f1ffde7e65b837cb7fb17d3 Mon Sep 17 00:00:00 2001 From: ZoeRichter Date: Wed, 14 Jan 2026 16:17:49 -0600 Subject: [PATCH 1/5] testing cone random pack --- openmc/model/triso.py | 264 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 264 insertions(+) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index f344b8edd42..aeadfe502ed 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 = seld.radius_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, + [(r_major-r_minor)/(2*z), + (r_major-r_minor)/(2*z)*(z - z0) + r_minor - r]]] + 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*r_avg/int(2*r_avg/h) + self._cell_length[j] = 2*r_avg/int(2*r_avg/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 > rlp: + 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. From d82dad814f40a9b4e4df812c6118ce4e29f61171 Mon Sep 17 00:00:00 2001 From: ZoeRichter Date: Thu, 15 Jan 2026 11:53:45 -0600 Subject: [PATCH 2/5] fixed rl equation for cones --- openmc/model/triso.py | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index aeadfe502ed..88bc6bbc78f 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -753,14 +753,15 @@ def limits(self): z = self.length/2 r = self.sphere_radius r_major = self.radius_major - r_minor = seld.radius_minor + 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, - [(r_major-r_minor)/(2*z), - (r_major-r_minor)/(2*z)*(z - z0) + r_minor - r]]] + [1/m, + r_minor - (z0 - z)/m - (r/m)*(m**2 + 1)**0.5]]] return self._limits @limits.setter @@ -771,11 +772,11 @@ def limits(self, limits): 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) + #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*r_avg/int(2*r_avg/h) - self._cell_length[j] = 2*r_avg/int(2*r_avg/h) + 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 @@ -886,7 +887,7 @@ def repel_spheres(self, p, q, d, d_new): 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 > rlp: + if r > rl: p[i] = (p[i] - c[i])*rl/r + c[i] p[j] = (p[j] - c[j])*rl/r + c[j] @@ -897,7 +898,6 @@ def repel_spheres(self, p, q, d, d_new): 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. @@ -1558,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) @@ -1593,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 From e17d4ab24c5bbce475a0174983c32456c3885ffd Mon Sep 17 00:00:00 2001 From: ZoeRichter Date: Thu, 15 Jan 2026 14:57:04 -0600 Subject: [PATCH 3/5] added unit tests --- tests/unit_tests/test_model_triso.py | 67 ++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index c1f8c039e57..d1fc1f8a65b 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': 1/3*pi*(1**2 + 1*2 + 2**2)}, + {'shape': 'y_cone', 'volume': 1/3*pi*(1**2 + 1*2 + 2**2)}, + {'shape': 'z_cone', 'volume': 1/3*pi*(1**2 + 1*2 + 2**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=-1, y0=1, z0=2, r2=1) + min_x = openmc.XPlane(0) + 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=-1, z0=2, r2=1) + min_y = openmc.YPlane(0) + 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=-1, r2=1) + min_z = openmc.ZPlane(0) + 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 + _RADIUS + 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 > 0 or x_min == pytest.approx(0) + + +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 + _RADIUS + 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(0) + + +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 + _RADIUS + 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(0) + def test_contained_sphere(centers_sphere): """Make sure all spheres are entirely contained within the domain.""" From ba3388c5d979e8bc653b6e68f5e858ff1fac3e7d Mon Sep 17 00:00:00 2001 From: ZoeRichter Date: Fri, 16 Jan 2026 15:40:53 -0600 Subject: [PATCH 4/5] rewrote limit equation to match theory eqn better, fixed sign error --- openmc/model/triso.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 88bc6bbc78f..0fc12ac4118 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -761,7 +761,7 @@ def limits(self): self._limits = [[z0 - z + r], [z0 + z - r, [1/m, - r_minor - (z0 - z)/m - (r/m)*(m**2 + 1)**0.5]]] + -(r/m)*(m**2 + 1)**0.5 - (z0-z)/m + r_minor]]] return self._limits @limits.setter From 9cd86e476be18bd06c8a9aa2bcf98ecc7c1d7608 Mon Sep 17 00:00:00 2001 From: ZoeRichter Date: Fri, 16 Jan 2026 15:54:22 -0600 Subject: [PATCH 5/5] changed unit test to catch bug --- tests/unit_tests/test_model_triso.py | 30 ++++++++++++++-------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index d1fc1f8a65b..4a16e264bab 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -17,9 +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': 1/3*pi*(1**2 + 1*2 + 2**2)}, - {'shape': 'y_cone', 'volume': 1/3*pi*(1**2 + 1*2 + 2**2)}, - {'shape': 'z_cone', 'volume': 1/3*pi*(1**2 + 1*2 + 2**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)} ] @@ -79,8 +79,8 @@ def centers_z_cylinder(): @pytest.fixture(scope='module') def centers_x_cone(): - cone = openmc.XCone(x0=-1, y0=1, z0=2, r2=1) - min_x = openmc.XPlane(0) + 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, @@ -89,8 +89,8 @@ def centers_x_cone(): @pytest.fixture(scope='module') def centers_y_cone(): - cone = openmc.YCone(x0=1, y0=-1, z0=2, r2=1) - min_y = openmc.YPlane(0) + 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, @@ -99,8 +99,8 @@ def centers_y_cone(): @pytest.fixture(scope='module') def centers_z_cone(): - cone = openmc.ZCone(x0=1, y0=2, z0=-1, r2=1) - min_z = openmc.ZPlane(0) + 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, @@ -189,36 +189,36 @@ 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 + _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 > 0 or x_min == pytest.approx(0) + 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 + _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(0) + 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 + _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(0) + assert z_min > 0 or z_min == pytest.approx(-1) def test_contained_sphere(centers_sphere):