1#!/usr/bin/env python
2
3from math import pi
4
5import numpy as np
6from numpy.linalg import norm
7import openmc
8import openmc.model
9import pytest
10import scipy.spatial
11
12
13_RADIUS = 0.1
14_PACKING_FRACTION = 0.35
15_PARAMS = [
16    {'shape': 'rectangular_prism', 'volume': 1**3},
17    {'shape': 'x_cylinder', 'volume': 1*pi*1**2},
18    {'shape': 'y_cylinder', 'volume': 1*pi*1**2},
19    {'shape': 'z_cylinder', 'volume': 1*pi*1**2},
20    {'shape': 'sphere', 'volume': 4/3*pi*1**3},
21    {'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)}
22]
23
24
25@pytest.fixture(scope='module', params=_PARAMS)
26def container(request):
27    return request.param
28
29
30@pytest.fixture(scope='module')
31def centers(request, container):
32    return request.getfixturevalue('centers_' + container['shape'])
33
34
35@pytest.fixture(scope='module')
36def centers_rectangular_prism():
37    min_x = openmc.XPlane(0)
38    max_x = openmc.XPlane(1)
39    min_y = openmc.YPlane(0)
40    max_y = openmc.YPlane(1)
41    min_z = openmc.ZPlane(0)
42    max_z = openmc.ZPlane(1)
43    region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z
44    return openmc.model.pack_spheres(radius=_RADIUS, region=region,
45        pf=_PACKING_FRACTION, initial_pf=0.2)
46
47
48@pytest.fixture(scope='module')
49def centers_x_cylinder():
50    cylinder = openmc.XCylinder(r=1, y0=1, z0=2)
51    min_x = openmc.XPlane(0)
52    max_x = openmc.XPlane(1)
53    region = +min_x & -max_x & -cylinder
54    return openmc.model.pack_spheres(radius=_RADIUS, region=region,
55        pf=_PACKING_FRACTION, initial_pf=0.2)
56
57
58@pytest.fixture(scope='module')
59def centers_y_cylinder():
60    cylinder = openmc.YCylinder(r=1, x0=1, z0=2)
61    min_y = openmc.YPlane(0)
62    max_y = openmc.YPlane(1)
63    region = +min_y & -max_y & -cylinder
64    return openmc.model.pack_spheres(radius=_RADIUS, region=region,
65        pf=_PACKING_FRACTION, initial_pf=0.2)
66
67
68@pytest.fixture(scope='module')
69def centers_z_cylinder():
70    cylinder = openmc.ZCylinder(r=1, x0=1, y0=2)
71    min_z = openmc.ZPlane(0)
72    max_z = openmc.ZPlane(1)
73    region = +min_z & -max_z & -cylinder
74    return openmc.model.pack_spheres(radius=_RADIUS, region=region,
75        pf=_PACKING_FRACTION, initial_pf=0.2)
76
77
78@pytest.fixture(scope='module')
79def centers_sphere():
80    sphere = openmc.Sphere(r=1, x0=1, y0=2, z0=3)
81    region = -sphere
82    return openmc.model.pack_spheres(radius=_RADIUS, region=region,
83        pf=_PACKING_FRACTION, initial_pf=0.2)
84
85
86@pytest.fixture(scope='module')
87def centers_spherical_shell():
88    sphere = openmc.Sphere(r=1, x0=1, y0=2, z0=3)
89    inner_sphere = openmc.Sphere(r=0.5, x0=1, y0=2, z0=3)
90    region = -sphere & +inner_sphere
91    return openmc.model.pack_spheres(radius=_RADIUS, region=region,
92        pf=_PACKING_FRACTION, initial_pf=0.2)
93
94
95@pytest.fixture(scope='module')
96def triso_universe():
97    sphere = openmc.Sphere(r=_RADIUS)
98    cell = openmc.Cell(region=-sphere)
99    univ = openmc.Universe(cells=[cell])
100    return univ
101
102
103def test_overlap(centers):
104    """Check that none of the spheres in the packed configuration overlap."""
105    # Create KD tree for quick nearest neighbor search
106    tree = scipy.spatial.cKDTree(centers)
107
108    # Find distance to nearest neighbor for all spheres
109    d = tree.query(centers, k=2)[0]
110
111    # Get the smallest distance between any two spheres
112    d_min = min(d[:, 1])
113    assert d_min > 2*_RADIUS or d_min == pytest.approx(2*_RADIUS)
114
115
116def test_contained_rectangular_prism(centers_rectangular_prism):
117    """Make sure all spheres are entirely contained within the domain."""
118    d_max = np.amax(centers_rectangular_prism) + _RADIUS
119    d_min = np.amin(centers_rectangular_prism) - _RADIUS
120    assert d_max < 1 or d_max == pytest.approx(1)
121    assert d_min > 0 or d_min == pytest.approx(0)
122
123
124def test_contained_x_cylinder(centers_x_cylinder):
125    """Make sure all spheres are entirely contained within the domain."""
126    d = np.linalg.norm(centers_x_cylinder[:,[1,2]] - [1, 2], axis=1)
127    r_max = max(d) + _RADIUS
128    x_max = max(centers_x_cylinder[:,0]) + _RADIUS
129    x_min = min(centers_x_cylinder[:,0]) - _RADIUS
130    assert r_max < 1 or r_max == pytest.approx(1)
131    assert x_max < 1 or x_max == pytest.approx(1)
132    assert x_min > 0 or x_min == pytest.approx(0)
133
134
135def test_contained_y_cylinder(centers_y_cylinder):
136    """Make sure all spheres are entirely contained within the domain."""
137    d = np.linalg.norm(centers_y_cylinder[:,[0,2]] - [1, 2], axis=1)
138    r_max = max(d) + _RADIUS
139    y_max = max(centers_y_cylinder[:,1]) + _RADIUS
140    y_min = min(centers_y_cylinder[:,1]) - _RADIUS
141    assert r_max < 1 or r_max == pytest.approx(1)
142    assert y_max < 1 or y_max == pytest.approx(1)
143    assert y_min > 0 or y_min == pytest.approx(0)
144
145
146def test_contained_z_cylinder(centers_z_cylinder):
147    """Make sure all spheres are entirely contained within the domain."""
148    d = np.linalg.norm(centers_z_cylinder[:,[0,1]] - [1, 2], axis=1)
149    r_max = max(d) + _RADIUS
150    z_max = max(centers_z_cylinder[:,2]) + _RADIUS
151    z_min = min(centers_z_cylinder[:,2]) - _RADIUS
152    assert r_max < 1 or r_max == pytest.approx(1)
153    assert z_max < 1 or z_max == pytest.approx(1)
154    assert z_min > 0 or z_min == pytest.approx(0)
155
156
157def test_contained_sphere(centers_sphere):
158    """Make sure all spheres are entirely contained within the domain."""
159    d = np.linalg.norm(centers_sphere - [1, 2, 3], axis=1)
160    r_max = max(d) + _RADIUS
161    assert r_max < 1 or r_max == pytest.approx(1)
162
163
164def test_contained_spherical_shell(centers_spherical_shell):
165    """Make sure all spheres are entirely contained within the domain."""
166    d = np.linalg.norm(centers_spherical_shell - [1, 2, 3], axis=1)
167    r_max = max(d) + _RADIUS
168    r_min = min(d) - _RADIUS
169    assert r_max < 1 or r_max == pytest.approx(1)
170    assert r_min > 0.5 or r_min == pytest.approx(0.5)
171
172
173def test_packing_fraction(container, centers):
174    """Check that the actual PF is close to the requested PF."""
175    pf = len(centers) * 4/3 * pi *_RADIUS**3 / container['volume']
176    assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2)
177
178
179def test_num_spheres():
180    """Check that the function returns the correct number of spheres"""
181    centers = openmc.model.pack_spheres(
182        radius=_RADIUS, region=-openmc.Sphere(r=1), num_spheres=50
183    )
184    assert len(centers) == 50
185
186
187def test_triso_lattice(triso_universe, centers_rectangular_prism):
188    trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c)
189              for c in centers_rectangular_prism]
190
191    lower_left = np.array((0, 0, 0))
192    upper_right = np.array((1, 1, 1))
193    shape = (3, 3, 3)
194    pitch = (upper_right - lower_left)/shape
195    background = openmc.Material()
196
197    lattice = openmc.model.create_triso_lattice(
198        trisos, lower_left, pitch, shape, background
199    )
200
201
202def test_container_input(triso_universe):
203    # Invalid container shape
204    with pytest.raises(ValueError):
205        centers = openmc.model.pack_spheres(
206            radius=_RADIUS, region=+openmc.Sphere(r=1), num_spheres=100
207        )
208
209
210def test_packing_fraction_input():
211    # Provide neither packing fraction nor number of spheres
212    with pytest.raises(ValueError):
213        centers = openmc.model.pack_spheres(
214            radius=_RADIUS, region=-openmc.Sphere(r=1)
215        )
216
217    # Specify a packing fraction that is too high for CRP
218    with pytest.raises(ValueError):
219        centers = openmc.model.pack_spheres(
220            radius=_RADIUS, region=-openmc.Sphere(r=1), pf=1
221        )
222
223    # Specify a packing fraction that is too high for RSP
224    with pytest.raises(ValueError):
225        centers = openmc.model.pack_spheres(
226            radius=_RADIUS, region=-openmc.Sphere(r=1), pf=0.5, initial_pf=0.4
227        )
228