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