From 98b6dc15a94e6911e73c40fd3d0b9a9163aa4513 Mon Sep 17 00:00:00 2001
From: amandalund <alund@anl.gov>
Date: Sun, 25 Mar 2018 13:47:47 -0500
Subject: [PATCH 1/4] Fix that ensures TRISO particles are completely contained
 within the domain

---
 openmc/model/triso.py | 95 +++++++++++++++++++++++++++----------------
 1 file changed, 59 insertions(+), 36 deletions(-)

diff --git a/openmc/model/triso.py b/openmc/model/triso.py
index 5fe51b6644..b5ee74f90b 100644
--- a/openmc/model/triso.py
+++ b/openmc/model/triso.py
@@ -113,8 +113,7 @@ class _Domain(metaclass=ABCMeta):
         Length in x-, y-, and z- directions of each cell in mesh overlaid on
         domain.
     limits : list of float
-        Minimum and maximum position in x-, y-, and z-directions where particle
-        center can be placed.
+        Constraint on where particle center can be placed.
     volume : float
         Volume of the container.
 
@@ -158,8 +157,6 @@ class _Domain(metaclass=ABCMeta):
             raise ValueError('Unable to set domain center to {} since it must '
                              'be of length 3'.format(center))
         self._center = [float(x) for x in center]
-        self._limits = None
-        self._cell_length = None
 
     def mesh_cell(self, p):
         """Calculate the index of the cell in a mesh overlaid on the domain in
@@ -211,6 +208,19 @@ class _Domain(metaclass=ABCMeta):
         """
         pass
 
+    @abstractmethod
+    def enforce_boundary(self, particle):
+        """Update the position of the particle if necessary to ensure that the
+        particle remains entirely within the domain.
+
+        Parameters
+        ----------
+        particle : numpy.ndarray
+            Cartesian coordinates of particle center.
+
+        """
+        pass
+
 
 class _CubicDomain(_Domain):
     """Cubic container in which to pack particles.
@@ -238,7 +248,7 @@ class _CubicDomain(_Domain):
         Length in x-, y-, and z- directions of each cell in mesh overlaid on
         domain.
     limits : list of float
-        Minimum and maximum position in x-, y-, and z-directions where particle
+        Maximum distance from center in x-, y-, or z-direction where particle
         center can be placed.
     volume : float
         Volume of the container.
@@ -256,9 +266,7 @@ class _CubicDomain(_Domain):
     @property
     def limits(self):
         if self._limits is None:
-            xlim = self.length/2 - self.particle_radius
-            self._limits = [[x - xlim for x in self.center],
-                            [x + xlim for x in self.center]]
+            self._limits = [self.length/2 - self.particle_radius]
         return self._limits
 
     @property
@@ -284,9 +292,14 @@ class _CubicDomain(_Domain):
         self._limits = limits
 
     def random_point(self):
-        return [uniform(self.limits[0][0], self.limits[1][0]),
-                uniform(self.limits[0][1], self.limits[1][1]),
-                uniform(self.limits[0][2], self.limits[1][2])]
+        x_max = self.limits[0]
+        return [uniform(-x_max, x_max),
+                uniform(-x_max, x_max),
+                uniform(-x_max, x_max)]
+
+    def enforce_boundary(self, particle):
+        x_max = self.limits[0]
+        particle[:] = np.clip(particle, -x_max, x_max)
 
 
 class _CylindricalDomain(_Domain):
@@ -317,8 +330,8 @@ class _CylindricalDomain(_Domain):
         Length in x-, y-, and z- directions of each cell in mesh overlaid on
         domain.
     limits : list of float
-        Minimum and maximum position in x-, y-, and z-directions where particle
-        center can be placed.
+        Maximum radial distance and maximum distance from center in z-direction
+        where particle center can be placed.
     volume : float
         Volume of the container.
 
@@ -340,12 +353,8 @@ class _CylindricalDomain(_Domain):
     @property
     def limits(self):
         if self._limits is None:
-            xlim = self.length/2 - self.particle_radius
-            rlim = self.radius - self.particle_radius
-            self._limits = [[self.center[0] - rlim, self.center[1] - rlim,
-                             self.center[2] - xlim],
-                            [self.center[0] + rlim, self.center[1] + rlim,
-                             self.center[2] + xlim]]
+            self._limits = [self.radius - self.particle_radius,
+                            self.length/2 - self.particle_radius]
         return self._limits
 
     @property
@@ -377,10 +386,20 @@ class _CylindricalDomain(_Domain):
         self._limits = limits
 
     def random_point(self):
-        r = sqrt(uniform(0, (self.radius - self.particle_radius)**2))
+        r_max = self.limits[0]
+        z_max = self.limits[1]
+        r = sqrt(uniform(0, r_max**2))
         t = uniform(0, 2*pi)
-        return [r*cos(t) + self.center[0], r*sin(t) + self.center[1],
-                uniform(self.limits[0][2], self.limits[1][2])]
+        return [r*cos(t), r*sin(t), uniform(-z_max, z_max)]
+
+    def enforce_boundary(self, particle):
+        r_max = self.limits[0]
+        z_max = self.limits[1]
+        d = sqrt(particle[0]**2 + particle[1]**2)
+        if d > r_max:
+            particle[0] *= r_max / d
+            particle[1] *= r_max / d
+        particle[2] = np.clip(particle[2], -z_max, z_max)
 
 
 class _SphericalDomain(_Domain):
@@ -407,8 +426,7 @@ class _SphericalDomain(_Domain):
         Length in x-, y-, and z- directions of each cell in mesh overlaid on
         domain.
     limits : list of float
-        Minimum and maximum position in x-, y-, and z-directions where particle
-        center can be placed.
+        Maximum radial distance where particle center can be placed.
     volume : float
         Volume of the container.
 
@@ -425,9 +443,7 @@ class _SphericalDomain(_Domain):
     @property
     def limits(self):
         if self._limits is None:
-            rlim = self.radius - self.particle_radius
-            self._limits = [[x - rlim for x in self.center],
-                            [x + rlim for x in self.center]]
+            self._limits = [self.radius - self.particle_radius]
         return self._limits
 
     @property
@@ -453,10 +469,18 @@ class _SphericalDomain(_Domain):
         self._limits = limits
 
     def random_point(self):
+        r_max = self.limits[0]
         x = (gauss(0, 1), gauss(0, 1), gauss(0, 1))
-        r = (uniform(0, (self.radius - self.particle_radius)**3)**(1/3) /
-             sqrt(x[0]**2 + x[1]**2 + x[2]**2))
-        return [r*x[i] + self.center[i] for i in range(3)]
+        r = (uniform(0, r_max**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2))
+        return [r*s for s in x]
+
+    def enforce_boundary(self, particle):
+        r_max = self.limits[0]
+        d = sqrt(particle[0]**2 + particle[1]**2 + particle[2]**2)
+        if d > r_max:
+            particle[0] *= r_max / d
+            particle[1] *= r_max / d
+            particle[2] *= r_max / d
 
 
 def create_triso_lattice(trisos, lower_left, pitch, shape, background):
@@ -767,9 +791,9 @@ def _close_random_pack(domain, particles, contraction_rate):
         particles[i] += r*v
         particles[j] -= r*v
 
-        # Apply reflective boundary conditions
-        particles[i] = particles[i].clip(domain.limits[0], domain.limits[1])
-        particles[j] = particles[j].clip(domain.limits[0], domain.limits[1])
+        # Apply boundary conditions
+        domain.enforce_boundary(particles[i])
+        domain.enforce_boundary(particles[j])
 
         update_mesh(i)
         update_mesh(j)
@@ -1008,8 +1032,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
     # Recalculate the limits for the initial random sequential packing using
     # the desired final particle radius to ensure particles 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]]]
+    domain.limits = [x + initial_radius - radius for x in domain.limits]
 
     # Generate non-overlapping particles for an initial inner radius using
     # random sequential packing algorithm
@@ -1024,5 +1047,5 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
 
     trisos = []
     for p in particles:
-        trisos.append(TRISO(radius, fill, p))
+        trisos.append(TRISO(radius, fill, [x + c for x, c in zip(p, domain.center)]))
     return trisos

From 00cb9149082e274dedacade9a0b171b467ab004c Mon Sep 17 00:00:00 2001
From: amandalund <alund@anl.gov>
Date: Mon, 26 Mar 2018 17:14:26 -0500
Subject: [PATCH 2/4] Add unit tests for TRISO

---
 openmc/model/triso.py                | 213 ++++++++++++++-------------
 tests/unit_tests/test_model_triso.py |  75 ++++++++++
 2 files changed, 187 insertions(+), 101 deletions(-)
 create mode 100644 tests/unit_tests/test_model_triso.py

diff --git a/openmc/model/triso.py b/openmc/model/triso.py
index b5ee74f90b..77488b891b 100644
--- a/openmc/model/triso.py
+++ b/openmc/model/triso.py
@@ -209,14 +209,21 @@ class _Domain(metaclass=ABCMeta):
         pass
 
     @abstractmethod
-    def enforce_boundary(self, particle):
-        """Update the position of the particle if necessary to ensure that the
-        particle remains entirely within the domain.
+    def repel_particles(self, p, q, d, d_new):
+        """Move particles 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
         ----------
-        particle : numpy.ndarray
+        p, q : numpy.ndarray
             Cartesian coordinates of particle center.
+        d : float
+            distance between centers of particles i and j.
+        d_new : float
+            final distance between centers of particles i and j.
 
         """
         pass
@@ -297,9 +304,20 @@ class _CubicDomain(_Domain):
                 uniform(-x_max, x_max),
                 uniform(-x_max, x_max)]
 
-    def enforce_boundary(self, particle):
+    def repel_particles(self, p, q, d, d_new):
+        # Moving each particle distance 'r' away from the other along the line
+        # joining the particle centers will ensure their final distance is
+        # equal to the outer diameter
+        r = (d_new - d)/2
+
+        v = (p - q)/d
+        p += r*v
+        q -= r*v
+
+        # Apply boundary conditions
         x_max = self.limits[0]
-        particle[:] = np.clip(particle, -x_max, x_max)
+        p[:] = np.clip(p, -x_max, x_max)
+        q[:] = np.clip(q, -x_max, x_max)
 
 
 class _CylindricalDomain(_Domain):
@@ -392,14 +410,29 @@ class _CylindricalDomain(_Domain):
         t = uniform(0, 2*pi)
         return [r*cos(t), r*sin(t), uniform(-z_max, z_max)]
 
-    def enforce_boundary(self, particle):
+    def repel_particles(self, p, q, d, d_new):
+        # Moving each particle distance 'r' away from the other along the line
+        # joining the particle centers will ensure their final distance is
+        # equal to the outer diameter
+        r = (d_new - d)/2
+
+        v = (p - q)/d
+        p += r*v
+        q -= r*v
+
+        # Apply boundary conditions
         r_max = self.limits[0]
         z_max = self.limits[1]
-        d = sqrt(particle[0]**2 + particle[1]**2)
+
+        d = sqrt(p[0]**2 + p[1]**2)
         if d > r_max:
-            particle[0] *= r_max / d
-            particle[1] *= r_max / d
-        particle[2] = np.clip(particle[2], -z_max, z_max)
+            p[0:2] *= r_max/d
+        p[2] = np.clip(p[2], -z_max, z_max)
+
+        d = sqrt(q[0]**2 + q[1]**2)
+        if d > r_max:
+            q[0:2] *= r_max/d
+        q[2] = np.clip(q[2], -z_max, z_max)
 
 
 class _SphericalDomain(_Domain):
@@ -474,13 +507,26 @@ class _SphericalDomain(_Domain):
         r = (uniform(0, r_max**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2))
         return [r*s for s in x]
 
-    def enforce_boundary(self, particle):
+    def repel_particles(self, p, q, d, d_new):
+        # Moving each particle distance 'r' away from the other along the line
+        # joining the particle centers will ensure their final distance is
+        # equal to the outer diameter
+        r = (d_new - d)/2
+
+        v = (p - q)/d
+        p += r*v
+        q -= r*v
+
+        # Apply boundary conditions
         r_max = self.limits[0]
-        d = sqrt(particle[0]**2 + particle[1]**2 + particle[2]**2)
+
+        d = sqrt(p[0]**2 + p[1]**2 + p[2]**2)
         if d > r_max:
-            particle[0] *= r_max / d
-            particle[1] *= r_max / d
-            particle[2] *= r_max / d
+            p *= r_max/d
+
+        d = sqrt(q[0]**2 + q[1]**2 + q[2]**2)
+        if d > r_max:
+            q *= r_max/d
 
 
 def create_triso_lattice(trisos, lower_left, pitch, shape, background):
@@ -711,13 +757,8 @@ def _close_random_pack(domain, particles, contraction_rate):
         for d, i, j in r:
             add_rod(d, i, j)
 
-        # Inner diameter is set initially to the shortest center-to-center
-        # distance between any two particles
-        if rods:
-            inner_diameter[0] = rods[0][0]
-
-    def update_mesh(i):
-        """Update which mesh cells the particle is in based on new particle
+    def update_mesh(indices):
+        """Update which mesh cells the particles are in based on new particle
         center coordinates.
 
         'mesh'/'mesh_map' is a two way dictionary used to look up which
@@ -727,22 +768,23 @@ def _close_random_pack(domain, particles, contraction_rate):
 
         Parameters
         ----------
-        i : int
-            Index of particle in particles array.
+        indices : List of int
+            Indices of particles in particles array.
 
         """
 
-        # Determine which mesh cells the particle is in and remove the
-        # particle id from those cells
-        for idx in mesh_map[i]:
-            mesh[idx].remove(i)
-        del mesh_map[i]
+        for i in indices:
+            # Determine which mesh cells the particle is in and remove the
+            # particle 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 particle's
-        # center and add this particle to the list of particles in those cells
-        for idx in domain.nearby_mesh_cells(particles[i]):
-            mesh[idx].add(i)
-            mesh_map[i].add(idx)
+            # Determine which mesh cells are within one diameter of particle's
+            # center and add this particle to the list of particles in those cells
+            for idx in domain.nearby_mesh_cells(particles[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:
@@ -753,50 +795,19 @@ def _close_random_pack(domain, particles, contraction_rate):
 
             j = floor(-log10(pf_out - pf_in)).
 
+        Returns
+        -------
+        float
+            New outer diameter
+
         """
 
-        inner_pf = (4/3 * pi * (inner_diameter[0]/2)**3 * n_particles /
-                    domain.volume)
-        outer_pf = (4/3 * pi * (outer_diameter[0]/2)**3 * n_particles /
-                    domain.volume)
+        inner_pf = 4/3*pi*(inner_diameter/2)**3*n_particles/domain.volume
+        outer_pf = 4/3*pi*(outer_diameter/2)**3*n_particles/domain.volume
 
         j = floor(-log10(outer_pf - inner_pf))
-        outer_diameter[0] = (outer_diameter[0] - 0.5**j * contraction_rate *
-                             initial_outer_diameter / n_particles)
-
-
-    def repel_particles(i, j, d):
-        """Move particles p and q apart according to the following
-        transformation (accounting for reflective 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
-        ----------
-        i, j : int
-            Index of particles in particles array.
-        d : float
-            distance between centers of particles i and j.
-
-        """
-
-        # Moving each particle distance 'r' away from the other along the line
-        # joining the particle centers will ensure their final distance is equal
-        # to the outer diameter
-        r = (outer_diameter[0] - d)/2
-
-        v = (particles[i] - particles[j])/d
-        particles[i] += r*v
-        particles[j] -= r*v
-
-        # Apply boundary conditions
-        domain.enforce_boundary(particles[i])
-        domain.enforce_boundary(particles[j])
-
-        update_mesh(i)
-        update_mesh(j)
+        return (outer_diameter - 0.5**j * contraction_rate *
+                initial_outer_diameter / n_particles)
 
     def nearest(i):
         """Find index of nearest neighbor of particle i.
@@ -827,33 +838,25 @@ def _close_random_pack(domain, particles, contraction_rate):
         else:
             return None, None
 
-    def update_rod_list(i, j):
-        """Update the rod list with the new nearest neighbors of particles i
-        and j since their overlap was eliminated.
+    def update_rod_list(indices):
+        """Update the rod list with the new nearest neighbors of particles in
+        list since their overlap was eliminated.
 
         Parameters
         ----------
-        i, j : int
-            Index of particles in particles array.
+        indices : List of int
+            Indices of particles in particles array.
 
         """
 
         # If the nearest neighbor k of particle 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:
-            remove_rod(k)
-            add_rod(d_ik, i, k)
-        l, d_jl = nearest(j)
-        if l and nearest(l)[0] == j:
-            remove_rod(l)
-            add_rod(d_jl, j, l)
-
-        # Set inner diameter to the shortest distance between two particle
-        # centers
-        if rods:
-            inner_diameter[0] = rods[0][0]
+        for i in indices:
+            k, d_ik = nearest(i)
+            if k and nearest(k)[0] == i:
+                remove_rod(k)
+                add_rod(d_ik, i, k)
 
     n_particles = len(particles)
     diameter = 2*domain.particle_radius
@@ -865,8 +868,8 @@ def _close_random_pack(domain, particles, contraction_rate):
     initial_outer_diameter = 2*(domain.volume/(n_particles*4/3*pi))**(1/3)
 
     # Inner and outer diameter of particles will change during packing
-    outer_diameter = [initial_outer_diameter]
-    inner_diameter = [0]
+    outer_diameter = initial_outer_diameter
+    inner_diameter = 0.
 
     rods = []
     rods_map = {}
@@ -880,14 +883,22 @@ def _close_random_pack(domain, particles, contraction_rate):
 
     while True:
         create_rod_list()
-        if inner_diameter[0] >= diameter:
+        if rods:
+            # Set inner diameter to the shortest center-to-center distance
+            # between any two particles
+            inner_diameter = rods[0][0]
+        if inner_diameter >= diameter:
             break
         while True:
             d, i, j = pop_rod()
-            reduce_outer_diameter()
-            repel_particles(i, j, d)
-            update_rod_list(i, j)
-            if inner_diameter[0] >= diameter or not rods:
+            outer_diameter = reduce_outer_diameter()
+            domain.repel_particles(particles[i], particles[j], d, outer_diameter)
+            update_mesh([i, j])
+            update_rod_list([i, j])
+            if not rods:
+                break
+            inner_diameter = rods[0][0]
+            if inner_diameter >= diameter:
                 break
 
 
@@ -957,7 +968,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
     to speed up the nearest neighbor search by only searching for a particle's
     neighbors within that mesh cell.
 
-    In CRP, each particle is assigned two diameters, and inner and an outer,
+    In CRP, each particle 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 particles and defines the pf. At each iteration the worst overlap
diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py
new file mode 100644
index 0000000000..c6f717d0b2
--- /dev/null
+++ b/tests/unit_tests/test_model_triso.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+
+from math import pi
+
+import numpy as np
+from numpy.linalg import norm
+import openmc
+import openmc.model
+import pytest
+import scipy.spatial
+
+
+_PACKING_FRACTION = 0.35
+_RADIUS = 1.
+_PARAMS = [{'shape' : 'cube', 'length' : 20., 'radius' : 0.},
+           {'shape' : 'cylinder', 'length' : 10., 'radius' : 10.},
+           {'shape' : 'sphere', 'length' : 0., 'radius' : 10.}]
+
+
+@pytest.fixture(scope='module', params=_PARAMS)
+def domain(request):
+    return request.param
+
+
+@pytest.fixture(scope='module')
+def trisos(domain):
+    return openmc.model.pack_trisos(
+        radius=_RADIUS,
+        fill=openmc.Universe(),
+        domain_shape=domain['shape'],
+        domain_length=domain['length'],
+        domain_radius=domain['radius'],
+        domain_center=[0., 0., 0.],
+        initial_packing_fraction=0.2,
+        packing_fraction=_PACKING_FRACTION
+    )
+
+
+def test_overlap(trisos):
+    """Check that no TRISO particles overlap."""
+    tree = scipy.spatial.cKDTree([t.center for t in trisos])
+    r = min(tree.query([t.center for t in trisos], k=2)[0][:,1])/2
+    assert r > _RADIUS or r == pytest.approx(_RADIUS)
+
+
+def test_contained(trisos, domain):
+    """Make sure all particles are entirely contained within the domain."""
+    if domain['shape'] == 'cube':
+        x = 2*(max(np.hstack([abs(t.center) for t in trisos])) + _RADIUS)
+        assert x < domain['length'] or x == pytest.approx(domain['length'])
+
+    elif domain['shape'] == 'cylinder':
+        r = max([norm(t.center[0:2]) for t in trisos]) + _RADIUS
+        z = 2*(max([abs(t.center[2]) for t in trisos]) + _RADIUS)
+        assert r < domain['radius'] or r == pytest.approx(domain['radius'])
+        assert z < domain['length'] or z == pytest.approx(domain['length'])
+
+    elif domain['shape'] == 'sphere':
+        r = max([norm(t.center) for t in trisos]) + _RADIUS
+        assert r < domain['radius'] or r == pytest.approx(domain['radius'])
+
+
+def test_packing_fraction(trisos, domain):
+    """Check that the actual PF is close to the requested PF."""
+    if domain['shape'] == 'cube':
+        volume = domain['length']**3
+
+    elif domain['shape'] == 'cylinder':
+        volume = domain['length']*pi*domain['radius']**2
+
+    elif domain['shape'] == 'sphere':
+        volume = 4/3*pi*domain['radius']**3
+
+    pf = len(trisos)*4/3*pi*_RADIUS**3/volume
+    assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2)

From 8838f14720f23c198bc112a2a0fc957917b3379a Mon Sep 17 00:00:00 2001
From: amandalund <alund@anl.gov>
Date: Fri, 30 Mar 2018 20:00:28 -0500
Subject: [PATCH 3/4] Added some more TRISO unit tests

---
 openmc/model/triso.py                |  12 ++-
 tests/unit_tests/test_model_triso.py | 144 ++++++++++++++++++++++-----
 2 files changed, 127 insertions(+), 29 deletions(-)

diff --git a/openmc/model/triso.py b/openmc/model/triso.py
index 77488b891b..7040d02cab 100644
--- a/openmc/model/triso.py
+++ b/openmc/model/triso.py
@@ -314,7 +314,9 @@ class _CubicDomain(_Domain):
         p += r*v
         q -= r*v
 
-        # Apply boundary conditions
+        # Enforce the rigid boundary by moving each particle back along the
+        # surface normal until it is completely within the container if it
+        # overlaps the surface
         x_max = self.limits[0]
         p[:] = np.clip(p, -x_max, x_max)
         q[:] = np.clip(q, -x_max, x_max)
@@ -420,7 +422,9 @@ class _CylindricalDomain(_Domain):
         p += r*v
         q -= r*v
 
-        # Apply boundary conditions
+        # Enforce the rigid boundary by moving each particle back along the
+        # surface normal until it is completely within the container if it
+        # overlaps the surface
         r_max = self.limits[0]
         z_max = self.limits[1]
 
@@ -517,7 +521,9 @@ class _SphericalDomain(_Domain):
         p += r*v
         q -= r*v
 
-        # Apply boundary conditions
+        # Enforce the rigid boundary by moving each particle back along the
+        # surface normal until it is completely within the container if it
+        # overlaps the surface
         r_max = self.limits[0]
 
         d = sqrt(p[0]**2 + p[1]**2 + p[2]**2)
diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py
index c6f717d0b2..bfd7319d21 100644
--- a/tests/unit_tests/test_model_triso.py
+++ b/tests/unit_tests/test_model_triso.py
@@ -11,49 +11,75 @@ import scipy.spatial
 
 
 _PACKING_FRACTION = 0.35
-_RADIUS = 1.
-_PARAMS = [{'shape' : 'cube', 'length' : 20., 'radius' : 0.},
-           {'shape' : 'cylinder', 'length' : 10., 'radius' : 10.},
-           {'shape' : 'sphere', 'length' : 0., 'radius' : 10.}]
+_RADIUS = 4.25e-2
 
 
-@pytest.fixture(scope='module', params=_PARAMS)
+@pytest.fixture(scope='module',
+                params=[{'shape': 'cube',
+                         'length': 0.75,
+                         'radius': 0.,
+                         'volume': 0.75**3},
+                        {'shape': 'cylinder',
+                         'length': 0.5,
+                         'radius': 0.5,
+                         'volume': 0.5*pi*0.5**2},
+                        {'shape': 'sphere',
+                         'length': 0.,
+                         'radius': 0.5,
+                         'volume': 4/3*pi*0.5**3}])
 def domain(request):
     return request.param
 
 
 @pytest.fixture(scope='module')
-def trisos(domain):
-    return openmc.model.pack_trisos(
+def triso_universe():
+    sphere = openmc.Sphere(R=_RADIUS)
+    cell = openmc.Cell(region=-sphere)
+    univ = openmc.Universe(cells=[cell])
+    return univ
+
+
+@pytest.fixture(scope='module')
+def trisos(domain, triso_universe):
+    trisos = openmc.model.pack_trisos(
         radius=_RADIUS,
-        fill=openmc.Universe(),
+        fill=triso_universe,
         domain_shape=domain['shape'],
         domain_length=domain['length'],
         domain_radius=domain['radius'],
-        domain_center=[0., 0., 0.],
+        domain_center=(0., 0., 0.),
         initial_packing_fraction=0.2,
         packing_fraction=_PACKING_FRACTION
     )
+    return trisos
 
 
 def test_overlap(trisos):
     """Check that no TRISO particles overlap."""
-    tree = scipy.spatial.cKDTree([t.center for t in trisos])
-    r = min(tree.query([t.center for t in trisos], k=2)[0][:,1])/2
-    assert r > _RADIUS or r == pytest.approx(_RADIUS)
+    centers = [t.center for t in trisos]
+
+    # Create KD tree for quick nearest neighbor search
+    tree = scipy.spatial.cKDTree(centers)
+
+    # Find distance to nearest neighbor for all particles
+    d = tree.query(centers, k=2)[0]
+
+    # Get the smallest distance between any two particles
+    d_min = min(d[:,1])
+    assert d_min > 2*_RADIUS or d_min == pytest.approx(2*_RADIUS)
 
 
 def test_contained(trisos, domain):
     """Make sure all particles are entirely contained within the domain."""
     if domain['shape'] == 'cube':
-        x = 2*(max(np.hstack([abs(t.center) for t in trisos])) + _RADIUS)
-        assert x < domain['length'] or x == pytest.approx(domain['length'])
+        x = max(np.hstack([abs(t.center) for t in trisos])) + _RADIUS
+        assert x < 0.5*domain['length'] or x == pytest.approx(0.5*domain['length'])
 
     elif domain['shape'] == 'cylinder':
         r = max([norm(t.center[0:2]) for t in trisos]) + _RADIUS
-        z = 2*(max([abs(t.center[2]) for t in trisos]) + _RADIUS)
+        z = max([abs(t.center[2]) for t in trisos]) + _RADIUS
         assert r < domain['radius'] or r == pytest.approx(domain['radius'])
-        assert z < domain['length'] or z == pytest.approx(domain['length'])
+        assert z < 0.5*domain['length'] or z == pytest.approx(0.5*domain['length'])
 
     elif domain['shape'] == 'sphere':
         r = max([norm(t.center) for t in trisos]) + _RADIUS
@@ -62,14 +88,80 @@ def test_contained(trisos, domain):
 
 def test_packing_fraction(trisos, domain):
     """Check that the actual PF is close to the requested PF."""
-    if domain['shape'] == 'cube':
-        volume = domain['length']**3
-
-    elif domain['shape'] == 'cylinder':
-        volume = domain['length']*pi*domain['radius']**2
-
-    elif domain['shape'] == 'sphere':
-        volume = 4/3*pi*domain['radius']**3
-
-    pf = len(trisos)*4/3*pi*_RADIUS**3/volume
+    pf = len(trisos)*4/3*pi*_RADIUS**3/domain['volume']
     assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2)
+
+
+def test_n_particles(triso_universe):
+    """Check that the function returns the correct number of particles"""
+    trisos = openmc.model.pack_trisos(
+        radius=_RADIUS, fill=triso_universe, domain_shape='cube',
+        domain_length=1.0, n_particles=800
+    )
+    assert len(trisos) == 800
+    
+
+def test_triso_lattice(triso_universe):
+    trisos = openmc.model.pack_trisos(
+        radius=_RADIUS, fill=triso_universe, domain_shape='cube',
+        domain_length=1.0, domain_center=(0., 0., 0.), packing_fraction=0.2
+    )
+
+    lower_left = np.array((-.5, -.5, -.5))
+    upper_right = np.array((.5, .5, .5))
+    shape = (3, 3, 3)
+    pitch = (upper_right - lower_left)/shape
+    background = openmc.Material()
+
+    lattice = openmc.model.create_triso_lattice(
+        trisos, lower_left, pitch, shape, background
+    )
+
+
+def test_domain_input(triso_universe):
+    # Invalid domain shape
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, n_particles=100,
+            domain_shape='circle'
+        )
+    # Don't specify domain length on a cube
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, n_particles=100,
+            domain_shape='cube'
+        )
+    # Don't specify domain radius on a sphere
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, n_particles=100,
+            domain_shape='sphere'
+        )
+    
+
+def test_packing_fraction_input(triso_universe):
+    # Provide neither packing fraction nor number of particles
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, domain_shape='cube',
+            domain_length=10
+        )
+    # Provide both packing fraction and number of particles
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, domain_shape='cube',
+            domain_length=10, n_particles=100, packing_fraction=0.2
+        )
+    # Specify a packing fraction that is too high for CRP
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, domain_shape='cube',
+            domain_length=10, packing_fraction=1
+        )
+    # Specify a packing fraction that is too high for RSP
+    with pytest.raises(ValueError):
+        trisos = openmc.model.pack_trisos(
+            radius=1, fill=triso_universe, domain_shape='cube',
+            domain_length=10, packing_fraction=0.5,
+            initial_packing_fraction=0.4
+        )

From 4a64e2a93139574a0d913b9b4c8266264a15bea2 Mon Sep 17 00:00:00 2001
From: amandalund <alund@anl.gov>
Date: Sat, 14 Apr 2018 14:10:35 -0500
Subject: [PATCH 4/4] Address #991 comments; add comments; add stopping
 condition and warning on early convergence

---
 openmc/model/triso.py                | 150 ++++++++++++++++-----------
 tests/unit_tests/test_model_triso.py |  22 ++--
 2 files changed, 97 insertions(+), 75 deletions(-)

diff --git a/openmc/model/triso.py b/openmc/model/triso.py
index 7040d02cab..1e4b4aa9b1 100644
--- a/openmc/model/triso.py
+++ b/openmc/model/triso.py
@@ -305,14 +305,14 @@ class _CubicDomain(_Domain):
                 uniform(-x_max, x_max)]
 
     def repel_particles(self, p, q, d, d_new):
-        # Moving each particle distance 'r' away from the other along the line
+        # Moving each particle distance 's' away from the other along the line
         # joining the particle centers will ensure their final distance is
         # equal to the outer diameter
-        r = (d_new - d)/2
+        s = (d_new - d)/2
 
         v = (p - q)/d
-        p += r*v
-        q -= r*v
+        p += s*v
+        q -= s*v
 
         # Enforce the rigid boundary by moving each particle back along the
         # surface normal until it is completely within the container if it
@@ -413,14 +413,14 @@ class _CylindricalDomain(_Domain):
         return [r*cos(t), r*sin(t), uniform(-z_max, z_max)]
 
     def repel_particles(self, p, q, d, d_new):
-        # Moving each particle distance 'r' away from the other along the line
+        # Moving each particle distance 's' away from the other along the line
         # joining the particle centers will ensure their final distance is
         # equal to the outer diameter
-        r = (d_new - d)/2
+        s = (d_new - d)/2
 
         v = (p - q)/d
-        p += r*v
-        q -= r*v
+        p += s*v
+        q -= s*v
 
         # Enforce the rigid boundary by moving each particle back along the
         # surface normal until it is completely within the container if it
@@ -428,14 +428,14 @@ class _CylindricalDomain(_Domain):
         r_max = self.limits[0]
         z_max = self.limits[1]
 
-        d = sqrt(p[0]**2 + p[1]**2)
-        if d > r_max:
-            p[0:2] *= r_max/d
+        r = sqrt(p[0]**2 + p[1]**2)
+        if r > r_max:
+            p[0:2] *= r_max/r
         p[2] = np.clip(p[2], -z_max, z_max)
 
-        d = sqrt(q[0]**2 + q[1]**2)
-        if d > r_max:
-            q[0:2] *= r_max/d
+        r = sqrt(q[0]**2 + q[1]**2)
+        if r > r_max:
+            q[0:2] *= r_max/r
         q[2] = np.clip(q[2], -z_max, z_max)
 
 
@@ -512,27 +512,27 @@ class _SphericalDomain(_Domain):
         return [r*s for s in x]
 
     def repel_particles(self, p, q, d, d_new):
-        # Moving each particle distance 'r' away from the other along the line
+        # Moving each particle distance 's' away from the other along the line
         # joining the particle centers will ensure their final distance is
         # equal to the outer diameter
-        r = (d_new - d)/2
+        s = (d_new - d)/2
 
         v = (p - q)/d
-        p += r*v
-        q -= r*v
+        p += s*v
+        q -= s*v
 
         # Enforce the rigid boundary by moving each particle back along the
         # surface normal until it is completely within the container if it
         # overlaps the surface
         r_max = self.limits[0]
 
-        d = sqrt(p[0]**2 + p[1]**2 + p[2]**2)
-        if d > r_max:
-            p *= r_max/d
+        r = sqrt(p[0]**2 + p[1]**2 + p[2]**2)
+        if r > r_max:
+            p *= r_max/r
 
-        d = sqrt(q[0]**2 + q[1]**2 + q[2]**2)
-        if d > r_max:
-            q *= r_max/d
+        r = sqrt(q[0]**2 + q[1]**2 + q[2]**2)
+        if r > r_max:
+            q *= r_max/r
 
 
 def create_triso_lattice(trisos, lower_left, pitch, shape, background):
@@ -712,6 +712,7 @@ def _close_random_pack(domain, particles, contraction_rate):
                 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 particle centers).
@@ -736,8 +737,8 @@ def _close_random_pack(domain, particles, contraction_rate):
         # Find distance to nearest neighbor and index of nearest neighbor for
         # all particles
         d, n = tree.query(particles, k=2)
-        d = d[:,1]
-        n = n[:,1]
+        d = d[:, 1]
+        n = n[:, 1]
 
         # Array of particle indices, indices of nearest neighbors, and
         # distances to nearest neighbors
@@ -746,8 +747,8 @@ def _close_random_pack(domain, particles, contraction_rate):
         # Sort along second column and swap first and second columns to create
         # array of nearest neighbor indices, indices of particles they are
         # nearest neighbors of, and distances between them
-        b = a[a[:,1].argsort()]
-        b[:,[0, 1]] = b[:,[1, 0]]
+        b = a[a[:, 1].argsort()]
+        b[:, [0, 1]] = b[:, [1, 0]]
 
         # Find the intersection between 'a' and 'b': a list of particles who
         # are each other's nearest neighbors and the distance between them
@@ -761,10 +762,11 @@ def _close_random_pack(domain, particles, contraction_rate):
         del rods[:]
         rods_map.clear()
         for d, i, j in r:
-            add_rod(d, i, j)
+            if d < outer_diameter and not np.isclose(d, outer_diameter, atol=1.0e-14):
+                add_rod(d, i, j)
 
-    def update_mesh(indices):
-        """Update which mesh cells the particles are in based on new particle
+    def update_mesh(i):
+        """Update which mesh cells the particle is in based on new particle
         center coordinates.
 
         'mesh'/'mesh_map' is a two way dictionary used to look up which
@@ -774,23 +776,22 @@ def _close_random_pack(domain, particles, contraction_rate):
 
         Parameters
         ----------
-        indices : List of int
-            Indices of particles in particles array.
+        i : int
+            Index of particle in particles array.
 
         """
 
-        for i in indices:
-            # Determine which mesh cells the particle is in and remove the
-            # particle id from those cells
-            for idx in mesh_map[i]:
-                mesh[idx].remove(i)
-            del mesh_map[i]
+        # Determine which mesh cells the particle is in and remove the
+        # particle 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 particle's
-            # center and add this particle to the list of particles in those cells
-            for idx in domain.nearby_mesh_cells(particles[i]):
-                mesh[idx].add(i)
-                mesh_map[i].add(idx)
+        # Determine which mesh cells are within one diameter of particle's
+        # center and add this particle to the list of particles in those cells
+        for idx in domain.nearby_mesh_cells(particles[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:
@@ -844,25 +845,25 @@ def _close_random_pack(domain, particles, contraction_rate):
         else:
             return None, None
 
-    def update_rod_list(indices):
-        """Update the rod list with the new nearest neighbors of particles in
-        list since their overlap was eliminated.
+    def update_rod_list(i):
+        """Update the rod list with the new nearest neighbors of particle since
+        its overlap was eliminated.
 
         Parameters
         ----------
-        indices : List of int
-            Indices of particles in particles array.
+        i : int
+            Index of particle in particles array.
 
         """
 
         # If the nearest neighbor k of particle 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
-        for i in indices:
-            k, d_ik = nearest(i)
-            if k and nearest(k)[0] == i:
-                remove_rod(k)
-                add_rod(d_ik, i, k)
+        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)
 
     n_particles = len(particles)
     diameter = 2*domain.particle_radius
@@ -877,41 +878,68 @@ def _close_random_pack(domain, particles, contraction_rate):
     outer_diameter = initial_outer_diameter
     inner_diameter = 0.
 
+    # List of rods arranged in a heap and mapping of particle ids to rods
     rods = []
     rods_map = {}
+
+    # Initialize two-way dictionary that identifies which particles are near a
+    # given mesh cell and which mesh cells a particle is near
     mesh = defaultdict(set)
     mesh_map = defaultdict(set)
-
     for i in range(n_particles):
         for idx in domain.nearby_mesh_cells(particles[i]):
             mesh[idx].add(i)
             mesh_map[i].add(idx)
 
     while True:
+        # Rebuild the sorted list of rods according to the current particle
+        # configuration
         create_rod_list()
+
+        # Set the inner diameter to the shortest center-to-center distance
+        # between any two particles
         if rods:
-            # Set inner diameter to the shortest center-to-center distance
-            # between any two particles
             inner_diameter = rods[0][0]
+
+        # Reached the desired particle radius
         if inner_diameter >= diameter:
             break
+
+        # The algorithm converged before reaching the desired particle 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 particle 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 '
+                          'particle radius; some particles 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_particles(particles[i], particles[j], d, outer_diameter)
-            update_mesh([i, j])
-            update_rod_list([i, j])
+            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:
+            if inner_diameter >= diameter or inner_diameter >= outer_diameter:
                 break
 
 
 def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
                 domain_radius=None, domain_center=[0., 0., 0.],
                 n_particles=None, packing_fraction=None,
-                initial_packing_fraction=0.3, contraction_rate=1/400, seed=1):
+                initial_packing_fraction=0.3, contraction_rate=1.e-3, seed=1):
     """Generate a random, non-overlapping configuration of TRISO particles
     within a container.
 
diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py
index bfd7319d21..19c4f9081e 100644
--- a/tests/unit_tests/test_model_triso.py
+++ b/tests/unit_tests/test_model_triso.py
@@ -12,21 +12,15 @@ import scipy.spatial
 
 _PACKING_FRACTION = 0.35
 _RADIUS = 4.25e-2
+domain_params = [
+    {'shape': 'cube', 'length': 0.75, 'radius': 0., 'volume': 0.75**3},
+    {'shape': 'cylinder', 'length': 0.5, 'radius': 0.5, 'volume': 0.5*pi*0.5**2},
+    {'shape': 'sphere', 'length': 0., 'radius': 0.5, 'volume': 4/3*pi*0.5**3}
+]
 
 
-@pytest.fixture(scope='module',
-                params=[{'shape': 'cube',
-                         'length': 0.75,
-                         'radius': 0.,
-                         'volume': 0.75**3},
-                        {'shape': 'cylinder',
-                         'length': 0.5,
-                         'radius': 0.5,
-                         'volume': 0.5*pi*0.5**2},
-                        {'shape': 'sphere',
-                         'length': 0.,
-                         'radius': 0.5,
-                         'volume': 4/3*pi*0.5**3}])
+@pytest.fixture(scope='module', params=domain_params,
+                ids=['cube', 'cylinder', 'sphere'])
 def domain(request):
     return request.param
 
@@ -65,7 +59,7 @@ def test_overlap(trisos):
     d = tree.query(centers, k=2)[0]
 
     # Get the smallest distance between any two particles
-    d_min = min(d[:,1])
+    d_min = min(d[:, 1])
     assert d_min > 2*_RADIUS or d_min == pytest.approx(2*_RADIUS)