From f3cfc181f47d959158d5e871209c45b9aa810f33 Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 16 Aug 2016 20:14:34 -0500 Subject: [PATCH 01/21] 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/21] 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/21] 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/21] 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 aaa80eb98926a271b59df6b992a90c48df88d637 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 23 Aug 2016 06:26:20 -0500 Subject: [PATCH 05/21] Update title() subroutine with ASCII logo --- src/output.F90 | 60 +++++++++++++++++++++++++++++++++----------------- 1 file changed, 40 insertions(+), 20 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 9e23f11d8..29950a0bc 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -38,43 +38,63 @@ contains use omp_lib #endif - write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & - ' .d88888b. 888b d888 .d8888b.', & - ' d88P" "Y88b 8888b d8888 d88P Y88b', & - ' 888 888 88888b.d88888 888 888', & - ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & - ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & - ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & - ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & - ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & - '__________________888______________________________________________________', & - ' 888', & - ' 888' + write(UNIT=OUTPUT_UNIT, FMT='(/29(A/))') & + ' %%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ################### %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ###################### %%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ####################### %%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ######################## %%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ######################### %%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ########################## %%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' ############################ %%%%%%%%%%%%%%%%%%%%%%%%', & + ' ############################# %%%%%%%%%%%%%%%%%%%%%%%', & + ' ############################## %%%%%%%%%%%%%%%%%%%%%', & + ' ############################# %%%%%%%%%%%%%%%%%%%%', & + ' ########################### %%%%%%%%%%%%%%%%%%%%', & + ' ######################### %%%%%%%%%%%%%%%%%%%%%', & + ' ###################### %%%%%%%%%%%%%%%%%%%%%', & + ' #################### %%%%%%%%%%%%%%%%%%%%%', & + ' ################# %%%%%%%%%%%%%%%%%%%%', & + ' ############## %%%%%%%%%%%%%%%%%%%', & + ' ########### %%%%%%%%%%%%%%%%%%', & + ' ####### %%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%' + ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Copyright: 2011-2016 Massachusetts Institute of Technology' + ' | The OpenMC Monte Carlo Code' write(UNIT=OUTPUT_UNIT, FMT=*) & - ' License: http://openmc.readthedocs.io/en/latest/license.html' - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') & + ' Copyright | 2011-2016 Massachusetts Institute of Technology' + write(UNIT=OUTPUT_UNIT, FMT=*) & + ' License | http://openmc.readthedocs.io/en/latest/license.html' + write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I1,".",I1)') & VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1 + write(UNIT=OUTPUT_UNIT, FMT='(10X,"Git SHA1 | ",A)') GIT_SHA1 #endif ! Write the date and time - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') & - time_stamp() + write(UNIT=OUTPUT_UNIT, FMT='(9X,"Date/Time | ",A)') time_stamp() #ifdef MPI ! Write number of processors - write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') & + write(UNIT=OUTPUT_UNIT, FMT='(5X,"MPI Processes | ",A)') & trim(to_str(n_procs)) #endif #ifdef _OPENMP ! Write number of OpenMP threads - write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') & + write(UNIT=OUTPUT_UNIT, FMT='(4X,"OpenMP Threads | ",A)') & trim(to_str(omp_get_max_threads())) #endif From 75a7986bb792115295a55ee6f25b40dec962b76e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 23 Aug 2016 06:26:20 -0500 Subject: [PATCH 06/21] Make ASCII logo a little smaller --- src/output.F90 | 55 ++++++++++++++++++++++---------------------------- 1 file changed, 24 insertions(+), 31 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 29950a0bc..dc007aab6 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -38,37 +38,30 @@ contains use omp_lib #endif - write(UNIT=OUTPUT_UNIT, FMT='(/29(A/))') & - ' %%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ################### %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ###################### %%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ####################### %%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ######################## %%%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ######################### %%%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ########################## %%%%%%%%%%%%%%%%%%%%%%%%%%', & - ' ############################ %%%%%%%%%%%%%%%%%%%%%%%%', & - ' ############################# %%%%%%%%%%%%%%%%%%%%%%%', & - ' ############################## %%%%%%%%%%%%%%%%%%%%%', & - ' ############################# %%%%%%%%%%%%%%%%%%%%', & - ' ########################### %%%%%%%%%%%%%%%%%%%%', & - ' ######################### %%%%%%%%%%%%%%%%%%%%%', & - ' ###################### %%%%%%%%%%%%%%%%%%%%%', & - ' #################### %%%%%%%%%%%%%%%%%%%%%', & - ' ################# %%%%%%%%%%%%%%%%%%%%', & - ' ############## %%%%%%%%%%%%%%%%%%%', & - ' ########### %%%%%%%%%%%%%%%%%%', & - ' ####### %%%%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%%%' - + write(UNIT=OUTPUT_UNIT, FMT='(/23(A/))') & + ' %%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%%%%%%%%%%%%%', & + ' ############### %%%%%%%%%%%%%%%%%%%%%%%%', & + ' ################## %%%%%%%%%%%%%%%%%%%%%%%', & + ' ################### %%%%%%%%%%%%%%%%%%%%%%%', & + ' #################### %%%%%%%%%%%%%%%%%%%%%%', & + ' ##################### %%%%%%%%%%%%%%%%%%%%%', & + ' ###################### %%%%%%%%%%%%%%%%%%%%', & + ' ####################### %%%%%%%%%%%%%%%%%%', & + ' ######################## %%%%%%%%%%%%%%%%%', & + ' ###################### %%%%%%%%%%%%%%%%%', & + ' #################### %%%%%%%%%%%%%%%%%', & + ' ################# %%%%%%%%%%%%%%%%%', & + ' ############## %%%%%%%%%%%%%%%%', & + ' ########### %%%%%%%%%%%%%%%', & + ' ####### %%%%%%%%%%%%%%', & + ' %%%%%%%%%%%%' ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & From 6b64cc7783e86fd3c2bc993921ca40a9e06dd8df Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 23 Aug 2016 10:23:52 -0500 Subject: [PATCH 07/21] Add SVG version of logo (thanks @samuelshaner!) --- docs/source/_images/openmc_logo.svg | 60 +++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 docs/source/_images/openmc_logo.svg diff --git a/docs/source/_images/openmc_logo.svg b/docs/source/_images/openmc_logo.svg new file mode 100644 index 000000000..a7352b79a --- /dev/null +++ b/docs/source/_images/openmc_logo.svg @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + From c91480e56c46f65f0e8ca3a4cc5e3035df524a4b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Aug 2016 10:10:52 -0500 Subject: [PATCH 08/21] Will this one satisfy @nelsonag? We shall see. --- src/output.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index dc007aab6..ee5aeeaf1 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -54,14 +54,14 @@ contains ' ##################### %%%%%%%%%%%%%%%%%%%%%', & ' ###################### %%%%%%%%%%%%%%%%%%%%', & ' ####################### %%%%%%%%%%%%%%%%%%', & - ' ######################## %%%%%%%%%%%%%%%%%', & + ' ####################### %%%%%%%%%%%%%%%%%', & ' ###################### %%%%%%%%%%%%%%%%%', & ' #################### %%%%%%%%%%%%%%%%%', & ' ################# %%%%%%%%%%%%%%%%%', & - ' ############## %%%%%%%%%%%%%%%%', & - ' ########### %%%%%%%%%%%%%%%', & - ' ####### %%%%%%%%%%%%%%', & - ' %%%%%%%%%%%%' + ' ############### %%%%%%%%%%%%%%%%', & + ' ############ %%%%%%%%%%%%%%%', & + ' ######## %%%%%%%%%%%%%%', & + ' %%%%%%%%%%%' ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & From 50042b4fe3a12b7a98359b236c073b482e2646f3 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 25 Aug 2016 22:39:40 -0500 Subject: [PATCH 09/21] 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 10/21] 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 11/21] 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 12/21] 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 6eaf6e04d0b739aab40f1f382761fb7b8b220cc2 Mon Sep 17 00:00:00 2001 From: jingang Date: Sun, 28 Aug 2016 17:28:23 -0400 Subject: [PATCH 13/21] openmc.plot can now set 'meshlines' and 'level' --- docs/source/usersguide/input.rst | 4 +- openmc/plots.py | 102 +++++++++++++++++++++++++++++-- 2 files changed, 99 insertions(+), 7 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 4bbc96cf8..2ab823e6e 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -2105,8 +2105,8 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: None :meshlines: - The ``meshlines`` sub-element allows for plotting the boundaries of - a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per + The ``meshlines`` sub-element allows for plotting the boundaries of a + regular mesh on top of a plot. Only one ``meshlines`` element is allowed per ``plot`` element, and it must contain as attributes or sub-elements a mesh type and a linewidth. Optionally, a color may be specified for the overlay: diff --git a/openmc/plots.py b/openmc/plots.py index 73b51da5e..4dacabf28 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -67,6 +67,16 @@ class Plot(object): col_spec : dict Dictionary indicating that certain cells/materials (keys) should be colored with a specific RGB (values) + level : int + Universe depth to plot at + meshlines_type : {'tally', 'entropy', 'ufs', 'cmfd'} + The type of the mesh to be plotted + meshlines_id : int + ID for the mesh specified on ``tallies.xml`` that should be plotted + meshlines_linewidth : int + Pixels of linewidth to specify for the mesh boundaries + meshlines_color : Iterable of int + Color for the meshlines boundaries by RGB. """ @@ -85,6 +95,11 @@ class Plot(object): self._mask_components = None self._mask_background = None self._col_spec = None + self._level = None + self._meshlines_type = None + self._meshlines_id = None + self._meshlines_linewidth = None + self._meshlines_color = None @property def id(self): @@ -138,6 +153,26 @@ class Plot(object): def col_spec(self): return self._col_spec + @property + def level(self): + return self._level + + @property + def meshlines_type(self): + return self._meshlines_type + + @property + def meshlines_id(self): + return self._meshlines_id + + @property + def meshlines_linewidth(self): + return self._meshlines_linewidth + + @property + def meshlines_color(self): + return self._meshlines_color + @id.setter def id(self, plot_id): if plot_id is None: @@ -231,9 +266,9 @@ class Plot(object): @mask_components.setter def mask_components(self, mask_components): - cv.check_type('plot mask_components', mask_components, Iterable, Integral) + cv.check_type('plot mask components', mask_components, Iterable, Integral) for component in mask_components: - cv.check_greater_than('plot mask_components', component, 0, True) + cv.check_greater_than('plot mask components', component, 0, True) self._mask_components = mask_components @mask_background.setter @@ -245,6 +280,40 @@ class Plot(object): cv.check_less_than('plot mask background', rgb, 256) self._mask_background = mask_background + @level.setter + def level(self, plot_level): + cv.check_type('plot level', plot_level, Integral) + cv.check_greater_than('plot level', plot_level, 0, equality=True) + self._level = plot_level + + @meshlines_type.setter + def meshlines_type(self, meshlines_type): + cv.check_type('plot meshlines type', meshlines_type, basestring) + cv.check_value('plot meshlines type', meshlines_type, + ['tally', 'entropy', 'ufs', 'cmfd']) + self._meshlines_type = meshlines_type + + @meshlines_id.setter + def meshlines_id(self, mesh_id): + cv.check_type('plot meshlines id', mesh_id, Integral) + cv.check_greater_than('plot meshlines id', mesh_id, 0, equality=True) + self._meshlines_id = mesh_id + + @meshlines_linewidth.setter + def meshlines_linewidth(self, linewidth): + cv.check_type('plot mesh linewidth', linewidth, Integral) + cv.check_greater_than('plot mesh linewidth', linewidth, 0, equality=True) + self._meshlines_linewidth = linewidth + + @meshlines_color.setter + def meshlines_color(self, color): + cv.check_type('plot meshlines color', color, Iterable, Integral) + cv.check_length('plot meshlines color', color, 3) + for rgb in color: + cv.check_greater_than('plot meshlines color', rgb, 0, True) + cv.check_less_than('plot meshlines color', rgb, 256) + self._meshlines_color = color + def __repr__(self): string = 'Plot\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) @@ -256,11 +325,20 @@ class Plot(object): string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin) string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin) string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color) - string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t', + string += '{0: <16}{1}{2}\n'.format('\tMask components', '=\t', self._mask_components) - string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t', + string += '{0: <16}{1}{2}\n'.format('\tMask background', '=\t', self._mask_background) string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec) + string += '{0: <16}{1}{2}\n'.format('\tLevel', '=\t', self._level) + string += '{0: <16}{1}{2}\n'.format('\tMeshlines type', '=\t', + self._meshlines_type) + string += '{0: <16}{1}{2}\n'.format('\tMeshlines id', '=\t', + self._meshlines_id) + string += '{0: <16}{1}{2}\n'.format('\tMeshlines color', '=\t', + self._meshlines_color) + string += '{0: <16}{1}{2}\n'.format('\tMeshlines linewidth', '=\t', + self._meshlines_linewidth) return string def colorize(self, geometry, seed=1): @@ -382,7 +460,7 @@ class Plot(object): subelement = ET.SubElement(element, "pixels") subelement.text = ' '.join(map(str, self._pixels)) - if self._mask_background is not None: + if self._background is not None: subelement = ET.SubElement(element, "background") subelement.text = ' '.join(map(str, self._background)) @@ -400,6 +478,20 @@ class Plot(object): subelement.set("background", ' '.join(map( str, self._mask_background))) + if self._level is not None: + subelement = ET.SubElement(element, "level") + subelement.text = ' '.join(str(self._level)) + + if self._meshlines_type is not None: + subelement = ET.SubElement(element, "meshlines") + subelement.set("meshtype", self._meshlines_type) + if self._meshlines_id is not None: + subelement.set("id", str(self._meshlines_id)) + if self._meshlines_linewidth is not None: + subelement.set("linewidth", str(self._meshlines_linewidth)) + if self._meshlines_color is not None: + subelement.set("color", ' '.join(map(str, self._meshlines_color))) + return element From b01a8972504edb17b20d6cd91fc9cf776448389b Mon Sep 17 00:00:00 2001 From: jingang Date: Mon, 29 Aug 2016 11:02:29 -0400 Subject: [PATCH 14/21] make meshlines a dict as a single attribute of plot --- openmc/plots.py | 118 ++++++++++++++++++++---------------------------- 1 file changed, 50 insertions(+), 68 deletions(-) diff --git a/openmc/plots.py b/openmc/plots.py index 4dacabf28..cc5c0d44b 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -69,14 +69,9 @@ class Plot(object): colored with a specific RGB (values) level : int Universe depth to plot at - meshlines_type : {'tally', 'entropy', 'ufs', 'cmfd'} - The type of the mesh to be plotted - meshlines_id : int - ID for the mesh specified on ``tallies.xml`` that should be plotted - meshlines_linewidth : int - Pixels of linewidth to specify for the mesh boundaries - meshlines_color : Iterable of int - Color for the meshlines boundaries by RGB. + meshlines : dict + Dictionary defining type, id, linewidth and color of a regular mesh + to be plotted on top of a plot """ @@ -91,15 +86,12 @@ class Plot(object): self._color = 'cell' self._type = 'slice' self._basis = 'xy' - self._background = [0, 0, 0] + self._background = None self._mask_components = None self._mask_background = None self._col_spec = None self._level = None - self._meshlines_type = None - self._meshlines_id = None - self._meshlines_linewidth = None - self._meshlines_color = None + self._meshlines = None @property def id(self): @@ -158,20 +150,8 @@ class Plot(object): return self._level @property - def meshlines_type(self): - return self._meshlines_type - - @property - def meshlines_id(self): - return self._meshlines_id - - @property - def meshlines_linewidth(self): - return self._meshlines_linewidth - - @property - def meshlines_color(self): - return self._meshlines_color + def meshlines(self): + return self._meshlines @id.setter def id(self, plot_id): @@ -286,33 +266,38 @@ class Plot(object): cv.check_greater_than('plot level', plot_level, 0, equality=True) self._level = plot_level - @meshlines_type.setter - def meshlines_type(self, meshlines_type): - cv.check_type('plot meshlines type', meshlines_type, basestring) - cv.check_value('plot meshlines type', meshlines_type, - ['tally', 'entropy', 'ufs', 'cmfd']) - self._meshlines_type = meshlines_type + @meshlines.setter + def meshlines(self, meshlines): + cv.check_type('plot meshlines', meshlines, dict) + if 'type' not in meshlines: + msg = 'Unable to set on plot the meshlines "{0}" which ' \ + 'does not have a "type" key'.format(meshlines) + raise ValueError(msg) - @meshlines_id.setter - def meshlines_id(self, mesh_id): - cv.check_type('plot meshlines id', mesh_id, Integral) - cv.check_greater_than('plot meshlines id', mesh_id, 0, equality=True) - self._meshlines_id = mesh_id + elif meshlines['type'] not in ['tally', 'entropy', 'ufs', 'cmfd']: + msg = 'Unable to set the meshlines with ' \ + 'type "{0}"'.format(meshlines['type']) + raise ValueError(msg) - @meshlines_linewidth.setter - def meshlines_linewidth(self, linewidth): - cv.check_type('plot mesh linewidth', linewidth, Integral) - cv.check_greater_than('plot mesh linewidth', linewidth, 0, equality=True) - self._meshlines_linewidth = linewidth + if 'id' in meshlines: + cv.check_type('plot meshlines id', meshlines['id'], Integral) + cv.check_greater_than('plot meshlines id', meshlines['id'], 0, + equality=True) - @meshlines_color.setter - def meshlines_color(self, color): - cv.check_type('plot meshlines color', color, Iterable, Integral) - cv.check_length('plot meshlines color', color, 3) - for rgb in color: - cv.check_greater_than('plot meshlines color', rgb, 0, True) - cv.check_less_than('plot meshlines color', rgb, 256) - self._meshlines_color = color + if 'linewidth' in meshlines: + cv.check_type('plot mesh linewidth', meshlines['linewidth'], Integral) + cv.check_greater_than('plot mesh linewidth', meshlines['linewidth'], + 0, equality=True) + + if 'color' in meshlines: + cv.check_type('plot meshlines color', meshlines['color'], Iterable, + Integral) + cv.check_length('plot meshlines color', meshlines['color'], 3) + for rgb in meshlines['color']: + cv.check_greater_than('plot meshlines color', rgb, 0, True) + cv.check_less_than('plot meshlines color', rgb, 256) + + self._meshlines = meshlines def __repr__(self): string = 'Plot\n' @@ -325,20 +310,16 @@ class Plot(object): string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin) string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin) string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color) + string += '{0: <16}{1}{2}\n'.format('\tBackground', '=\t', + self._background) string += '{0: <16}{1}{2}\n'.format('\tMask components', '=\t', self._mask_components) string += '{0: <16}{1}{2}\n'.format('\tMask background', '=\t', self._mask_background) string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec) string += '{0: <16}{1}{2}\n'.format('\tLevel', '=\t', self._level) - string += '{0: <16}{1}{2}\n'.format('\tMeshlines type', '=\t', - self._meshlines_type) - string += '{0: <16}{1}{2}\n'.format('\tMeshlines id', '=\t', - self._meshlines_id) - string += '{0: <16}{1}{2}\n'.format('\tMeshlines color', '=\t', - self._meshlines_color) - string += '{0: <16}{1}{2}\n'.format('\tMeshlines linewidth', '=\t', - self._meshlines_linewidth) + string += '{0: <16}{1}{2}\n'.format('\tMeshlines', '=\t', + self._meshlines) return string def colorize(self, geometry, seed=1): @@ -480,17 +461,18 @@ class Plot(object): if self._level is not None: subelement = ET.SubElement(element, "level") - subelement.text = ' '.join(str(self._level)) + subelement.text = str(self._level) - if self._meshlines_type is not None: + if self._meshlines is not None: subelement = ET.SubElement(element, "meshlines") - subelement.set("meshtype", self._meshlines_type) - if self._meshlines_id is not None: - subelement.set("id", str(self._meshlines_id)) - if self._meshlines_linewidth is not None: - subelement.set("linewidth", str(self._meshlines_linewidth)) - if self._meshlines_color is not None: - subelement.set("color", ' '.join(map(str, self._meshlines_color))) + subelement.set("meshtype", self._meshlines['type']) + if self._meshlines['id'] is not None: + subelement.set("id", str(self._meshlines['id'])) + if self._meshlines['linewidth'] is not None: + subelement.set("linewidth", str(self._meshlines['linewidth'])) + if self._meshlines['color'] is not None: + subelement.set("color", ' '.join(map( + str, self._meshlines['color']))) return element From 813efa3b66dd9b89378621bd0057c2520287b152 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 29 Aug 2016 10:45:34 -0500 Subject: [PATCH 15/21] 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 16/21] 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 17/21] 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 From 0fbe8e72f1dcf09178e5647e93797c2c866f3fa3 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 30 Aug 2016 12:04:34 -0600 Subject: [PATCH 18/21] Bug fix for MGXS Pandas DataFrames with selected nuclides --- openmc/mgxs/mgxs.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index de7475308..d0bc07987 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -1496,12 +1496,14 @@ class MGXS(object): # If the user requested a specific set of nuclides elif self.by_nuclide and nuclides != 'all': + query_nuclides = nuclides xs_tally = self.xs_tally.get_slice(nuclides=nuclides) df = xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) # If the user requested all nuclides, keep nuclide column in dataframe else: + query_nuclides = self.get_nuclides() df = self.xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) @@ -1513,7 +1515,7 @@ class MGXS(object): # Override energy groups bounds with indices all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int) - all_groups = np.repeat(all_groups, self.num_nuclides) + all_groups = np.repeat(all_groups, len(query_nuclides)) if 'energy low [MeV]' in df and 'energyout low [MeV]' in df: df.rename(columns={'energy low [MeV]': 'group in'}, inplace=True) From b80b0fae835e920e1fc2a7686265cc37cd200cb1 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 30 Aug 2016 14:29:19 -0400 Subject: [PATCH 19/21] Fixed isue with MGXS Pandas DataFrame with total nuclides --- openmc/mgxs/mgxs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index d0bc07987..dd5c2056c 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -337,7 +337,7 @@ class MGXS(object): if self.by_nuclide: return self.get_nuclides() else: - return 'sum' + return ['sum'] @property def loaded_sp(self): @@ -1483,7 +1483,7 @@ class MGXS(object): if self.by_nuclide and nuclides == 'sum': # Use tally summation to sum across all nuclides - query_nuclides = self.get_nuclides() + query_nuclides = [nuclides] xs_tally = self.xs_tally.summation(nuclides=query_nuclides) df = xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) @@ -1503,7 +1503,7 @@ class MGXS(object): # If the user requested all nuclides, keep nuclide column in dataframe else: - query_nuclides = self.get_nuclides() + query_nuclides = self.nuclides df = self.xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) From 627057d87717411c94475e53ac097a33a3dddadf Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 30 Aug 2016 16:03:29 -0400 Subject: [PATCH 20/21] Fixed query_nuclides for sum in MGXS Pandas DataFrames --- openmc/mgxs/mgxs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index dd5c2056c..543caf088 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -1484,7 +1484,7 @@ class MGXS(object): # Use tally summation to sum across all nuclides query_nuclides = [nuclides] - xs_tally = self.xs_tally.summation(nuclides=query_nuclides) + xs_tally = self.xs_tally.summation(nuclides=self.get_nuclides()) df = xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) From 96acab8ba1f8e34200a72a66fd768d9bcaed15aa Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 30 Aug 2016 21:40:50 -0400 Subject: [PATCH 21/21] Fixed bug in nuclide summation for MGXS Pandas DataFrame --- openmc/mgxs/mgxs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 543caf088..1d2e9098d 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -1490,9 +1490,9 @@ class MGXS(object): # Remove nuclide column since it is homogeneous and redundant if self.domain_type == 'mesh': - df.drop('nuclide', axis=1, level=0, inplace=True) + df.drop('sum(nuclide)', axis=1, level=0, inplace=True) else: - df.drop('nuclide', axis=1, inplace=True) + df.drop('sum(nuclide)', axis=1, inplace=True) # If the user requested a specific set of nuclides elif self.by_nuclide and nuclides != 'all':