diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 03ca6f5e03..a99f136bb3 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -882,21 +882,21 @@ def _create_container(region, sphere_radius): if p1.surface.x0 > p2.surface.x0: p1, p2 = p2, p1 length = p2.surface.x0 - p1.surface.x0 - center = (length/2, cyl.surface.y0, cyl.surface.z0) + center = (p1.surface.x0 + length/2, cyl.surface.y0, cyl.surface.z0) # Calculate the parameters for a cylinder 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 = (cyl.surface.x0, length/2, cyl.surface.z0) + center = (cyl.surface.x0, p1.surface.y0 + length/2, cyl.surface.z0) # Calculate the parameters for a cylinder along the z-axis else: if p1.surface.z0 > p2.surface.z0: p1, p2 = p2, p1 length = p2.surface.z0 - p1.surface.z0 - center = (cyl.surface.x0, cyl.surface.y0, length/2) + center = (cyl.surface.x0, cyl.surface.y0, p1.surface.z0 + length/2) # Make sure the half-spaces are on the correct side of the surfaces if cyl.side != '-' or p1.side != '+' or p2.side != '-': @@ -1424,23 +1424,3 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, _close_random_pack(domain, spheres, contraction_rate) return spheres + domain.center - -def create_trisos(outer_radius, fill, centers): - """Create TRISO particles at the given coordinates. - - Parameters - ---------- - outer_radius : float - Outer radius of TRISO particle - fill : openmc.Universe - Universe which contains all layers of the TRISO particle - centers : Iterable of float - Cartesian coordinates of the centers of the TRISO particles in cm - - Returns - ------- - list of openmc.model.TRISO - List of TRISO particles in the domain. - - """ - return [TRISO(outer_radius, fill, c) for c in centers] diff --git a/tests/regression_tests/triso/inputs_true.dat b/tests/regression_tests/triso/inputs_true.dat index a131876285..7e5007eca2 100644 --- a/tests/regression_tests/triso/inputs_true.dat +++ b/tests/regression_tests/triso/inputs_true.dat @@ -5,7 +5,7 @@ - + @@ -220,12 +220,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index dfc700c058..2bbecb4a92 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.683226E+00 7.383559E-02 +1.683227E+00 7.383566E-02 diff --git a/tests/regression_tests/triso/test.py b/tests/regression_tests/triso/test.py index 174bba94fd..9ede6dbb03 100644 --- a/tests/regression_tests/triso/test.py +++ b/tests/regression_tests/triso/test.py @@ -53,19 +53,21 @@ class TRISOTestHarness(PyAPITestHarness): c5 = openmc.Cell(fill=opyc, region=+spheres[3]) inner_univ = openmc.Universe(cells=[c1, c2, c3, c4, c5]) - outer_radius = 422.5*1e-4 - trisos = openmc.model.pack_trisos( - radius=outer_radius, fill=inner_univ, domain_shape='cube', - domain_length=1., domain_center=(0., 0., 0.), n_particles=100) - - # Define box to contain lattice + # Define box to contain lattice and to pack TRISO particles in min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective') max_x = openmc.XPlane(x0=0.5, boundary_type='reflective') min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective') max_y = openmc.YPlane(y0=0.5, boundary_type='reflective') min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective') max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective') - box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z) + box_region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + box = openmc.Cell(region=box_region) + + outer_radius = 422.5*1e-4 + centers = openmc.model.pack_spheres(radius=outer_radius, + region=box_region, num_spheres=100) + trisos = [openmc.model.TRISO(outer_radius, inner_univ, c) + for c in centers] # Create lattice ll, ur = box.region.bounding_box diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index 19c4f9081e..3c261ed6f6 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -10,19 +10,54 @@ import pytest import scipy.spatial +_RADIUS = 0.1 _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} -] +_SHAPES = ['rectangular_prism', 'cylinder', 'sphere', 'spherical_shell'] +_VOLUMES = [1.**3, 1.*pi*1.**2, 4/3*pi*1.**3, 4/3*pi*(1.**3 - 0.5**3)] -@pytest.fixture(scope='module', params=domain_params, - ids=['cube', 'cylinder', 'sphere']) -def domain(request): - return request.param +@pytest.fixture(scope='module') +def container(request): + return request.getfixturevalue(request.param) + + +@pytest.fixture(scope='module') +def centers(request): + container = request.getfixturevalue(request.param) + return openmc.model.pack_spheres(radius=_RADIUS, region=container, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def rectangular_prism(): + min_x = openmc.XPlane(x0=-0.5) + max_x = openmc.XPlane(x0=0.5) + min_y = openmc.YPlane(y0=-0.5) + max_y = openmc.YPlane(y0=0.5) + min_z = openmc.ZPlane(z0=-0.5) + max_z = openmc.ZPlane(z0=0.5) + return +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + + +@pytest.fixture(scope='module') +def cylinder(): + cylinder = openmc.ZCylinder(R=1.) + min_z = openmc.ZPlane(z0=-0.5) + max_z = openmc.ZPlane(z0=0.5) + return +min_z & -max_z & -cylinder + + +@pytest.fixture(scope='module') +def sphere(): + sphere = openmc.Sphere(R=1.) + return -sphere + + +@pytest.fixture(scope='module') +def spherical_shell(): + sphere = openmc.Sphere(R=1.) + inner_sphere = openmc.Sphere(R=0.5) + return -sphere & +inner_sphere @pytest.fixture(scope='module') @@ -33,73 +68,68 @@ def triso_universe(): return univ -@pytest.fixture(scope='module') -def trisos(domain, triso_universe): - trisos = openmc.model.pack_trisos( - radius=_RADIUS, - fill=triso_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 - ) - return trisos - - -def test_overlap(trisos): - """Check that no TRISO particles overlap.""" - centers = [t.center for t in trisos] - +@pytest.mark.parametrize('centers', _SHAPES, indirect=True) +def test_overlap(centers): + """Check that none of the spheres in the packed configuration overlap.""" # Create KD tree for quick nearest neighbor search tree = scipy.spatial.cKDTree(centers) - # Find distance to nearest neighbor for all particles + # Find distance to nearest neighbor for all spheres d = tree.query(centers, k=2)[0] - # Get the smallest distance between any two particles + # Get the smallest distance between any two spheres 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 = 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']) +@pytest.mark.parametrize('centers,shape', zip(_SHAPES, _SHAPES), + indirect=['centers']) +def test_contained(centers, shape): + """Make sure all spheres are entirely contained within the domain.""" + if shape == 'rectangular_prism': + x = np.amax(abs(centers)) + _RADIUS + assert x < 0.5 or x == pytest.approx(0.5) - elif domain['shape'] == 'cylinder': - r = max([norm(t.center[0: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 < 0.5*domain['length'] or z == pytest.approx(0.5*domain['length']) + elif shape == 'cylinder': + r = max(np.linalg.norm(centers[:,0:2], axis=1)) + _RADIUS + z = max(abs(centers[:,2])) + _RADIUS + assert r < 1. or r == pytest.approx(1.) + assert z < 0.5 or z == pytest.approx(0.5) - elif domain['shape'] == 'sphere': - r = max([norm(t.center) for t in trisos]) + _RADIUS - assert r < domain['radius'] or r == pytest.approx(domain['radius']) + elif shape == 'sphere': + r = max(np.linalg.norm(centers, axis=1)) + _RADIUS + assert r < 1. or r == pytest.approx(1.) + + elif shape == 'spherical_shell': + d = np.linalg.norm(centers, axis=1) + r_max = max(d) + _RADIUS + r_min = min(d) - _RADIUS + assert r_max < 1. or r_max == pytest.approx(1.) + assert r_min > 0.5 or r_min == pytest.approx(0.5) -def test_packing_fraction(trisos, domain): +@pytest.mark.parametrize('centers,volume', zip(_SHAPES, _VOLUMES), + indirect=['centers']) +def test_packing_fraction(centers, volume): """Check that the actual PF is close to the requested PF.""" - pf = len(trisos)*4/3*pi*_RADIUS**3/domain['volume'] + pf = len(centers) * 4/3 * pi *_RADIUS**3 / 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 +@pytest.mark.parametrize('container', _SHAPES, indirect=True) +def test_num_spheres(container): + """Check that the function returns the correct number of spheres""" + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=container, num_spheres=50 ) - assert len(trisos) == 800 + assert len(centers) == 50 -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 +def test_triso_lattice(triso_universe, rectangular_prism): + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=rectangular_prism, pf=0.2 ) + trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c) for c in centers] lower_left = np.array((-.5, -.5, -.5)) upper_right = np.array((.5, .5, .5)) @@ -112,50 +142,29 @@ def test_triso_lattice(triso_universe): ) -def test_domain_input(triso_universe): - # Invalid domain shape +def test_container_input(triso_universe): + # Invalid container shape with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='circle' + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=+openmc.Sphere(R=1.), num_spheres=100 ) - # 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 + +def test_packing_fraction_input(sphere): + # Provide neither packing fraction nor number of spheres 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 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=sphere ) + # 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 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=sphere, pf=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 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=sphere, pf=0.5, initial_pf=0.4 )