diff --git a/docs/source/_images/openmc.png b/docs/source/_images/openmc.png deleted file mode 100644 index 9f5e97cd6e..0000000000 Binary files a/docs/source/_images/openmc.png and /dev/null differ diff --git a/docs/source/_images/openmc200px.png b/docs/source/_images/openmc200px.png deleted file mode 100644 index 3997c6baa0..0000000000 Binary files a/docs/source/_images/openmc200px.png and /dev/null differ diff --git a/docs/source/_images/openmc_logo.png b/docs/source/_images/openmc_logo.png new file mode 100644 index 0000000000..73d4765387 Binary files /dev/null and b/docs/source/_images/openmc_logo.png differ diff --git a/docs/source/_images/openmc_logo.svg b/docs/source/_images/openmc_logo.svg new file mode 100644 index 0000000000..a7352b79ab --- /dev/null +++ b/docs/source/_images/openmc_logo.svg @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + diff --git a/docs/source/_static/theme_overrides.css b/docs/source/_static/theme_overrides.css index bee03f4150..dea941814d 100644 --- a/docs/source/_static/theme_overrides.css +++ b/docs/source/_static/theme_overrides.css @@ -16,3 +16,7 @@ .wy-table, .rst-content table.docutils, .rst-content table.field-list { margin-bottom: 0px; } + +.wy-side-nav-search { + background-color: #343131; +} diff --git a/docs/source/conf.py b/docs/source/conf.py index 1baea2b03c..75e621bc17 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -129,7 +129,7 @@ if not on_rtd: html_theme = 'sphinx_rtd_theme' html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] -html_logo = '_images/openmc200px.png' +html_logo = '_images/openmc_logo.png' # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". diff --git a/docs/source/pythonapi/examples/mdgxs-part-i.rst b/docs/source/pythonapi/examples/mdgxs-part-i.rst new file mode 100644 index 0000000000..953dcf4700 --- /dev/null +++ b/docs/source/pythonapi/examples/mdgxs-part-i.rst @@ -0,0 +1,13 @@ +.. _notebook_mdgxs_part_i: + +========================== +MDGXS Part I: Introduction +========================== + +.. only:: html + + .. notebook:: mdgxs-part-i.ipynb + +.. only:: latex + + IPython notebooks must be viewed in the online HTML documentation. diff --git a/docs/source/pythonapi/examples/mdgxs-part-ii.rst b/docs/source/pythonapi/examples/mdgxs-part-ii.rst new file mode 100644 index 0000000000..a42eb766b7 --- /dev/null +++ b/docs/source/pythonapi/examples/mdgxs-part-ii.rst @@ -0,0 +1,13 @@ +.. _notebook_mdgxs_part_ii: + +================================ +MDGXS Part II: Advanced Features +================================ + +.. only:: html + + .. notebook:: mdgxs-part-ii.ipynb + +.. only:: latex + + IPython notebooks must be viewed in the online HTML documentation. diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 392c2c72df..14f4a2128d 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -334,6 +334,7 @@ Functions :nosignatures: openmc.model.create_triso_lattice + openmc.model.pack_trisos -------------------------------------------- :mod:`openmc.data` -- Nuclear Data Interface diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 41623b6a16..e9126f16b3 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -2135,8 +2135,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: diff --git a/openmc/mesh.py b/openmc/mesh.py index 58b9c7c0e5..7d7b483f73 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -187,15 +187,20 @@ class Mesh(object): of the mesh. For example the following code: - for mesh_index in mymesh.cell_generator(): - print mesh_index - will produce the following output for a 3-D 2x2x2 mesh in mymesh: - [1, 1, 1] - [1, 1, 2] - [1, 2, 1] - [1, 2, 2] - ... + .. code-block:: python + + for mesh_index in mymesh.cell_generator(): + print mesh_index + + will produce the following output for a 3-D 2x2x2 mesh in mymesh:: + + [1, 1, 1] + [1, 1, 2] + [1, 2, 1] + [1, 2, 2] + ... + """ diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index de7475308c..1d2e9098d3 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -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) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 89e0d8aa76..5525ea559c 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -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 diff --git a/openmc/plots.py b/openmc/plots.py index 73b51da5e8..cc5c0d44b3 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -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 diff --git a/src/output.F90 b/src/output.F90 index f8b0c5ab09..b8bd804259 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index 3e29529450..4f58f14fbb 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -1 +1 @@ -b59af664a4db28471fbf7a19514eb1d61c9c890d2bfa83d23ed186ad0f23d8fa2115fb56c18212476183bcab78622622d8394ba97a19e300e0ed84acecdce3ca \ No newline at end of file +033b09236ffceab5f0f7860d4926e0a66f328c6dd6fa48b28d6237278d42f9e06a28f2368f273773e417b8177705bc9a5e3ef67335734e5f17ff9843b72a2357 \ No newline at end of file diff --git a/tests/test_triso/results_true.dat b/tests/test_triso/results_true.dat index ea7da21edf..15107e8c84 100644 --- a/tests/test_triso/results_true.dat +++ b/tests/test_triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.662675E+00 1.475968E-02 +1.636336E+00 1.154000E-01 diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py index 04bd1322a9..386be00b36 100644 --- a/tests/test_triso/test_triso.py +++ b/tests/test_triso/test_triso.py @@ -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')