mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 05:05:30 -04:00
Merge pull request #1121 from amandalund/triso-refactor
Update TRISO sphere packing interface and add spherical shell domain
This commit is contained in:
Paul Romano
GitHub
commit
d7359f151c
9 changed files with 808 additions and 494 deletions
|
|
@ -22,7 +22,7 @@ there are many substantial benefits to using the Python API, including:
|
|||
- Ability to plot individual universes as geometry is being created
|
||||
- A :math:`k_\text{eff}` search function (:func:`openmc.search_for_keff`)
|
||||
- Random sphere packing for generating TRISO particle locations
|
||||
(:func:`openmc.model.pack_trisos`)
|
||||
(:func:`openmc.model.pack_spheres`)
|
||||
- Ability to create materials based on natural elements or uranium enrichment
|
||||
|
||||
For those new to Python, there are many good tutorials available online. We
|
||||
|
|
|
|||
|
|
@ -37,7 +37,7 @@ Functions
|
|||
:template: myfunction.rst
|
||||
|
||||
openmc.model.create_triso_lattice
|
||||
openmc.model.pack_trisos
|
||||
openmc.model.pack_spheres
|
||||
|
||||
Model Container
|
||||
---------------
|
||||
|
|
|
|||
File diff suppressed because one or more lines are too long
File diff suppressed because it is too large
Load diff
|
|
@ -475,7 +475,7 @@ class Complement(Region):
|
|||
>>> xr = openmc.XPlane(x0=10.0)
|
||||
>>> yl = openmc.YPlane(y0=-10.0)
|
||||
>>> yr = openmc.YPlane(y0=10.0)
|
||||
>>> inside_box = +xl & -xr & +yl & -yl
|
||||
>>> inside_box = +xl & -xr & +yl & -yr
|
||||
>>> outside_box = ~inside_box
|
||||
>>> type(outside_box)
|
||||
<class 'openmc.region.Complement'>
|
||||
|
|
|
|||
|
|
@ -62,7 +62,7 @@ class Surface(IDManagerMixin):
|
|||
self.boundary_type = boundary_type
|
||||
|
||||
# A dictionary of the quadratic surface coefficients
|
||||
# Key - coefficeint name
|
||||
# Key - coefficient name
|
||||
# Value - coefficient value
|
||||
self._coefficients = {}
|
||||
|
||||
|
|
@ -1661,7 +1661,7 @@ class Quadric(Surface):
|
|||
a, b, c, d, e, f, g, h, j, k : float, optional
|
||||
coefficients for the surface. All default to 0.
|
||||
name : str, optional
|
||||
Name of the sphere. If not specified, the name will be the empty string.
|
||||
Name of the surface. If not specified, the name will be the empty string.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@
|
|||
<cell id="15" material="15" region="10 -11" universe="9" />
|
||||
<cell id="16" material="16" region="11 -12" universe="9" />
|
||||
<cell id="17" material="17" region="12" universe="9" />
|
||||
<cell fill="10" id="134" region="122 -123 124 -125 126 -127" universe="0" />
|
||||
<cell fill="10" id="18" region="13 -14 15 -16 17 -18" universe="0" />
|
||||
<cell fill="9" id="135" region="-128" translation="-0.0014062011817633224 -0.015257747167296776 -0.09184766966024316" universe="35" />
|
||||
<cell fill="9" id="136" region="-129" translation="0.10909902639777394 -0.004179177767328068 -0.046240930185910245" universe="23" />
|
||||
<cell fill="9" id="137" region="-130" translation="0.13878336652768952 -0.06900710536877097 -0.03848821497485455" universe="18" />
|
||||
|
|
@ -220,12 +220,12 @@
|
|||
<surface coeffs="0.0 0.0 0.0 0.03125" id="10" type="sphere" />
|
||||
<surface coeffs="0.0 0.0 0.0 0.03475" id="11" type="sphere" />
|
||||
<surface coeffs="0.0 0.0 0.0 0.03825" id="12" type="sphere" />
|
||||
<surface boundary="reflective" coeffs="-0.5" id="122" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="0.5" id="123" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="-0.5" id="124" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="0.5" id="125" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="-0.5" id="126" type="z-plane" />
|
||||
<surface boundary="reflective" coeffs="0.5" id="127" type="z-plane" />
|
||||
<surface boundary="reflective" coeffs="-0.5" id="13" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="0.5" id="14" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="-0.5" id="15" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="0.5" id="16" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="-0.5" id="17" type="z-plane" />
|
||||
<surface boundary="reflective" coeffs="0.5" id="18" type="z-plane" />
|
||||
<surface coeffs="-0.0014062011817633224 -0.015257747167296776 -0.09184766966024316 0.04225" id="128" type="sphere" />
|
||||
<surface coeffs="0.10909902639777394 -0.004179177767328068 -0.046240930185910245 0.04225" id="129" type="sphere" />
|
||||
<surface coeffs="0.13878336652768952 -0.06900710536877097 -0.03848821497485455 0.04225" id="130" type="sphere" />
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -10,21 +10,88 @@ 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}
|
||||
_PARAMS = [
|
||||
{'shape': 'rectangular_prism', 'volume': 1**3},
|
||||
{'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': 'sphere', 'volume': 4/3*pi*1**3},
|
||||
{'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)}
|
||||
]
|
||||
|
||||
|
||||
@pytest.fixture(scope='module', params=domain_params,
|
||||
ids=['cube', 'cylinder', 'sphere'])
|
||||
def domain(request):
|
||||
@pytest.fixture(scope='module', params=_PARAMS)
|
||||
def container(request):
|
||||
return request.param
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers(request, container):
|
||||
return request.getfixturevalue('centers_' + container['shape'])
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers_rectangular_prism():
|
||||
min_x = openmc.XPlane(x0=0)
|
||||
max_x = openmc.XPlane(x0=1)
|
||||
min_y = openmc.YPlane(y0=0)
|
||||
max_y = openmc.YPlane(y0=1)
|
||||
min_z = openmc.ZPlane(z0=0)
|
||||
max_z = openmc.ZPlane(z0=1)
|
||||
region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z
|
||||
return openmc.model.pack_spheres(radius=_RADIUS, region=region,
|
||||
pf=_PACKING_FRACTION, initial_pf=0.2)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers_x_cylinder():
|
||||
cylinder = openmc.XCylinder(R=1, y0=1, z0=2)
|
||||
min_x = openmc.XPlane(x0=0)
|
||||
max_x = openmc.XPlane(x0=1)
|
||||
region = +min_x & -max_x & -cylinder
|
||||
return openmc.model.pack_spheres(radius=_RADIUS, region=region,
|
||||
pf=_PACKING_FRACTION, initial_pf=0.2)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers_y_cylinder():
|
||||
cylinder = openmc.YCylinder(R=1, x0=1, z0=2)
|
||||
min_y = openmc.YPlane(y0=0)
|
||||
max_y = openmc.YPlane(y0=1)
|
||||
region = +min_y & -max_y & -cylinder
|
||||
return openmc.model.pack_spheres(radius=_RADIUS, region=region,
|
||||
pf=_PACKING_FRACTION, initial_pf=0.2)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers_z_cylinder():
|
||||
cylinder = openmc.ZCylinder(R=1, x0=1, y0=2)
|
||||
min_z = openmc.ZPlane(z0=0)
|
||||
max_z = openmc.ZPlane(z0=1)
|
||||
region = +min_z & -max_z & -cylinder
|
||||
return openmc.model.pack_spheres(radius=_RADIUS, region=region,
|
||||
pf=_PACKING_FRACTION, initial_pf=0.2)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers_sphere():
|
||||
sphere = openmc.Sphere(R=1, x0=1, y0=2, z0=3)
|
||||
region = -sphere
|
||||
return openmc.model.pack_spheres(radius=_RADIUS, region=region,
|
||||
pf=_PACKING_FRACTION, initial_pf=0.2)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def centers_spherical_shell():
|
||||
sphere = openmc.Sphere(R=1, x0=1, y0=2, z0=3)
|
||||
inner_sphere = openmc.Sphere(R=0.5, x0=1, y0=2, z0=3)
|
||||
region = -sphere & +inner_sphere
|
||||
return openmc.model.pack_spheres(radius=_RADIUS, region=region,
|
||||
pf=_PACKING_FRACTION, initial_pf=0.2)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def triso_universe():
|
||||
sphere = openmc.Sphere(R=_RADIUS)
|
||||
|
|
@ -33,76 +100,96 @@ 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]
|
||||
|
||||
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'])
|
||||
|
||||
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 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_contained_rectangular_prism(centers_rectangular_prism):
|
||||
"""Make sure all spheres are entirely contained within the domain."""
|
||||
d_max = np.amax(centers_rectangular_prism) + _RADIUS
|
||||
d_min = np.amin(centers_rectangular_prism) - _RADIUS
|
||||
assert d_max < 1 or d_max == pytest.approx(1)
|
||||
assert d_min > 0 or d_min == pytest.approx(0)
|
||||
|
||||
|
||||
def test_packing_fraction(trisos, domain):
|
||||
def test_contained_x_cylinder(centers_x_cylinder):
|
||||
"""Make sure all spheres are entirely contained within the domain."""
|
||||
d = np.linalg.norm(centers_x_cylinder[:,[1,2]] - [1, 2], axis=1)
|
||||
r_max = max(d) + _RADIUS
|
||||
x_max = max(centers_x_cylinder[:,0]) + _RADIUS
|
||||
x_min = min(centers_x_cylinder[:,0]) - _RADIUS
|
||||
assert r_max < 1 or r_max == pytest.approx(1)
|
||||
assert x_max < 1 or x_max == pytest.approx(1)
|
||||
assert x_min > 0 or x_min == pytest.approx(0)
|
||||
|
||||
|
||||
def test_contained_y_cylinder(centers_y_cylinder):
|
||||
"""Make sure all spheres are entirely contained within the domain."""
|
||||
d = np.linalg.norm(centers_y_cylinder[:,[0,2]] - [1, 2], axis=1)
|
||||
r_max = max(d) + _RADIUS
|
||||
y_max = max(centers_y_cylinder[:,1]) + _RADIUS
|
||||
y_min = min(centers_y_cylinder[:,1]) - _RADIUS
|
||||
assert r_max < 1 or r_max == pytest.approx(1)
|
||||
assert y_max < 1 or y_max == pytest.approx(1)
|
||||
assert y_min > 0 or y_min == pytest.approx(0)
|
||||
|
||||
|
||||
def test_contained_z_cylinder(centers_z_cylinder):
|
||||
"""Make sure all spheres are entirely contained within the domain."""
|
||||
d = np.linalg.norm(centers_z_cylinder[:,[0,1]] - [1, 2], axis=1)
|
||||
r_max = max(d) + _RADIUS
|
||||
z_max = max(centers_z_cylinder[:,2]) + _RADIUS
|
||||
z_min = min(centers_z_cylinder[:,2]) - _RADIUS
|
||||
assert r_max < 1 or r_max == pytest.approx(1)
|
||||
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."""
|
||||
d = np.linalg.norm(centers_sphere - [1, 2, 3], axis=1)
|
||||
r_max = max(d) + _RADIUS
|
||||
assert r_max < 1 or r_max == pytest.approx(1)
|
||||
|
||||
|
||||
def test_contained_spherical_shell(centers_spherical_shell):
|
||||
"""Make sure all spheres are entirely contained within the domain."""
|
||||
d = np.linalg.norm(centers_spherical_shell - [1, 2, 3], 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(container, centers):
|
||||
"""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 / container['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
|
||||
def test_num_spheres():
|
||||
"""Check that the function returns the correct number of spheres"""
|
||||
centers = openmc.model.pack_spheres(
|
||||
radius=_RADIUS, region=-openmc.Sphere(R=1), 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, centers_rectangular_prism):
|
||||
trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c)
|
||||
for c in centers_rectangular_prism]
|
||||
|
||||
lower_left = np.array((-.5, -.5, -.5))
|
||||
upper_right = np.array((.5, .5, .5))
|
||||
lower_left = np.array((0, 0, 0))
|
||||
upper_right = np.array((1, 1, 1))
|
||||
shape = (3, 3, 3)
|
||||
pitch = (upper_right - lower_left)/shape
|
||||
background = openmc.Material()
|
||||
|
|
@ -112,50 +199,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():
|
||||
# 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=-openmc.Sphere(R=1)
|
||||
)
|
||||
|
||||
# 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=-openmc.Sphere(R=1), 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=-openmc.Sphere(R=1), pf=0.5, initial_pf=0.4
|
||||
)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue