From f3cfc181f47d959158d5e871209c45b9aa810f33 Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 16 Aug 2016 20:14:34 -0500 Subject: [PATCH 01/11] Add packing function for TRISO particles --- openmc/model/triso.py | 851 +++++++++++++++++++++++++++++- tests/test_triso/inputs_true.dat | 2 +- tests/test_triso/results_true.dat | 2 +- tests/test_triso/test_triso.py | 47 +- 4 files changed, 868 insertions(+), 34 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 89e0d8aa7..c1ed9336a 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -1,13 +1,22 @@ +from __future__ import division import copy -from collections import Iterable +from collections import Iterable, defaultdict from numbers import Real import warnings +import itertools +import scipy.spatial +from scipy.spatial.distance import cdist +import random +from random import uniform, gauss +from heapq import heappush, heappop +from math import pi, sin, cos, floor, log10 import numpy as np import openmc import openmc.checkvalue as cv + class TRISO(openmc.Cell): """Tristructural-isotopic (TRISO) micro fuel particle @@ -153,3 +162,843 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): lattice.outer = openmc.Universe(cells=[background_cell]) return lattice + + +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 is faster 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 performs better + 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 rondom, and placement + attempts for a particles are made until the particle is not overlapping any + others. This implementation of the algorithm uses a lattice over the domain + to speed up the nearest neighbor search by only searching for a particle's + neighbors within that lattice 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 and the outer diameter + is decreased. 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. + + """ + + def get_domain_volume(): + """Calculates the volume of the container in which the TRISO particles + are packed. + + Returns + ------- + float + Volume of the domain. + + """ + + if domain_shape is 'cube': + return domain_length**3 + elif domain_shape is 'cylinder': + return domain_length * pi * domain_radius**2 + elif domain_shape is 'sphere': + return 4/3 * pi * domain_radius**3 + + + def get_cell_length(radius): + """Calculates the length of a lattice element in x-, y-, and + z-directions. + + Parameters + ---------- + radius : float + Radius of the particle. + + Returns + ------- + tuple of float + Length of lattice cell in x-, y-, and z-directions. + + """ + + if domain_length: + m = domain_length/int(domain_length/(4*radius)) + if domain_radius: + n = 2*domain_radius/int(domain_radius/(2*radius)) + + if domain_shape is 'cube': + return (m, m, m) + elif domain_shape is 'cylinder': + return (n, n, m) + elif domain_shape is 'sphere': + return (n, n, n) + + + def get_boundary_extremes(): + """Calculates the minimum and maximum positions in x-, y-, and + z-directions where a particle center can be placed within the domain. + + Returns + ------- + llim, ulim : tuple of float + Minimum and maximum position in x-, y-, and z-directions where + particle center can be placed. + + """ + + if domain_length: + x_min = radius + x_max = domain_length - radius + if domain_radius: + r_min = radius - domain_radius + r_max = domain_radius - radius + + if domain_shape is 'cube': + return (x_min, x_min, x_min), (x_max, x_max, x_max) + elif domain_shape is 'cylinder': + return (r_min, r_min, x_min), (r_max, r_max, x_max) + elif domain_shape is 'sphere': + return (r_min, r_min, r_min), (r_max, r_max, r_max) + + + def get_particle_offset(): + """Calculates the offset in x-, y-, and z-directions of the particle + center based on the domain center + + Returns + ------- + tuple of float + Amount to offset particle center in x-, y-, and z-directions + + """ + + if domain_shape is 'cube': + return np.array(domain_center) - 3*(domain_length/2,) + elif domain_shape is 'cylinder': + return np.array(domain_center) - (0, 0, domain_length/2) + elif domain_shape is 'sphere': + return np.array(domain_center) + + + def inner_packing_fraction(): + """Calculates the true packing fraction of the particles based on the + inner diameter. + + Returns + ------- + float + Packing fraction calculated from inner diameter. + + """ + + return (4/3 * pi * (inner_diameter[0]/2)**3 * n_particles / + domain_volume) + + + def outer_packing_fraction(): + """Calculates the nominal packing fraction of the particles based on + the outer diameter. + + Returns + ------- + float + Packing fraction calculated from outer diameter. + + """ + + return (4/3 * pi * (outer_diameter[0]/2)**3 * n_particles / + domain_volume) + + + def random_point_cube(): + """Generate Cartesian coordinates of center of a particle that is + contained entirely within cubic domain with uniform probability. + + Returns + ------- + list of float + Cartesian coordinates of particle center. + + """ + + return [uniform(llim[0], ulim[0]), + uniform(llim[0], ulim[0]), + uniform(llim[0], ulim[0])] + + + def random_point_cylinder(): + """Generate Cartesian coordinates of center of a particle that is + contained entirely within cylindrical domain with uniform probability + (see http://mathworld.wolfram.com/DiskPointPicking.html for generating + random points on a disk). + + Returns + ------- + list of float + Cartesian coordinates of particle center. + + """ + + r = uniform(0, ulim[0]**2)**.5 + t = uniform(0, 2*pi) + return [r*cos(t), r*sin(t), uniform(llim[2], ulim[2])] + + + def random_point_sphere(): + """Generate Cartesian coordinates of center of a particle that is + contained entirely within spherical domain with uniform probability. + + Returns + ------- + list of float + Cartesian coordinates of particle center. + + """ + + x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) + r = (uniform(0, ulim[0]**3)**(1/3) / (x[0]**2 + x[1]**2 + x[2]**2)**.5) + return [r*i for i in x] + + + 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] = None + rod[2] = None + + + 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 is not None and j is not None: + 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.dstack(([i for i in range(len(n))], n, d))[0] + + # 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 = [x for x in {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 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)). + + """ + + j = floor(-log10(outer_packing_fraction() - inner_packing_fraction())) + outer_diameter[0] = (outer_diameter[0] - 0.5**j * + initial_outer_diameter * contraction_rate / + n_particles) + + + def update_mesh(i): + """Update which lattice 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 lattice cell and + which lattice 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 lattice 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 lattice cells are within one diameter of particle's + # center and add this particle to the list of particles in those cells + for idx in cell_list(particles[i], diameter): + mesh[idx].add(i) + mesh_map[i].add(idx) + + + def apply_boundary_conditions(i, j): + """Apply reflective boundary conditions to particles i and j. + + Parameters + ---------- + i, j : int + Index of particles in particles array. + + """ + + for k in range(3): + if particles[i][k] < llim[k]: + particles[i][k] = llim[k] + elif particles[i][k] > ulim[k]: + particles[i][k] = ulim[k] + if particles[j][k] < llim[k]: + particles[j][k] = llim[k] + elif particles[j][k] > ulim[k]: + particles[j][k] = ulim[k] + + + 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] = particles[i] + r*v + particles[j] = particles[j] - r*v + + # Apply reflective boundary conditions + apply_boundary_conditions(i, j) + + 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 + double + 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[cell_index(particles[i])]) + dists = 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] + + + def cell_index_cube(p, cl=None): + """Calculate the index of the lattice cell in which the given particle + center falls. + + Parameters + ---------- + p : list of float + Cartesian coordinates of particle center. + cl : list of float + Length of the lattice cells in x-, y-, and z-directions. + + Returns + ------- + tuple of int + Indices of lattice cell. + + """ + + if cl is None: + cl = cell_length + + return tuple(int(p[i]/cl[i]) for i in range(3)) + + + def cell_index_cylinder(p, cl=None): + """Calculate the index of the lattice cell in which the given particle + center falls. + + Parameters + ---------- + p : list of float + Cartesian coordinates of particle center. + cl : list of float + Length of the lattice cells in x-, y-, and z-directions. + + Returns + ------- + tuple of int + Indices of lattice cell. + + """ + + if cl is None: + cl = cell_length + + return tuple([int((p[0] + domain_radius)/cl[0]), + int((p[1] + domain_radius)/cl[1]), int(p[2]/cl[2])]) + + + def cell_index_sphere(p, cl=None): + """Calculate the index of the lattice cell in which the given particle + center falls. + + Parameters + ---------- + p : list of float + Cartesian coordinates of particle center. + cl : list of float + Length of the lattice cells in x-, y-, and z-directions. + + Returns + ------- + tuple of int + Indices of lattice cell. + + """ + + if cl is None: + cl = cell_length + + return tuple(int((p[i] + domain_radius)/cl[i]) for i in range(3)) + + + def cell_list_cube(p, d, cl=None): + """Return the indices of all cells within the given distance of the + point. + + Parameters + ---------- + p : list of float + Cartesian coordinates of particle center. + d : float + Find all lattice cells that are within a radius of length 'd' of + the particle center. + cl : list of float + Length of the lattice cells in x-, y-, and z-directions. + + Returns + ------- + list of tuple of int + Indices of lattice cells. + + """ + + if cl is None: + cl = cell_length + + r = [[a/cl[i] for a in [p[i]-d, p[i], p[i]+d] if a > 0 and + a < domain_length] for i in range(3)] + + return list(itertools.product(*({int(i) for i in j} for j in r))) + + + def cell_list_cylinder(p, d, cl=None): + """Return the indices of all cells within the given distance of the + point. + + Parameters + ---------- + p : list of float + Cartesian coordinates of particle center. + d : float + Find all lattice cells that are within a radius of length 'd' of + the particle center. + cl : list of float + Length of the lattice cells in x-, y-, and z-directions. + + Returns + ------- + list of tuple of int + Indices of lattice cells. + + """ + + if cl is None: + cl = cell_length + + x,y = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] + if a > -domain_radius and a < domain_radius] for i in range(2)] + + z = [a/cl[2] for a in [p[2]-d, p[2], p[2]+d] if a > 0 + and a < domain_length] + + return list(itertools.product(*({int(i) for i in j} for j in (x,y,z)))) + + + def cell_list_sphere(p, d, cl=None): + """Return the indices of all cells within the given distance of the + point. + + Parameters + ---------- + p : list of float + Cartesian coordinates of particle center. + d : float + Find all lattice cells that are within a radius of length 'd' of + the particle center. + cl : list of float + Length of the lattice cells in x-, y-, and z-directions. + + Returns + ------- + list of tuple of int + Indices of lattice cells. + + """ + + if cl is None: + cl = cell_length + + r = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] + if a > -domain_radius and a < domain_radius] for i in range(3)] + + return list(itertools.product(*({int(i) for i in j} for j in r))) + + + def random_sequential_pack(): + """Random sequential packing of particles whose radius is determined by + initial packing fraction. + + Returns + ------ + numpy.ndarray + Cartesian coordinates of centers of TRISO particles. + + """ + + # Set parameters for initial random sequential packing of particles. + r = (3/4*initial_packing_fraction*domain_volume/(pi*n_particles))**(1/3) + d = 2*r + sqd = d**2 + cl = get_cell_length(r) + + particles = [] + mesh = defaultdict(list) + + for i in range(n_particles): + # Randomly sample new center coordinates while there are any overlaps + while True: + p = random_point() + idx = cell_index(p, cl) + 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 cell_list(p, d, cl): + mesh[idx].append(p) + + return np.array(particles) + + + def close_random_pack(): + """Close random packing of particles using the Jodrey-Tory algorithm. + + """ + + for i in range(n_particles): + for idx in cell_list(particles[i], diameter): + 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 + + + # 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)) + + domain_volume = get_domain_volume() + llim, ulim = get_boundary_extremes() + offset = get_particle_offset() + + # 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: + packing_fraction = 4/3*pi*radius**3*n_particles / domain_volume + elif n_particles is None: + 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 + + # Set domain dependent functions + if domain_shape is 'cube': + random_point = random_point_cube + cell_list = cell_list_cube + cell_index = cell_index_cube + elif domain_shape is 'cylinder': + random_point = random_point_cylinder + cell_list = cell_list_cylinder + cell_index = cell_index_cylinder + elif domain_shape is 'sphere': + random_point = random_point_sphere + cell_list = cell_list_sphere + cell_index = cell_index_sphere + + random.seed(seed) + + # Generate non-overlapping particles for an initial inner radius using + # random sequential packing algorithm + particles = random_sequential_pack() + + # 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: + diameter = 2*radius + cell_length = get_cell_length(radius) + + # 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) + + close_random_pack() + + trisos = [] + for i in range(n_particles): + trisos.append(TRISO(radius, fill, particles[i] + offset)) + + return trisos diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index 04119a2dc..0df6042f6 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -1 +1 @@ -f33e6653b883200457df2ff2ba9cf715d5ddaa1296dd71d277c6f1d9d5b7831cc92aaf1e97509d26e5a93235cd9f775c0cfaa5ebc3dfe8fc71469bac166d362b \ No newline at end of file +2285ba99573743929cee590e2ba4d86becbf38d58af765498b73425f2fa3ccc3e0d20a59260d283a3ff39038d87e12142e0156b22152ccecbe2609a291d1d347 \ No newline at end of file diff --git a/tests/test_triso/results_true.dat b/tests/test_triso/results_true.dat index ea7da21ed..15107e8c8 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 9a8fb0c3b..8b6296ab5 100644 --- a/tests/test_triso/test_triso.py +++ b/tests/test_triso/test_triso.py @@ -19,35 +19,35 @@ class TRISOTestHarness(PyAPITestHarness): # Define TRISO matrials fuel = openmc.Material() fuel.set_density('g/cm3', 10.5) - fuel.add_nuclide('U235', 0.14154) - fuel.add_nuclide('U238', 0.85846) - fuel.add_nuclide('C0', 0.5) - fuel.add_nuclide('O16', 1.5) + fuel.add_nuclide('U-235', 0.14154) + fuel.add_nuclide('U-238', 0.85846) + fuel.add_nuclide('C-Nat', 0.5) + fuel.add_nuclide('O-16', 1.5) porous_carbon = openmc.Material() porous_carbon.set_density('g/cm3', 1.0) - porous_carbon.add_nuclide('C0', 1.0) - porous_carbon.add_s_alpha_beta('c_Graphite', '71t') + porous_carbon.add_nuclide('C-Nat', 1.0) + porous_carbon.add_s_alpha_beta('Graph', '71t') ipyc = openmc.Material() ipyc.set_density('g/cm3', 1.90) - ipyc.add_nuclide('C0', 1.0) - ipyc.add_s_alpha_beta('c_Graphite', '71t') + ipyc.add_nuclide('C-Nat', 1.0) + ipyc.add_s_alpha_beta('Graph', '71t') sic = openmc.Material() sic.set_density('g/cm3', 3.20) sic.add_element('Si', 1.0) - sic.add_nuclide('C0', 1.0) + sic.add_nuclide('C-Nat', 1.0) opyc = openmc.Material() opyc.set_density('g/cm3', 1.87) - opyc.add_nuclide('C0', 1.0) - opyc.add_s_alpha_beta('c_Graphite', '71t') + opyc.add_nuclide('C-Nat', 1.0) + opyc.add_s_alpha_beta('Graph', '71t') graphite = openmc.Material() graphite.set_density('g/cm3', 1.1995) - graphite.add_nuclide('C0', 1.0) - graphite.add_s_alpha_beta('c_Graphite', '71t') + graphite.add_nuclide('C-Nat', 1.0) + graphite.add_s_alpha_beta('Graph', '71t') # Create TRISO particles spheres = [openmc.Sphere(R=r*1e-4) @@ -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') From 15aaa15530904b83afd0f7c6444d10f98eca7196 Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 16 Aug 2016 20:25:07 -0500 Subject: [PATCH 02/11] Updated trisos test --- tests/test_triso/test_triso.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py index 8b6296ab5..da5deef00 100644 --- a/tests/test_triso/test_triso.py +++ b/tests/test_triso/test_triso.py @@ -19,35 +19,35 @@ class TRISOTestHarness(PyAPITestHarness): # Define TRISO matrials fuel = openmc.Material() fuel.set_density('g/cm3', 10.5) - fuel.add_nuclide('U-235', 0.14154) - fuel.add_nuclide('U-238', 0.85846) - fuel.add_nuclide('C-Nat', 0.5) - fuel.add_nuclide('O-16', 1.5) + fuel.add_nuclide('U235', 0.14154) + fuel.add_nuclide('U238', 0.85846) + fuel.add_nuclide('C0', 0.5) + fuel.add_nuclide('O16', 1.5) porous_carbon = openmc.Material() porous_carbon.set_density('g/cm3', 1.0) - porous_carbon.add_nuclide('C-Nat', 1.0) - porous_carbon.add_s_alpha_beta('Graph', '71t') + porous_carbon.add_nuclide('C0', 1.0) + porous_carbon.add_s_alpha_beta('c_Graphite', '71t') ipyc = openmc.Material() ipyc.set_density('g/cm3', 1.90) - ipyc.add_nuclide('C-Nat', 1.0) - ipyc.add_s_alpha_beta('Graph', '71t') + ipyc.add_nuclide('C0', 1.0) + ipyc.add_s_alpha_beta('c_Graphite', '71t') sic = openmc.Material() sic.set_density('g/cm3', 3.20) sic.add_element('Si', 1.0) - sic.add_nuclide('C-Nat', 1.0) + sic.add_nuclide('C0', 1.0) opyc = openmc.Material() opyc.set_density('g/cm3', 1.87) - opyc.add_nuclide('C-Nat', 1.0) - opyc.add_s_alpha_beta('Graph', '71t') + opyc.add_nuclide('C0', 1.0) + opyc.add_s_alpha_beta('c_Graphite', '71t') graphite = openmc.Material() graphite.set_density('g/cm3', 1.1995) - graphite.add_nuclide('C-Nat', 1.0) - graphite.add_s_alpha_beta('Graph', '71t') + graphite.add_nuclide('C0', 1.0) + graphite.add_s_alpha_beta('c_Graphite', '71t') # Create TRISO particles spheres = [openmc.Sphere(R=r*1e-4) From 760168a82e1e3ffb4e4024ededdd84249a901f51 Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 17 Aug 2016 16:21:10 -0500 Subject: [PATCH 03/11] Address #706 comments --- docs/source/pythonapi/index.rst | 1 + openmc/model/triso.py | 42 ++++++++++++++++----------------- 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 392c2c72d..14f4a2128 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/openmc/model/triso.py b/openmc/model/triso.py index c1ed9336a..763fc0ddd 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -4,14 +4,14 @@ from collections import Iterable, defaultdict from numbers import Real import warnings import itertools -import scipy.spatial -from scipy.spatial.distance import cdist import random from random import uniform, gauss from heapq import heappush, heappop from math import pi, sin, cos, floor, log10 import numpy as np +import scipy.spatial +from scipy.spatial.distance import cdist import openmc import openmc.checkvalue as cv @@ -433,8 +433,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ rod = [d, i, j] - rods_map[i] = j, rod - rods_map[j] = i, rod + rods_map[i] = (j, rod) + rods_map[j] = (i, rod) heappush(rods, rod) @@ -516,7 +516,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, # 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])) + for x in r])) # Clear priority queue and add rods del rods[:] @@ -617,7 +617,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, # 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; + r = (outer_diameter[0] - d)/2 v = (particles[i] - particles[j])/d particles[i] = particles[i] + r*v @@ -735,7 +735,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, cl = cell_length return tuple([int((p[0] + domain_radius)/cl[0]), - int((p[1] + domain_radius)/cl[1]), int(p[2]/cl[2])]) + int((p[1] + domain_radius)/cl[1]), int(p[2]/cl[2])]) def cell_index_sphere(p, cl=None): @@ -787,7 +787,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, cl = cell_length r = [[a/cl[i] for a in [p[i]-d, p[i], p[i]+d] if a > 0 and - a < domain_length] for i in range(3)] + a < domain_length] for i in range(3)] return list(itertools.product(*({int(i) for i in j} for j in r))) @@ -816,13 +816,13 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, if cl is None: cl = cell_length - x,y = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] - if a > -domain_radius and a < domain_radius] for i in range(2)] + x, y = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] + if a > -domain_radius and a < domain_radius] for i in range(2)] z = [a/cl[2] for a in [p[2]-d, p[2], p[2]+d] if a > 0 and a < domain_length] - return list(itertools.product(*({int(i) for i in j} for j in (x,y,z)))) + return list(itertools.product(*({int(i) for i in j} for j in (x, y, z)))) def cell_list_sphere(p, d, cl=None): @@ -850,7 +850,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, cl = cell_length r = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] - if a > -domain_radius and a < domain_radius] for i in range(3)] + if a > -domain_radius and a < domain_radius] for i in range(3)] return list(itertools.product(*({int(i) for i in j} for j in r))) @@ -878,13 +878,13 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, for i in range(n_particles): # Randomly sample new center coordinates while there are any overlaps while True: - p = random_point() - idx = cell_index(p, cl) - 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 + p = random_point() + idx = cell_index(p, cl) + 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 cell_list(p, d, cl): @@ -998,7 +998,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, close_random_pack() trisos = [] - for i in range(n_particles): - trisos.append(TRISO(radius, fill, particles[i] + offset)) + for p in particles: + trisos.append(TRISO(radius, fill, p + offset)) return trisos From a13d5a27480cb4d38a0c1774a07fd8b5696d530c Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 18 Aug 2016 10:54:02 -0500 Subject: [PATCH 04/11] Address #706 comments --- openmc/model/triso.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 763fc0ddd..33c955e3b 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -7,7 +7,7 @@ import itertools import random from random import uniform, gauss from heapq import heappush, heappop -from math import pi, sin, cos, floor, log10 +from math import pi, sin, cos, floor, log10, sqrt import numpy as np import scipy.spatial @@ -224,8 +224,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, 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 rondom, and placement - attempts for a particles are made until the particle is not overlapping any + 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 lattice over the domain to speed up the nearest neighbor search by only searching for a particle's neighbors within that lattice cell. @@ -333,7 +333,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ if domain_shape is 'cube': - return np.array(domain_center) - 3*(domain_length/2,) + return np.array(domain_center) - domain_length/2 elif domain_shape is 'cylinder': return np.array(domain_center) - (0, 0, domain_length/2) elif domain_shape is 'sphere': @@ -399,7 +399,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ - r = uniform(0, ulim[0]**2)**.5 + r = sqrt(uniform(0, ulim[0]**2)) t = uniform(0, 2*pi) return [r*cos(t), r*sin(t), uniform(llim[2], ulim[2])] @@ -416,7 +416,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) - r = (uniform(0, ulim[0]**3)**(1/3) / (x[0]**2 + x[1]**2 + x[2]**2)**.5) + r = (uniform(0, ulim[0]**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2)) return [r*i for i in x] @@ -620,8 +620,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, r = (outer_diameter[0] - d)/2 v = (particles[i] - particles[j])/d - particles[i] = particles[i] + r*v - particles[j] = particles[j] - r*v + particles[i] += r*v + particles[j] -= r*v # Apply reflective boundary conditions apply_boundary_conditions(i, j) @@ -643,7 +643,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, ------- int Index in particles array of nearest neighbor of i - double + float distance between i and nearest neighbor. """ @@ -734,8 +734,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, if cl is None: cl = cell_length - return tuple([int((p[0] + domain_radius)/cl[0]), - int((p[1] + domain_radius)/cl[1]), int(p[2]/cl[2])]) + return (int((p[0] + domain_radius)/cl[0]), + int((p[1] + domain_radius)/cl[1]), int(p[2]/cl[2])) def cell_index_sphere(p, cl=None): @@ -935,7 +935,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, # 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)): + (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: From 50042b4fe3a12b7a98359b236c073b482e2646f3 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 25 Aug 2016 22:39:40 -0500 Subject: [PATCH 05/11] Added domain classes; Moved packing algorithms into separate functions --- openmc/model/triso.py | 1121 +++++++++++++++--------------- tests/test_triso/inputs_true.dat | 2 +- 2 files changed, 570 insertions(+), 553 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 33c955e3b..d0d67b687 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -1,13 +1,14 @@ from __future__ import division import copy -from collections import Iterable, defaultdict -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 import scipy.spatial @@ -91,6 +92,363 @@ 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._particle_radius = None + self._center = None + 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 + + @property + def cell_length(self): + return self._cell_length + + @property + def limits(self): + return self._limits + + @abstractproperty + def volume(self): + pass + + @particle_radius.setter + def particle_radius(self, particle_radius): + self._particle_radius = float(particle_radius) + self.reset() + + @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.reset() + + @cell_length.setter + def cell_length(self, cell_length): + self._cell_length = cell_length + + @limits.setter + def limits(self, limits): + self._limits = limits + + 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 reset(self): + """Recalculate attributes that depend on input parameters if any of the + parameters are modified. + + """ + pass + + @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.]): + self._length = None + super(_CubicDomain, self).__init__(particle_radius, center) + self.length = length + + @property + def volume(self): + return self.length**3 + + @property + def length(self): + return self._length + + @length.setter + def length(self, length): + self._length = float(length) + self.reset() + + def reset(self): + if (self.particle_radius is not None and self.center is not None + and self.length is not None): + xlim = self.length/2 - self.particle_radius + self.limits = [[x - xlim for x in self.center], + [x + xlim for x in self.center]] + mesh_length = [self.length, self.length, self.length] + self.cell_length = [x/int(x/(4*self.particle_radius)) + for x in mesh_length] + + 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.]): + self._length = None + self._radius = None + super(_CylindricalDomain, self).__init__(particle_radius, center) + self.length = length + self.radius = radius + + @property + def volume(self): + return self.length * pi * self.radius**2 + + @property + def length(self): + return self._length + + @property + def radius(self): + return self._radius + + @length.setter + def length(self, length): + self._length = float(length) + self.reset() + + @radius.setter + def radius(self, radius): + self._radius = float(radius) + self.reset() + + def reset(self): + if (self.particle_radius is not None and self.center is not None + and self.length is not None and self.radius is not 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]] + 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] + + 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.]): + self._radius = None + super(_SphericalDomain, self).__init__(particle_radius, center) + self.radius = radius + + @property + def volume(self): + return 4/3 * pi * self.radius**3 + + @property + def radius(self): + return self._radius + + @radius.setter + def radius(self, radius): + self._radius = float(radius) + self.reset() + + def reset(self): + if (self.particle_radius is not None and self.center is not None + and self.radius is not None): + rlim = self.radius - self.particle_radius + self.limits = [[x - rlim for x in self.center], + [x + rlim for x in self.center]] + 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] + + 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. @@ -164,261 +522,58 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): return lattice -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. +def _random_sequential_pack(domain, n_particles): + """Random sequential packing of 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. + domain : openmc.model._Domain + Container in which to pack particles. 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. + Number of particles to pack. 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 is faster 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 performs better - 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 lattice over the domain - to speed up the nearest neighbor search by only searching for a particle's - neighbors within that lattice 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 and the outer diameter - is decreased. 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. + ------ + numpy.ndarray + Cartesian coordinates of centers of particles. """ - def get_domain_volume(): - """Calculates the volume of the container in which the TRISO particles - are packed. + sqd = (2*domain.particle_radius)**2 + particles = [] + mesh = defaultdict(list) - Returns - ------- - float - Volume of the domain. + 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) - if domain_shape is 'cube': - return domain_length**3 - elif domain_shape is 'cylinder': - return domain_length * pi * domain_radius**2 - elif domain_shape is 'sphere': - return 4/3 * pi * domain_radius**3 + return np.array(particles) - def get_cell_length(radius): - """Calculates the length of a lattice element in x-, y-, and - z-directions. +def _close_random_pack(domain, particles, contraction_rate): + """Close random packing of particles using the Jodrey-Tory algorithm. - Parameters - ---------- - radius : float - Radius of the particle. - - Returns - ------- - tuple of float - Length of lattice cell in x-, y-, and z-directions. - - """ - - if domain_length: - m = domain_length/int(domain_length/(4*radius)) - if domain_radius: - n = 2*domain_radius/int(domain_radius/(2*radius)) - - if domain_shape is 'cube': - return (m, m, m) - elif domain_shape is 'cylinder': - return (n, n, m) - elif domain_shape is 'sphere': - return (n, n, n) - - - def get_boundary_extremes(): - """Calculates the minimum and maximum positions in x-, y-, and - z-directions where a particle center can be placed within the domain. - - Returns - ------- - llim, ulim : tuple of float - Minimum and maximum position in x-, y-, and z-directions where - particle center can be placed. - - """ - - if domain_length: - x_min = radius - x_max = domain_length - radius - if domain_radius: - r_min = radius - domain_radius - r_max = domain_radius - radius - - if domain_shape is 'cube': - return (x_min, x_min, x_min), (x_max, x_max, x_max) - elif domain_shape is 'cylinder': - return (r_min, r_min, x_min), (r_max, r_max, x_max) - elif domain_shape is 'sphere': - return (r_min, r_min, r_min), (r_max, r_max, r_max) - - - def get_particle_offset(): - """Calculates the offset in x-, y-, and z-directions of the particle - center based on the domain center - - Returns - ------- - tuple of float - Amount to offset particle center in x-, y-, and z-directions - - """ - - if domain_shape is 'cube': - return np.array(domain_center) - domain_length/2 - elif domain_shape is 'cylinder': - return np.array(domain_center) - (0, 0, domain_length/2) - elif domain_shape is 'sphere': - return np.array(domain_center) - - - def inner_packing_fraction(): - """Calculates the true packing fraction of the particles based on the - inner diameter. - - Returns - ------- - float - Packing fraction calculated from inner diameter. - - """ - - return (4/3 * pi * (inner_diameter[0]/2)**3 * n_particles / - domain_volume) - - - def outer_packing_fraction(): - """Calculates the nominal packing fraction of the particles based on - the outer diameter. - - Returns - ------- - float - Packing fraction calculated from outer diameter. - - """ - - return (4/3 * pi * (outer_diameter[0]/2)**3 * n_particles / - domain_volume) - - - def random_point_cube(): - """Generate Cartesian coordinates of center of a particle that is - contained entirely within cubic domain with uniform probability. - - Returns - ------- - list of float - Cartesian coordinates of particle center. - - """ - - return [uniform(llim[0], ulim[0]), - uniform(llim[0], ulim[0]), - uniform(llim[0], ulim[0])] - - - def random_point_cylinder(): - """Generate Cartesian coordinates of center of a particle that is - contained entirely within cylindrical domain with uniform probability - (see http://mathworld.wolfram.com/DiskPointPicking.html for generating - random points on a disk). - - Returns - ------- - list of float - Cartesian coordinates of particle center. - - """ - - r = sqrt(uniform(0, ulim[0]**2)) - t = uniform(0, 2*pi) - return [r*cos(t), r*sin(t), uniform(llim[2], ulim[2])] - - - def random_point_sphere(): - """Generate Cartesian coordinates of center of a particle that is - contained entirely within spherical domain with uniform probability. - - Returns - ------- - list of float - Cartesian coordinates of particle center. - - """ - - x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) - r = (uniform(0, ulim[0]**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2)) - return [r*i for i in x] + 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. @@ -530,6 +685,35 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, 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: @@ -541,60 +725,14 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ - j = floor(-log10(outer_packing_fraction() - inner_packing_fraction())) - outer_diameter[0] = (outer_diameter[0] - 0.5**j * - initial_outer_diameter * contraction_rate / - n_particles) + 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) - - def update_mesh(i): - """Update which lattice 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 lattice cell and - which lattice 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 lattice 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 lattice cells are within one diameter of particle's - # center and add this particle to the list of particles in those cells - for idx in cell_list(particles[i], diameter): - mesh[idx].add(i) - mesh_map[i].add(idx) - - - def apply_boundary_conditions(i, j): - """Apply reflective boundary conditions to particles i and j. - - Parameters - ---------- - i, j : int - Index of particles in particles array. - - """ - - for k in range(3): - if particles[i][k] < llim[k]: - particles[i][k] = llim[k] - elif particles[i][k] > ulim[k]: - particles[i][k] = ulim[k] - if particles[j][k] < llim[k]: - particles[j][k] = llim[k] - elif particles[j][k] > ulim[k]: - particles[j][k] = ulim[k] + 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): @@ -624,7 +762,15 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, particles[j] -= r*v # Apply reflective boundary conditions - apply_boundary_conditions(i, j) + for k in range(3): + if particles[i][k] < domain.limits[0][k]: + particles[i][k] = domain.limits[0][k] + elif particles[i][k] > domain.limits[1][k]: + particles[i][k] = domain.limits[1][k] + if particles[j][k] < domain.limits[0][k]: + particles[j][k] = domain.limits[0][k] + elif particles[j][k] > domain.limits[1][k]: + particles[j][k] = domain.limits[1][k] update_mesh(i) update_mesh(j) @@ -651,7 +797,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, # 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[cell_index(particles[i])]) + idx = list(mesh[domain.mesh_cell(particles[i])]) dists = cdist([particles[i]], particles[idx])[0] if dists.size > 1: j = dists.argpartition(1)[1] @@ -689,233 +835,120 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, inner_diameter[0] = rods[0][0] - def cell_index_cube(p, cl=None): - """Calculate the index of the lattice cell in which the given particle - center falls. + n_particles = len(particles) + diameter = 2*domain.particle_radius - Parameters - ---------- - p : list of float - Cartesian coordinates of particle center. - cl : list of float - Length of the lattice cells in x-, y-, and z-directions. + # 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) - Returns - ------- - tuple of int - Indices of lattice cell. + # 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) - if cl is None: - cl = cell_length - - return tuple(int(p[i]/cl[i]) for i in range(3)) - - - def cell_index_cylinder(p, cl=None): - """Calculate the index of the lattice cell in which the given particle - center falls. - - Parameters - ---------- - p : list of float - Cartesian coordinates of particle center. - cl : list of float - Length of the lattice cells in x-, y-, and z-directions. - - Returns - ------- - tuple of int - Indices of lattice cell. - - """ - - if cl is None: - cl = cell_length - - return (int((p[0] + domain_radius)/cl[0]), - int((p[1] + domain_radius)/cl[1]), int(p[2]/cl[2])) - - - def cell_index_sphere(p, cl=None): - """Calculate the index of the lattice cell in which the given particle - center falls. - - Parameters - ---------- - p : list of float - Cartesian coordinates of particle center. - cl : list of float - Length of the lattice cells in x-, y-, and z-directions. - - Returns - ------- - tuple of int - Indices of lattice cell. - - """ - - if cl is None: - cl = cell_length - - return tuple(int((p[i] + domain_radius)/cl[i]) for i in range(3)) - - - def cell_list_cube(p, d, cl=None): - """Return the indices of all cells within the given distance of the - point. - - Parameters - ---------- - p : list of float - Cartesian coordinates of particle center. - d : float - Find all lattice cells that are within a radius of length 'd' of - the particle center. - cl : list of float - Length of the lattice cells in x-, y-, and z-directions. - - Returns - ------- - list of tuple of int - Indices of lattice cells. - - """ - - if cl is None: - cl = cell_length - - r = [[a/cl[i] for a in [p[i]-d, p[i], p[i]+d] if a > 0 and - a < domain_length] for i in range(3)] - - return list(itertools.product(*({int(i) for i in j} for j in r))) - - - def cell_list_cylinder(p, d, cl=None): - """Return the indices of all cells within the given distance of the - point. - - Parameters - ---------- - p : list of float - Cartesian coordinates of particle center. - d : float - Find all lattice cells that are within a radius of length 'd' of - the particle center. - cl : list of float - Length of the lattice cells in x-, y-, and z-directions. - - Returns - ------- - list of tuple of int - Indices of lattice cells. - - """ - - if cl is None: - cl = cell_length - - x, y = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] - if a > -domain_radius and a < domain_radius] for i in range(2)] - - z = [a/cl[2] for a in [p[2]-d, p[2], p[2]+d] if a > 0 - and a < domain_length] - - return list(itertools.product(*({int(i) for i in j} for j in (x, y, z)))) - - - def cell_list_sphere(p, d, cl=None): - """Return the indices of all cells within the given distance of the - point. - - Parameters - ---------- - p : list of float - Cartesian coordinates of particle center. - d : float - Find all lattice cells that are within a radius of length 'd' of - the particle center. - cl : list of float - Length of the lattice cells in x-, y-, and z-directions. - - Returns - ------- - list of tuple of int - Indices of lattice cells. - - """ - - if cl is None: - cl = cell_length - - r = [[(a + domain_radius)/cl[i] for a in [p[i]-d, p[i], p[i]+d] - if a > -domain_radius and a < domain_radius] for i in range(3)] - - return list(itertools.product(*({int(i) for i in j} for j in r))) - - - def random_sequential_pack(): - """Random sequential packing of particles whose radius is determined by - initial packing fraction. - - Returns - ------ - numpy.ndarray - Cartesian coordinates of centers of TRISO particles. - - """ - - # Set parameters for initial random sequential packing of particles. - r = (3/4*initial_packing_fraction*domain_volume/(pi*n_particles))**(1/3) - d = 2*r - sqd = d**2 - cl = get_cell_length(r) - - particles = [] - mesh = defaultdict(list) - - for i in range(n_particles): - # Randomly sample new center coordinates while there are any overlaps - while True: - p = random_point() - idx = cell_index(p, cl) - 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 cell_list(p, d, cl): - mesh[idx].append(p) - - return np.array(particles) - - - def close_random_pack(): - """Close random packing of particles using the Jodrey-Tory algorithm. - - """ - - for i in range(n_particles): - for idx in cell_list(particles[i], diameter): - mesh[idx].add(i) - mesh_map[i].add(idx) + 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: - create_rod_list() - if inner_diameter[0] >= diameter: + 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 - 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", ' @@ -928,9 +961,15 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, raise ValueError('"domain_radius" must be specified for {} domain ' 'geometry '.format(domain_shape)) - domain_volume = get_domain_volume() - llim, ulim = get_boundary_extremes() - offset = get_particle_offset() + 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. @@ -939,9 +978,9 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, raise ValueError('Exactly one of "n_particles" and "packing_fraction" ' 'must be specified.') elif packing_fraction is None: - packing_fraction = 4/3*pi*radius**3*n_particles / domain_volume + packing_fraction = 4/3*pi*radius**3*n_particles / domain.volume elif n_particles is None: - n_particles = int(packing_fraction*domain_volume // (4/3*pi*radius**3)) + 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: @@ -957,48 +996,26 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, if packing_fraction > 0.3: initial_packing_fraction = 0.3 - # Set domain dependent functions - if domain_shape is 'cube': - random_point = random_point_cube - cell_list = cell_list_cube - cell_index = cell_index_cube - elif domain_shape is 'cylinder': - random_point = random_point_cylinder - cell_list = cell_list_cylinder - cell_index = cell_index_cylinder - elif domain_shape is 'sphere': - random_point = random_point_sphere - cell_list = cell_list_sphere - cell_index = cell_index_sphere - random.seed(seed) + # Set parameters for initial random sequential packing of particles. + initial_radius = (3/4 * initial_packing_fraction * domain.volume / + (pi * n_particles))**(1/3) + domain.particle_radius = initial_radius + 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() + 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: - diameter = 2*radius - cell_length = get_cell_length(radius) - - # 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) - - close_random_pack() + domain.particle_radius = radius + _close_random_pack(domain, particles, contraction_rate) trisos = [] for p in particles: - trisos.append(TRISO(radius, fill, p + offset)) - + trisos.append(TRISO(radius, fill, p)) return trisos diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index 0df6042f6..aed053263 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -1 +1 @@ -2285ba99573743929cee590e2ba4d86becbf38d58af765498b73425f2fa3ccc3e0d20a59260d283a3ff39038d87e12142e0156b22152ccecbe2609a291d1d347 \ No newline at end of file +792d82b08d5fa6ac19df668d84cb90cbbe178e8769c87b732e9616ffeb70d8cfb4096607e58cda749cb800bb48139ac72f6983e84f38237d7d55711bdd1c6d4d \ No newline at end of file From 6d82c7b48b9a83f0b6ad7f73c57562e920954e97 Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 26 Aug 2016 12:08:13 -0500 Subject: [PATCH 06/11] Address #706 comments --- openmc/model/triso.py | 190 +++++++++++++++++++++--------------------- 1 file changed, 97 insertions(+), 93 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index d0d67b687..83be99a33 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -124,8 +124,6 @@ class _Domain(object): __metaclass__ = ABCMeta def __init__(self, particle_radius, center=[0., 0., 0.]): - self._particle_radius = None - self._center = None self._cell_length = None self._limits = None @@ -140,13 +138,13 @@ class _Domain(object): def center(self): return self._center - @property - def cell_length(self): - return self._cell_length - - @property + @abstractproperty def limits(self): - return self._limits + pass + + @abstractproperty + def cell_length(self): + pass @abstractproperty def volume(self): @@ -155,7 +153,8 @@ class _Domain(object): @particle_radius.setter def particle_radius(self, particle_radius): self._particle_radius = float(particle_radius) - self.reset() + self._limits = None + self._cell_length = None @center.setter def center(self, center): @@ -163,15 +162,8 @@ class _Domain(object): 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.reset() - - @cell_length.setter - def cell_length(self, cell_length): - self._cell_length = cell_length - - @limits.setter - def limits(self, limits): - self._limits = limits + 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 @@ -210,14 +202,6 @@ class _Domain(object): for i in range(3)] return list(itertools.product(*({int(x) for x in y} for y in r))) - @abstractmethod - def reset(self): - """Recalculate attributes that depend on input parameters if any of the - parameters are modified. - - """ - pass - @abstractmethod def random_point(self): """Generate Cartesian coordinates of center of a particle that is @@ -270,28 +254,39 @@ class _CubicDomain(_Domain): super(_CubicDomain, self).__init__(particle_radius, center) self.length = length - @property - def volume(self): - return self.length**3 - @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.reset() + self._limits = None + self._cell_length = None - def reset(self): - if (self.particle_radius is not None and self.center is not None - and self.length is not None): - xlim = self.length/2 - self.particle_radius - self.limits = [[x - xlim for x in self.center], - [x + xlim for x in self.center]] - mesh_length = [self.length, self.length, self.length] - self.cell_length = [x/int(x/(4*self.particle_radius)) - for x in mesh_length] + @limits.setter + def limits(self, limits): + self._limits = limits def random_point(self): return [uniform(self.limits[0][0], self.limits[1][0]), @@ -341,10 +336,6 @@ class _CylindricalDomain(_Domain): self.length = length self.radius = radius - @property - def volume(self): - return self.length * pi * self.radius**2 - @property def length(self): return self._length @@ -353,28 +344,44 @@ class _CylindricalDomain(_Domain): 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.reset() + self._limits = None + self._cell_length = None @radius.setter def radius(self, radius): self._radius = float(radius) - self.reset() + self._limits = None + self._cell_length = None - def reset(self): - if (self.particle_radius is not None and self.center is not None - and self.length is not None and self.radius is not 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]] - 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] + @limits.setter + def limits(self, limits): + self._limits = limits def random_point(self): r = sqrt(uniform(0, (self.radius - self.particle_radius)**2)) @@ -419,28 +426,39 @@ class _SphericalDomain(_Domain): super(_SphericalDomain, self).__init__(particle_radius, center) self.radius = radius - @property - def volume(self): - return 4/3 * pi * self.radius**3 - @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.reset() + self._limits = None + self._cell_length = None - def reset(self): - if (self.particle_radius is not None and self.center is not None - and self.radius is not None): - rlim = self.radius - self.particle_radius - self.limits = [[x - rlim for x in self.center], - [x + rlim for x in self.center]] - 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] + @limits.setter + def limits(self, limits): + self._limits = limits def random_point(self): x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) @@ -592,7 +610,6 @@ def _close_random_pack(domain, particles, contraction_rate): rods_map[j] = (i, rod) heappush(rods, rod) - def remove_rod(i): """Mark the rod containing particle i as removed. @@ -609,7 +626,6 @@ def _close_random_pack(domain, particles, contraction_rate): rod[1] = None rod[2] = None - def pop_rod(): """Remove and return the shortest rod. @@ -629,7 +645,6 @@ def _close_random_pack(domain, particles, contraction_rate): del rods_map[j] return d, i, j - def create_rod_list(): """Generate sorted list of rods (distances between particle centers). @@ -660,14 +675,15 @@ def _close_random_pack(domain, particles, contraction_rate): # distances to nearest neighbors a = np.dstack(([i for i in range(len(n))], n, d))[0] - # Array of nearest neighbor indices, indices of particles they are + # 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 = [x for x in {tuple(x) for x in a} & {tuple(x) for x in b}] + r = list([x for x in {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]))) @@ -684,7 +700,6 @@ def _close_random_pack(domain, particles, contraction_rate): 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. @@ -713,7 +728,6 @@ def _close_random_pack(domain, particles, contraction_rate): 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: @@ -725,9 +739,9 @@ def _close_random_pack(domain, particles, contraction_rate): """ - inner_pf = (4/3 * pi * (inner_diameter[0]/2)**3 * n_particles / + 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 / + outer_pf = (4/3 * pi * (outer_diameter[0]/2)**3 * n_particles / domain.volume) j = floor(-log10(outer_pf - inner_pf)) @@ -762,20 +776,12 @@ def _close_random_pack(domain, particles, contraction_rate): particles[j] -= r*v # Apply reflective boundary conditions - for k in range(3): - if particles[i][k] < domain.limits[0][k]: - particles[i][k] = domain.limits[0][k] - elif particles[i][k] > domain.limits[1][k]: - particles[i][k] = domain.limits[1][k] - if particles[j][k] < domain.limits[0][k]: - particles[j][k] = domain.limits[0][k] - elif particles[j][k] > domain.limits[1][k]: - particles[j][k] = domain.limits[1][k] + 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. @@ -805,7 +811,6 @@ def _close_random_pack(domain, particles, contraction_rate): 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. @@ -834,7 +839,6 @@ def _close_random_pack(domain, particles, contraction_rate): if rods: inner_diameter[0] = rods[0][0] - n_particles = len(particles) diameter = 2*domain.particle_radius From 19eca29642bf0cbe7d4bc024e9fc93ecce3fe71e Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 26 Aug 2016 16:19:12 -0500 Subject: [PATCH 07/11] Address #706 comments --- openmc/model/triso.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 83be99a33..b75ffd88f 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -673,7 +673,7 @@ def _close_random_pack(domain, particles, contraction_rate): # Array of particle indices, indices of nearest neighbors, and # distances to nearest neighbors - a = np.dstack(([i for i in range(len(n))], n, d))[0] + 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 From 56f85323b84d9b15e31fc9934bc217230723f843 Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 26 Aug 2016 17:13:09 -0500 Subject: [PATCH 08/11] Make sure input parameters are native Python types --- openmc/model/triso.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index b75ffd88f..e8ed04e3a 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -982,8 +982,10 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=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 From 813efa3b66dd9b89378621bd0057c2520287b152 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 29 Aug 2016 10:45:34 -0500 Subject: [PATCH 09/11] Address #706 comments --- openmc/model/triso.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index e8ed04e3a..76293051b 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -11,8 +11,11 @@ from math import pi, sin, cos, floor, log10, sqrt from abc import ABCMeta, abstractproperty, abstractmethod import numpy as np -import scipy.spatial -from scipy.spatial.distance import cdist +try: + import scipy.spatial + _SCIPY_AVAILABLE = True +except ImportError: + _SCIPY_AVAILABLE = False import openmc import openmc.checkvalue as cv @@ -250,7 +253,6 @@ class _CubicDomain(_Domain): """ def __init__(self, length, particle_radius, center=[0., 0., 0.]): - self._length = None super(_CubicDomain, self).__init__(particle_radius, center) self.length = length @@ -330,8 +332,6 @@ class _CylindricalDomain(_Domain): """ def __init__(self, length, radius, particle_radius, center=[0., 0., 0.]): - self._length = None - self._radius = None super(_CylindricalDomain, self).__init__(particle_radius, center) self.length = length self.radius = radius @@ -422,7 +422,6 @@ class _SphericalDomain(_Domain): """ def __init__(self, radius, particle_radius, center=[0., 0., 0.]): - self._radius = None super(_SphericalDomain, self).__init__(particle_radius, center) self.radius = radius @@ -435,7 +434,7 @@ class _SphericalDomain(_Domain): 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]] + [x + rlim for x in self.center]] return self._limits @property @@ -443,7 +442,7 @@ class _SphericalDomain(_Domain): 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] + for x in mesh_length] return self._cell_length @property @@ -683,7 +682,7 @@ def _close_random_pack(domain, particles, contraction_rate): # 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([x for x in {tuple(x) for x in a} & {tuple(x) for x in b}]) + 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]))) @@ -804,7 +803,7 @@ def _close_random_pack(domain, particles, contraction_rate): # will be itself. Using argpartition, the k-th nearest neighbor is # placed at index k. idx = list(mesh[domain.mesh_cell(particles[i])]) - dists = cdist([particles[i]], particles[idx])[0] + dists = scipy.spatial.distance.cdist([particles[i]], particles[idx])[0] if dists.size > 1: j = dists.argpartition(1)[1] return idx[j], dists[j] @@ -839,6 +838,10 @@ def _close_random_pack(domain, particles, contraction_rate): 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 @@ -1004,10 +1007,15 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, random.seed(seed) - # Set parameters for initial random sequential packing of particles. + # 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]]] @@ -1016,7 +1024,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, 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 + # 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) From ce6b893dd4889dc3de202a2bf96cb21e92fc2898 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 29 Aug 2016 14:24:52 -0500 Subject: [PATCH 10/11] Changed flag marking removed rods to integer to fix unorderable type error --- openmc/model/triso.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 76293051b..0c372aeef 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -622,8 +622,8 @@ def _close_random_pack(domain, particles, contraction_rate): if i in rods_map: j, rod = rods_map.pop(i) del rods_map[j] - rod[1] = None - rod[2] = None + rod[1] = removed + rod[2] = removed def pop_rod(): """Remove and return the shortest rod. @@ -639,7 +639,7 @@ def _close_random_pack(domain, particles, contraction_rate): while rods: d, i, j = heappop(rods) - if i is not None and j is not None: + if i is not removed and j is not removed: del rods_map[i] del rods_map[j] return d, i, j @@ -845,6 +845,9 @@ def _close_random_pack(domain, particles, contraction_rate): 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) From 4631c5c7df44395536cf2d5bef2f4f7976f6c404 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 29 Aug 2016 14:27:53 -0500 Subject: [PATCH 11/11] is not -> != --- openmc/model/triso.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 0c372aeef..5525ea559 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -639,7 +639,7 @@ def _close_random_pack(domain, particles, contraction_rate): while rods: d, i, j = heappop(rods) - if i is not removed and j is not removed: + if i != removed and j != removed: del rods_map[i] del rods_map[j] return d, i, j