Merge remote-tracking branch 'upstream/develop' into incoming-current

This commit is contained in:
Sam Shaner 2016-08-31 07:06:02 -04:00
commit 2357dee05f
10 changed files with 1075 additions and 56 deletions

View file

@ -0,0 +1,60 @@
<?xml version="1.0" encoding="utf-8"?>
<!-- Generator: Adobe Illustrator 16.0.0, SVG Export Plug-In . SVG Version: 6.00 Build 0) -->
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg version="1.1" id="Layer_1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px"
width="257.157px" height="60px" viewBox="0 0 257.157 60" enable-background="new 0 0 257.157 60" xml:space="preserve">
<g>
<g>
<path fill="#656565" d="M81.676,48.809c-2.67,0-5.128-0.469-7.367-1.409c-2.243-0.937-4.179-2.206-5.81-3.805
c-1.632-1.601-2.901-3.477-3.807-5.638c-0.908-2.158-1.36-4.475-1.36-6.947v-0.099c0-2.472,0.46-4.788,1.385-6.945
c0.922-2.159,2.199-4.057,3.832-5.688c1.63-1.631,3.576-2.916,5.834-3.856s4.724-1.409,7.392-1.409
c2.67,0,5.125,0.469,7.367,1.409s4.179,2.21,5.812,3.808c1.63,1.6,2.899,3.479,3.805,5.637c0.907,2.161,1.361,4.476,1.361,6.947
v0.099c0,2.472-0.463,4.79-1.385,6.945c-0.922,2.161-2.2,4.059-3.832,5.688c-1.63,1.629-3.576,2.918-5.834,3.854
C86.81,48.34,84.344,48.809,81.676,48.809z M81.773,41.788c1.517,0,2.918-0.279,4.203-0.839c1.286-0.561,2.382-1.336,3.288-2.325
c0.907-0.986,1.614-2.131,2.126-3.437c0.512-1.302,0.768-2.694,0.768-4.178v-0.099c0-1.481-0.256-2.883-0.768-4.202
c-0.512-1.318-1.235-2.475-2.175-3.461c-0.94-0.991-2.053-1.771-3.338-2.35c-1.285-0.576-2.687-0.864-4.201-0.864
c-1.552,0-2.958,0.279-4.228,0.838c-1.27,0.562-2.358,1.336-3.264,2.327c-0.907,0.987-1.616,2.132-2.125,3.436
c-0.511,1.303-0.768,2.693-0.768,4.178v0.099c0,1.483,0.257,2.884,0.768,4.205c0.51,1.318,1.234,2.47,2.175,3.46
c0.939,0.988,2.042,1.772,3.313,2.347C78.815,41.499,80.224,41.788,81.773,41.788z"/>
<path fill="#656565" d="M102.145,21.715h7.516v3.808c0.924-1.253,2.035-2.283,3.338-3.09c1.301-0.809,2.942-1.211,4.92-1.211
c1.549,0,3.048,0.297,4.498,0.889c1.453,0.594,2.737,1.475,3.859,2.645c1.119,1.17,2.019,2.604,2.694,4.301
c0.675,1.7,1.013,3.652,1.013,5.859v0.1c0,2.209-0.338,4.162-1.013,5.858c-0.675,1.697-1.565,3.133-2.67,4.301
c-1.105,1.173-2.381,2.054-3.832,2.648c-1.452,0.594-2.966,0.889-4.548,0.889c-2.012,0-3.667-0.397-4.971-1.187
c-1.302-0.79-2.396-1.712-3.287-2.768v11.371h-7.516V21.715z M115.989,42.332c0.889,0,1.722-0.172,2.498-0.518
c0.773-0.348,1.458-0.843,2.051-1.484c0.592-0.641,1.064-1.408,1.409-2.298c0.347-0.89,0.52-1.896,0.52-3.018v-0.1
c0-1.087-0.173-2.082-0.52-2.99c-0.345-0.907-0.816-1.682-1.409-2.325c-0.593-0.642-1.278-1.137-2.051-1.482
c-0.776-0.345-1.609-0.52-2.498-0.52s-1.722,0.175-2.496,0.52c-0.775,0.346-1.451,0.841-2.028,1.482
c-0.577,0.644-1.037,1.418-1.385,2.325c-0.345,0.908-0.518,1.903-0.518,2.99v0.1c0,1.087,0.173,2.086,0.518,2.991
c0.348,0.906,0.808,1.684,1.385,2.324c0.577,0.642,1.253,1.137,2.028,1.484C114.267,42.16,115.101,42.332,115.989,42.332z"/>
<path fill="#656565" d="M145.62,48.809c-1.979,0-3.814-0.329-5.514-0.986c-1.696-0.661-3.162-1.6-4.399-2.816
c-1.237-1.223-2.199-2.666-2.891-4.331c-0.693-1.661-1.041-3.517-1.041-5.559v-0.102c0-1.877,0.321-3.658,0.964-5.34
c0.645-1.682,1.543-3.147,2.695-4.4c1.154-1.251,2.53-2.241,4.129-2.967c1.6-0.725,3.37-1.086,5.316-1.086
c2.206,0,4.118,0.394,5.733,1.187c1.617,0.791,2.958,1.854,4.03,3.188c1.072,1.335,1.864,2.867,2.374,4.598
c0.511,1.731,0.766,3.537,0.766,5.417c0,0.295-0.008,0.61-0.024,0.938c-0.016,0.331-0.04,0.676-0.074,1.038h-18.442
c0.364,1.716,1.112,3.009,2.25,3.881c1.138,0.873,2.547,1.313,4.228,1.313c1.251,0,2.373-0.216,3.363-0.646
c0.988-0.427,2.01-1.118,3.064-2.076l4.304,3.81c-1.254,1.547-2.771,2.76-4.552,3.633C150.12,48.372,148.026,48.809,145.62,48.809
z M150.464,32.89c-0.227-1.681-0.821-3.042-1.777-4.08c-0.956-1.037-2.226-1.557-3.807-1.557c-1.582,0-2.861,0.512-3.832,1.53
c-0.973,1.023-1.607,2.393-1.904,4.106H150.464z"/>
<path fill="#656565" d="M159.415,21.715h7.518v3.787c0.428-0.562,0.896-1.103,1.408-1.616c0.512-0.517,1.08-0.971,1.706-1.371
c0.625-0.397,1.318-0.712,2.078-0.946c0.757-0.233,1.613-0.347,2.569-0.347c2.867,0,5.084,0.873,6.65,2.619
c1.565,1.746,2.35,4.152,2.35,7.219v17.156h-7.518V33.471c0-1.776-0.395-3.118-1.186-4.021c-0.792-0.903-1.913-1.355-3.361-1.355
c-1.451,0-2.598,0.452-3.437,1.355c-0.841,0.903-1.261,2.245-1.261,4.021v14.745h-7.518V21.715z"/>
<path fill="#A31F34" d="M187.056,13.607h8.209l9.095,14.634l9.099-14.634h8.207v34.608h-7.513V25.619l-9.742,14.784h-0.196
l-9.645-14.634v22.446h-7.514V13.607z"/>
<path fill="#A31F34" d="M242.235,48.809c-2.537,0-4.893-0.459-7.071-1.385c-2.175-0.922-4.052-2.18-5.636-3.78
c-1.583-1.601-2.819-3.485-3.709-5.66c-0.888-2.179-1.333-4.501-1.333-6.974v-0.099c0-2.472,0.445-4.788,1.333-6.945
c0.89-2.159,2.126-4.057,3.709-5.688c1.584-1.631,3.477-2.916,5.687-3.856c2.207-0.94,4.647-1.409,7.317-1.409
c1.613,0,3.091,0.132,4.425,0.396c1.336,0.265,2.547,0.625,3.635,1.089c1.088,0.461,2.093,1.021,3.016,1.682
c0.923,0.657,1.781,1.385,2.57,2.175l-4.846,5.586c-1.349-1.218-2.728-2.175-4.128-2.866c-1.401-0.692-2.974-1.038-4.721-1.038
c-1.453,0-2.793,0.279-4.03,0.838c-1.235,0.562-2.299,1.336-3.188,2.327c-0.89,0.987-1.582,2.132-2.077,3.436
c-0.494,1.303-0.742,2.693-0.742,4.178v0.099c0,1.483,0.248,2.884,0.742,4.205c0.495,1.318,1.178,2.47,2.053,3.46
c0.872,0.988,1.929,1.772,3.164,2.347c1.235,0.576,2.594,0.865,4.079,0.865c1.979,0,3.65-0.362,5.02-1.086
c1.366-0.726,2.726-1.714,4.079-2.968l4.844,4.895c-0.89,0.958-1.814,1.814-2.768,2.571c-0.958,0.758-2.003,1.411-3.142,1.953
c-1.136,0.544-2.379,0.959-3.733,1.236C245.433,48.669,243.915,48.809,242.235,48.809z"/>
</g>
<path fill="#A31F34" d="M30.614,1.024c-10.866,0-20.325,5.991-25.284,14.845h20.495l13.167,22.804l-11.609,20.11
c1.062,0.119,2.139,0.192,3.231,0.192c16.003,0,28.975-12.972,28.975-28.976C59.589,13.997,46.617,1.024,30.614,1.024z"/>
<path fill="#656565" d="M16.416,55.243h5.807l9.567-16.57l-9.567-16.566H3.093l-0.607,1.052C1.954,25.355,1.639,27.638,1.639,30
C1.639,40.84,7.603,50.274,16.416,55.243z"/>
</g>
</svg>

After

Width:  |  Height:  |  Size: 5.9 KiB

View file

@ -334,6 +334,7 @@ Functions
:nosignatures:
openmc.model.create_triso_lattice
openmc.model.pack_trisos
--------------------------------------------
:mod:`openmc.data` -- Nuclear Data Interface

View file

@ -2105,8 +2105,8 @@ attributes or sub-elements. These are not used in "voxel" plots:
*Default*: None
:meshlines:
The ``meshlines`` sub-element allows for plotting the boundaries of
a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per
The ``meshlines`` sub-element allows for plotting the boundaries of a
regular mesh on top of a plot. Only one ``meshlines`` element is allowed per
``plot`` element, and it must contain as attributes or sub-elements a mesh
type and a linewidth. Optionally, a color may be specified for the overlay:

View file

@ -337,7 +337,7 @@ class MGXS(object):
if self.by_nuclide:
return self.get_nuclides()
else:
return 'sum'
return ['sum']
@property
def loaded_sp(self):
@ -1483,25 +1483,27 @@ class MGXS(object):
if self.by_nuclide and nuclides == 'sum':
# Use tally summation to sum across all nuclides
query_nuclides = self.get_nuclides()
xs_tally = self.xs_tally.summation(nuclides=query_nuclides)
query_nuclides = [nuclides]
xs_tally = self.xs_tally.summation(nuclides=self.get_nuclides())
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
# Remove nuclide column since it is homogeneous and redundant
if self.domain_type == 'mesh':
df.drop('nuclide', axis=1, level=0, inplace=True)
df.drop('sum(nuclide)', axis=1, level=0, inplace=True)
else:
df.drop('nuclide', axis=1, inplace=True)
df.drop('sum(nuclide)', axis=1, inplace=True)
# If the user requested a specific set of nuclides
elif self.by_nuclide and nuclides != 'all':
query_nuclides = nuclides
xs_tally = self.xs_tally.get_slice(nuclides=nuclides)
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
# If the user requested all nuclides, keep nuclide column in dataframe
else:
query_nuclides = self.nuclides
df = self.xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
@ -1513,7 +1515,7 @@ class MGXS(object):
# Override energy groups bounds with indices
all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
all_groups = np.repeat(all_groups, self.num_nuclides)
all_groups = np.repeat(all_groups, len(query_nuclides))
if 'energy low [MeV]' in df and 'energyout low [MeV]' in df:
df.rename(columns={'energy low [MeV]': 'group in'},
inplace=True)

View file

@ -1,13 +1,26 @@
from __future__ import division
import copy
from collections import Iterable
from numbers import Real
import warnings
import itertools
import random
from collections import Iterable, defaultdict
from numbers import Real
from random import uniform, gauss
from heapq import heappush, heappop
from math import pi, sin, cos, floor, log10, sqrt
from abc import ABCMeta, abstractproperty, abstractmethod
import numpy as np
try:
import scipy.spatial
_SCIPY_AVAILABLE = True
except ImportError:
_SCIPY_AVAILABLE = False
import openmc
import openmc.checkvalue as cv
class TRISO(openmc.Cell):
"""Tristructural-isotopic (TRISO) micro fuel particle
@ -82,6 +95,377 @@ class TRISO(openmc.Cell):
k_min:k_max+1, j_min:j_max+1, i_min:i_max+1]))
class _Domain(object):
"""Container in which to pack particles.
Parameters
----------
particle_radius : float
Radius of particles to be packed in container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
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.
volume : float
Volume of the container.
"""
__metaclass__ = ABCMeta
def __init__(self, particle_radius, center=[0., 0., 0.]):
self._cell_length = None
self._limits = None
self.particle_radius = particle_radius
self.center = center
@property
def particle_radius(self):
return self._particle_radius
@property
def center(self):
return self._center
@abstractproperty
def limits(self):
pass
@abstractproperty
def cell_length(self):
pass
@abstractproperty
def volume(self):
pass
@particle_radius.setter
def particle_radius(self, particle_radius):
self._particle_radius = float(particle_radius)
self._limits = None
self._cell_length = None
@center.setter
def center(self, center):
if np.asarray(center).size != 3:
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
which the given particle center falls.
Parameters
----------
p : Iterable of float
Cartesian coordinates of particle center.
Returns
-------
tuple of int
Indices of mesh cell.
"""
return tuple(int(p[i]/self.cell_length[i]) for i in range(3))
def nearby_mesh_cells(self, p):
"""Calculates the indices of all cells in a mesh overlaid on the domain
within one diameter of the given particle.
Parameters
----------
p : Iterable of float
Cartesian coordinates of particle center.
Returns
-------
list of tuple of int
Indices of mesh cells.
"""
d = 2*self.particle_radius
r = [[a/self.cell_length[i] for a in [p[i]-d, p[i], p[i]+d]]
for i in range(3)]
return list(itertools.product(*({int(x) for x in y} for y in r)))
@abstractmethod
def random_point(self):
"""Generate Cartesian coordinates of center of a particle that is
contained entirely within the domain with uniform probability.
Returns
-------
list of float
Cartesian coordinates of particle center.
"""
pass
class _CubicDomain(_Domain):
"""Cubic container in which to pack particles.
Parameters
----------
length : float
Length of each side of the cubic container.
particle_radius : float
Radius of particles to be packed in container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
length : float
Length of each side of the cubic container.
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
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.
volume : float
Volume of the container.
"""
def __init__(self, length, particle_radius, center=[0., 0., 0.]):
super(_CubicDomain, self).__init__(particle_radius, center)
self.length = length
@property
def length(self):
return self._length
@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]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [self.length, self.length, self.length]
self._cell_length = [x/int(x/(4*self.particle_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return self.length**3
@length.setter
def length(self, length):
self._length = float(length)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
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])]
class _CylindricalDomain(_Domain):
"""Cylindrical container in which to pack particles.
Parameters
----------
length : float
Length along z-axis of the cylindrical container.
radius : float
Radius of the cylindrical container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
length : float
Length along z-axis of the cylindrical container.
radius : float
Radius of the cylindrical container.
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
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.
volume : float
Volume of the container.
"""
def __init__(self, length, radius, particle_radius, center=[0., 0., 0.]):
super(_CylindricalDomain, self).__init__(particle_radius, center)
self.length = length
self.radius = radius
@property
def length(self):
return self._length
@property
def radius(self):
return self._radius
@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]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [2*self.radius, 2*self.radius, self.length]
self._cell_length = [x/int(x/(4*self.particle_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return self.length * pi * self.radius**2
@length.setter
def length(self, length):
self._length = float(length)
self._limits = None
self._cell_length = None
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
self._limits = limits
def random_point(self):
r = sqrt(uniform(0, (self.radius - self.particle_radius)**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])]
class _SphericalDomain(_Domain):
"""Spherical container in which to pack particles.
Parameters
----------
radius : float
Radius of the spherical container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
radius : float
Radius of the spherical container.
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
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.
volume : float
Volume of the container.
"""
def __init__(self, radius, particle_radius, center=[0., 0., 0.]):
super(_SphericalDomain, self).__init__(particle_radius, center)
self.radius = radius
@property
def radius(self):
return self._radius
@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]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [2*self.radius, 2*self.radius, 2*self.radius]
self._cell_length = [x/int(x/(4*self.particle_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return 4/3 * pi * self.radius**3
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
self._limits = limits
def random_point(self):
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)]
def create_triso_lattice(trisos, lower_left, pitch, shape, background):
"""Create a lattice containing TRISO particles for optimized tracking.
@ -153,3 +537,503 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
lattice.outer = openmc.Universe(cells=[background_cell])
return lattice
def _random_sequential_pack(domain, n_particles):
"""Random sequential packing of particles within a container.
Parameters
----------
domain : openmc.model._Domain
Container in which to pack particles.
n_particles : int
Number of particles to pack.
Returns
------
numpy.ndarray
Cartesian coordinates of centers of particles.
"""
sqd = (2*domain.particle_radius)**2
particles = []
mesh = defaultdict(list)
for i in range(n_particles):
# Randomly sample new center coordinates while there are any overlaps
while True:
p = domain.random_point()
idx = domain.mesh_cell(p)
if any((p[0]-q[0])**2 + (p[1]-q[1])**2 + (p[2]-q[2])**2 < sqd
for q in mesh[idx]):
continue
else:
break
particles.append(p)
for idx in domain.nearby_mesh_cells(p):
mesh[idx].append(p)
return np.array(particles)
def _close_random_pack(domain, particles, contraction_rate):
"""Close random packing of particles using the Jodrey-Tory algorithm.
Parameters
----------
domain : openmc.model._Domain
Container in which to pack particles.
particles : numpy.ndarray
Initial Cartesian coordinates of centers of particles.
contraction_rate : float
Contraction rate of outer diameter.
"""
def add_rod(d, i, j):
"""Add a new rod to the priority queue.
Parameters
----------
d : float
distance between centers of particles i and j.
i, j : int
Index of particles in particles array.
"""
rod = [d, i, j]
rods_map[i] = (j, rod)
rods_map[j] = (i, rod)
heappush(rods, rod)
def remove_rod(i):
"""Mark the rod containing particle i as removed.
Parameters
----------
i : int
Index of particle in particles array.
"""
if i in rods_map:
j, rod = rods_map.pop(i)
del rods_map[j]
rod[1] = removed
rod[2] = removed
def pop_rod():
"""Remove and return the shortest rod.
Returns
-------
d : float
distance between centers of particles i and j.
i, j : int
Index of particles in particles array.
"""
while rods:
d, i, j = heappop(rods)
if i != removed and j != removed:
del rods_map[i]
del rods_map[j]
return d, i, j
def create_rod_list():
"""Generate sorted list of rods (distances between particle centers).
Rods are arranged in a heap where each element contains the rod length
and the particle indices. A rod between particles p and q is only
included if the distance between p and q could not be changed by the
elimination of a greater overlap, i.e. q has no nearer neighbors than p.
A mapping of particle ids to rods is maintained in 'rods_map'. Each key
in the dict is the id of a particle that is in the rod list, and the
value is the id of its nearest neighbor and the rod that contains them.
The dict is used to find rods in the priority queue and to mark removed
rods so rods can be "removed" without breaking the heap structure
invariant.
"""
# Create KD tree for quick nearest neighbor search
tree = scipy.spatial.cKDTree(particles)
# 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]
# Array of particle indices, indices of nearest neighbors, and
# distances to nearest neighbors
a = np.vstack((list(range(n.size)), n, d)).T
# 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]]
# Find the intersection between 'a' and 'b': a list of particles who
# are each other's nearest neighbors and the distance between them
r = list({tuple(x) for x in a} & {tuple(x) for x in b})
# Remove duplicate rods and sort by distance
r = map(list, set([(x[2], int(min(x[0:2])), int(max(x[0:2])))
for x in r]))
# Clear priority queue and add rods
del rods[:]
rods_map.clear()
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
center coordinates.
'mesh'/'mesh_map' is a two way dictionary used to look up which
particles are located within one diameter of a given mesh cell and
which mesh cells a given particle center is within one diameter of.
This is used to speed up the nearest neighbor search.
Parameters
----------
i : int
Index of particle 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]
# 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:
d_out^(i+1) = d_out^(i) - (1/2)^(j) * d_out0 * k / n,
where k is the contraction rate, n is the number of particles, and
j = floor(-log10(pf_out - pf_in)).
"""
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)
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 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])
update_mesh(i)
update_mesh(j)
def nearest(i):
"""Find index of nearest neighbor of particle i.
Parameters
----------
i : int
Index in particles array of particle for which to find nearest
neighbor.
Returns
-------
int
Index in particles array of nearest neighbor of i
float
distance between i and nearest neighbor.
"""
# Need the second nearest neighbor of i since the nearest neighbor
# will be itself. Using argpartition, the k-th nearest neighbor is
# placed at index k.
idx = list(mesh[domain.mesh_cell(particles[i])])
dists = scipy.spatial.distance.cdist([particles[i]], particles[idx])[0]
if dists.size > 1:
j = dists.argpartition(1)[1]
return idx[j], dists[j]
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.
Parameters
----------
i, j : int
Index 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]
if not _SCIPY_AVAILABLE:
raise ImportError('SciPy must be installed to perform '
'close random packing.')
n_particles = len(particles)
diameter = 2*domain.particle_radius
# Flag for marking rods that have been removed from priority queue
removed = -1
# Outer diameter initially set to arbitrary value that yields pf of 1
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]
rods = []
rods_map = {}
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:
create_rod_list()
if inner_diameter[0] >= 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:
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):
"""Generate a random, non-overlapping configuration of TRISO particles
within a container.
Parameters
----------
radius : float
Outer radius of TRISO particles.
fill : openmc.Universe
Universe which contains all layers of the TRISO particle.
domain_shape : {'cube', 'cylinder', or 'sphere'}
Geometry of the container in which the TRISO particles are packed.
domain_length : float
Length of the container (if cube or cylinder).
domain_radius : float
Radius of the container (if cylinder or sphere).
domain_center : Iterable of float
Cartesian coordinates of the center of the container.
n_particles : int
Number of TRISO particles to pack in the domain. Exactly one of
'n_particles' and 'packing_fraction' should be specified -- the other
will be calculated.
packing_fraction : float
Packing fraction of particles. Exactly one of 'n_particles' and
'packing_fraction' should be specified -- the other will be calculated.
initial_packing_fraction : float, optional
Packing fraction used to initialize the configuration of particles in
the domain. Default value is 0.3. It is not recommended to set the
initial packing fraction much higher than 0.3 as the random sequential
packing algorithm becomes prohibitively slow as it approaches its limit
(~0.38).
contraction_rate : float, optional
Contraction rate of outer diameter. This can affect the speed of the
close random packing algorithm. Default value is 1/400.
seed : int, optional
RNG seed.
Returns
-------
trisos : list of openmc.model.TRISO
List of TRISO particles in the domain.
Notes
-----
The particle configuration is generated using a combination of random
sequential packing (RSP) and close random packing (CRP). RSP performs
better than CRP for lower packing fractions (pf), but it becomes
prohibitively slow as it approaches its packing limit (~0.38). CRP can
achieve higher pf of up to ~0.64 and scales better with increasing pf.
If the desired pf is below some threshold for which RSP will be faster than
CRP ('initial_packing_fraction'), only RSP is used. If a higher pf is
required, particles with a radius smaller than the desired final radius
(and therefore with a smaller pf) are initialized within the domain using
RSP. This initial configuration of particles is then used as a starting
point for CRP using Jodrey and Tory's algorithm [1]_.
In RSP, particle centers are placed one by one at random, and placement
attempts for a particle are made until the particle is not overlapping any
others. This implementation of the algorithm uses a mesh over the domain
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,
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
between particles based on outer diameter is eliminated by moving the
particles apart along the line joining their centers. Iterations continue
until the two diameters converge or until the desired pf is reached.
References
----------
.. [1] W. S. Jodrey and E. M. Tory, "Computer simulation of close random
packing of equal spheres", Phys. Rev. A 32 (1985) 2347-2351.
"""
# Check for valid container geometry and dimensions
if domain_shape not in ['cube', 'cylinder', 'sphere']:
raise ValueError('Unable to set domain_shape to "{}". Only "cube", '
'"cylinder", and "sphere" are '
'supported."'.format(domain_shape))
if not domain_length and domain_shape in ['cube', 'cylinder']:
raise ValueError('"domain_length" must be specified for {} domain '
'geometry '.format(domain_shape))
if not domain_radius and domain_shape in ['cylinder', 'sphere']:
raise ValueError('"domain_radius" must be specified for {} domain '
'geometry '.format(domain_shape))
if domain_shape is 'cube':
domain = _CubicDomain(length=domain_length, particle_radius=radius,
center=domain_center)
elif domain_shape is 'cylinder':
domain = _CylindricalDomain(length=domain_length, radius=domain_radius,
particle_radius=radius, center=domain_center)
elif domain_shape is 'sphere':
domain = _SphericalDomain(radius=domain_radius, particle_radius=radius,
center=domain_center)
# Calculate the packing fraction if the number of particles is specified;
# otherwise, calculate the number of particles from the packing fraction.
if ((n_particles is None and packing_fraction is None) or
(n_particles is not None and packing_fraction is not None)):
raise ValueError('Exactly one of "n_particles" and "packing_fraction" '
'must be specified.')
elif packing_fraction is None:
n_particles = int(n_particles)
packing_fraction = 4/3*pi*radius**3*n_particles / domain.volume
elif n_particles is None:
packing_fraction = float(packing_fraction)
n_particles = int(packing_fraction*domain.volume // (4/3*pi*radius**3))
# Check for valid packing fractions for each algorithm
if packing_fraction >= 0.64:
raise ValueError('Packing fraction of {} is greater than the '
'packing fraction limit for close random '
'packing (0.64)'.format(packing_fraction))
if initial_packing_fraction >= 0.38:
raise ValueError('Initial packing fraction of {} is greater than the '
'packing fraction limit for random sequential'
'packing (0.38)'.format(initial_packing_fraction))
if initial_packing_fraction > packing_fraction:
initial_packing_fraction = packing_fraction
if packing_fraction > 0.3:
initial_packing_fraction = 0.3
random.seed(seed)
# Calculate the particle radius used in the initial random sequential
# packing from the initial packing fraction
initial_radius = (3/4 * initial_packing_fraction * domain.volume /
(pi * n_particles))**(1/3)
domain.particle_radius = initial_radius
# 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]]]
# Generate non-overlapping particles for an initial inner radius using
# random sequential packing algorithm
particles = _random_sequential_pack(domain, n_particles)
# Use the particle configuration produced in random sequential packing as a
# starting point for close random pack with the desired final particle
# radius
if initial_packing_fraction != packing_fraction:
domain.particle_radius = radius
_close_random_pack(domain, particles, contraction_rate)
trisos = []
for p in particles:
trisos.append(TRISO(radius, fill, p))
return trisos

View file

@ -67,6 +67,11 @@ class Plot(object):
col_spec : dict
Dictionary indicating that certain cells/materials (keys) should be
colored with a specific RGB (values)
level : int
Universe depth to plot at
meshlines : dict
Dictionary defining type, id, linewidth and color of a regular mesh
to be plotted on top of a plot
"""
@ -81,10 +86,12 @@ class Plot(object):
self._color = 'cell'
self._type = 'slice'
self._basis = 'xy'
self._background = [0, 0, 0]
self._background = None
self._mask_components = None
self._mask_background = None
self._col_spec = None
self._level = None
self._meshlines = None
@property
def id(self):
@ -138,6 +145,14 @@ class Plot(object):
def col_spec(self):
return self._col_spec
@property
def level(self):
return self._level
@property
def meshlines(self):
return self._meshlines
@id.setter
def id(self, plot_id):
if plot_id is None:
@ -231,9 +246,9 @@ class Plot(object):
@mask_components.setter
def mask_components(self, mask_components):
cv.check_type('plot mask_components', mask_components, Iterable, Integral)
cv.check_type('plot mask components', mask_components, Iterable, Integral)
for component in mask_components:
cv.check_greater_than('plot mask_components', component, 0, True)
cv.check_greater_than('plot mask components', component, 0, True)
self._mask_components = mask_components
@mask_background.setter
@ -245,6 +260,45 @@ class Plot(object):
cv.check_less_than('plot mask background', rgb, 256)
self._mask_background = mask_background
@level.setter
def level(self, plot_level):
cv.check_type('plot level', plot_level, Integral)
cv.check_greater_than('plot level', plot_level, 0, equality=True)
self._level = plot_level
@meshlines.setter
def meshlines(self, meshlines):
cv.check_type('plot meshlines', meshlines, dict)
if 'type' not in meshlines:
msg = 'Unable to set on plot the meshlines "{0}" which ' \
'does not have a "type" key'.format(meshlines)
raise ValueError(msg)
elif meshlines['type'] not in ['tally', 'entropy', 'ufs', 'cmfd']:
msg = 'Unable to set the meshlines with ' \
'type "{0}"'.format(meshlines['type'])
raise ValueError(msg)
if 'id' in meshlines:
cv.check_type('plot meshlines id', meshlines['id'], Integral)
cv.check_greater_than('plot meshlines id', meshlines['id'], 0,
equality=True)
if 'linewidth' in meshlines:
cv.check_type('plot mesh linewidth', meshlines['linewidth'], Integral)
cv.check_greater_than('plot mesh linewidth', meshlines['linewidth'],
0, equality=True)
if 'color' in meshlines:
cv.check_type('plot meshlines color', meshlines['color'], Iterable,
Integral)
cv.check_length('plot meshlines color', meshlines['color'], 3)
for rgb in meshlines['color']:
cv.check_greater_than('plot meshlines color', rgb, 0, True)
cv.check_less_than('plot meshlines color', rgb, 256)
self._meshlines = meshlines
def __repr__(self):
string = 'Plot\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@ -256,11 +310,16 @@ class Plot(object):
string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin)
string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin)
string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color)
string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t',
string += '{0: <16}{1}{2}\n'.format('\tBackground', '=\t',
self._background)
string += '{0: <16}{1}{2}\n'.format('\tMask components', '=\t',
self._mask_components)
string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t',
string += '{0: <16}{1}{2}\n'.format('\tMask background', '=\t',
self._mask_background)
string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec)
string += '{0: <16}{1}{2}\n'.format('\tLevel', '=\t', self._level)
string += '{0: <16}{1}{2}\n'.format('\tMeshlines', '=\t',
self._meshlines)
return string
def colorize(self, geometry, seed=1):
@ -382,7 +441,7 @@ class Plot(object):
subelement = ET.SubElement(element, "pixels")
subelement.text = ' '.join(map(str, self._pixels))
if self._mask_background is not None:
if self._background is not None:
subelement = ET.SubElement(element, "background")
subelement.text = ' '.join(map(str, self._background))
@ -400,6 +459,21 @@ class Plot(object):
subelement.set("background", ' '.join(map(
str, self._mask_background)))
if self._level is not None:
subelement = ET.SubElement(element, "level")
subelement.text = str(self._level)
if self._meshlines is not None:
subelement = ET.SubElement(element, "meshlines")
subelement.set("meshtype", self._meshlines['type'])
if self._meshlines['id'] is not None:
subelement.set("id", str(self._meshlines['id']))
if self._meshlines['linewidth'] is not None:
subelement.set("linewidth", str(self._meshlines['linewidth']))
if self._meshlines['color'] is not None:
subelement.set("color", ' '.join(map(
str, self._meshlines['color'])))
return element

View file

@ -38,43 +38,56 @@ contains
use omp_lib
#endif
write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
' .d88888b. 888b d888 .d8888b.', &
' d88P" "Y88b 8888b d8888 d88P Y88b', &
' 888 888 88888b.d88888 888 888', &
' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', &
' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', &
' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', &
' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', &
' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', &
'__________________888______________________________________________________', &
' 888', &
' 888'
write(UNIT=OUTPUT_UNIT, FMT='(/23(A/))') &
' %%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' ############### %%%%%%%%%%%%%%%%%%%%%%%%', &
' ################## %%%%%%%%%%%%%%%%%%%%%%%', &
' ################### %%%%%%%%%%%%%%%%%%%%%%%', &
' #################### %%%%%%%%%%%%%%%%%%%%%%', &
' ##################### %%%%%%%%%%%%%%%%%%%%%', &
' ###################### %%%%%%%%%%%%%%%%%%%%', &
' ####################### %%%%%%%%%%%%%%%%%%', &
' ####################### %%%%%%%%%%%%%%%%%', &
' ###################### %%%%%%%%%%%%%%%%%', &
' #################### %%%%%%%%%%%%%%%%%', &
' ################# %%%%%%%%%%%%%%%%%', &
' ############### %%%%%%%%%%%%%%%%', &
' ############ %%%%%%%%%%%%%%%', &
' ######## %%%%%%%%%%%%%%', &
' %%%%%%%%%%%'
! Write version information
write(UNIT=OUTPUT_UNIT, FMT=*) &
' Copyright: 2011-2016 Massachusetts Institute of Technology'
' | The OpenMC Monte Carlo Code'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License: http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') &
' Copyright | 2011-2016 Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License | http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I1,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(10X,"Git SHA1 | ",A)') GIT_SHA1
#endif
! Write the date and time
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') &
time_stamp()
write(UNIT=OUTPUT_UNIT, FMT='(9X,"Date/Time | ",A)') time_stamp()
#ifdef MPI
! Write number of processors
write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') &
write(UNIT=OUTPUT_UNIT, FMT='(5X,"MPI Processes | ",A)') &
trim(to_str(n_procs))
#endif
#ifdef _OPENMP
! Write number of OpenMP threads
write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') &
write(UNIT=OUTPUT_UNIT, FMT='(4X,"OpenMP Threads | ",A)') &
trim(to_str(omp_get_max_threads()))
#endif

View file

@ -1 +1 @@
f33e6653b883200457df2ff2ba9cf715d5ddaa1296dd71d277c6f1d9d5b7831cc92aaf1e97509d26e5a93235cd9f775c0cfaa5ebc3dfe8fc71469bac166d362b
792d82b08d5fa6ac19df668d84cb90cbbe178e8769c87b732e9616ffeb70d8cfb4096607e58cda749cb800bb48139ac72f6983e84f38237d7d55711bdd1c6d4d

View file

@ -1,2 +1,2 @@
k-combined:
1.662675E+00 1.475968E-02
1.636336E+00 1.154000E-01

View file

@ -60,24 +60,9 @@ class TRISOTestHarness(PyAPITestHarness):
inner_univ = openmc.Universe(cells=[c1, c2, c3, c4, c5])
outer_radius = 422.5*1e-4
trisos = []
random.seed(1)
for i in range(100):
# Randomly sample location
lim = 0.5 - outer_radius*1.001
x = random.uniform(-lim, lim)
y = random.uniform(-lim, lim)
z = random.uniform(-lim, lim)
t = openmc.model.TRISO(outer_radius, inner_univ, (x, y, z))
# Make sure TRISO doesn't overlap with another
for tp in trisos:
xp, yp, zp = tp.center
distance = sqrt((x - xp)**2 + (y - yp)**2 + (z - zp)**2)
if distance <= 2*outer_radius:
break
else:
trisos.append(t)
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
min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')