From f3cfc181f47d959158d5e871209c45b9aa810f33 Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 16 Aug 2016 20:14:34 -0500 Subject: [PATCH 01/23] 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 89e0d8aa76..c1ed9336ab 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 04119a2dc8..0df6042f6a 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 ea7da21edf..15107e8c84 100644 --- a/tests/test_triso/results_true.dat +++ b/tests/test_triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.662675E+00 1.475968E-02 +1.636336E+00 1.154000E-01 diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py index 9a8fb0c3be..8b6296ab5b 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/23] 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 8b6296ab5b..da5deef00d 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/23] 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 392c2c72df..14f4a2128d 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -334,6 +334,7 @@ Functions :nosignatures: openmc.model.create_triso_lattice + openmc.model.pack_trisos -------------------------------------------- :mod:`openmc.data` -- Nuclear Data Interface diff --git a/openmc/model/triso.py b/openmc/model/triso.py index c1ed9336ab..763fc0ddda 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/23] 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 763fc0ddda..33c955e3be 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 7d07b485b9a79f04fc2d34c2e1d27795befa5be4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 20 Aug 2016 16:32:42 -0500 Subject: [PATCH 05/23] New OpenMC logo!! --- docs/source/_images/openmc.png | Bin 8862 -> 0 bytes docs/source/_images/openmc200px.png | Bin 6919 -> 0 bytes docs/source/_images/openmc_logo.png | Bin 0 -> 15452 bytes docs/source/_static/theme_overrides.css | 4 ++++ docs/source/conf.py | 2 +- 5 files changed, 5 insertions(+), 1 deletion(-) delete mode 100644 docs/source/_images/openmc.png delete mode 100644 docs/source/_images/openmc200px.png create mode 100644 docs/source/_images/openmc_logo.png diff --git a/docs/source/_images/openmc.png b/docs/source/_images/openmc.png deleted file mode 100644 index 9f5e97cd6e844052caf6e6095a5088922eb89113..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8862 zcmaJ{cRX9~_ZP88@iB@TAGG$~wQH}~sz!?1yNFS0wf5drYKvL4c4MWbs8+425h7Jo zC{n~8-~0LV_s{RXUPWh=%qUFLlUfMU_}8O&naA!fosaA`nJKq{qz5QK$?Q^aG;S5rt=798Q=~Ja}IJN z3JVJp_wx4%c5!~{CLR#v0pC+$CnEZn$UysnRru@O(y+d}I;E(=IhjKYI?5hi1a-PC zpUn^ljSpo2Jy2EAnPz0Y8weU5f@UC3@xxKDFA;j%)<_^UHuhQLom9alPia^zYuOnpMIM9DZ~ZIq_VwrcFt zVf~zO^Tc{Yts2zzHPdzs;UqnjJ)jD)yBIC;lo;>tB=lKRCi zyXLE2IJUy_P5X>c$OVJ(XOj2;L(BK~KD~$a|C-i%AIy3mMEaKXR2?fZLtx2rQf9k} z*RD?JpRSIV8R=2(of`4{Pw}-_>ZG1kg zzqDrO3e>~YOv5eafDOe(IA+xrv9Xn6H6ELpa!5b}vfllniG`_xo_D#=CQ~qT`A8NTTUfKq|NU%k*0JAf&aR$O?)4m$5^! znqQ?jkG>xW;vz2fx167nXi*v9TH(&W6DO=n)upVk^znOOp_SdIx9%@k(A=8c^7Y3>GQ}H$?e6W7So4uNY9OsAxI3 zLUqA!&biOJeo<_SVmbu-=@UZ0N0cMRRXM4mt~EdlE;<}YhhLJ!+syeJGXV-6pm?tI z6B|QpLI=JHpd?uEK=<{JvQNf}y{}OAA5k|^0!v>{e`_HzLwyw&4F~;_cy(FbtNXXa zXtoz+XfE}X7cz;&iA)UZBV_P6K)E;#I4PFXbH(vB*R~XP3wV*0!Y~yuw>pqx%6)7< zN&x88#ZuKI1=imU`VEp1vp%2B?QXa?kxqk)$eT_~A{gzwM17RRm*07*wvzBL=KJ9{ zNxTHkN%=R3JF_83;i(<#>3nG8L`fGlDh6B*LW<3m5#6E%$J`q>UgcHJPy`WVP@F;_ zQhKs$#Mar_zNj@nkxi{F&)8VK3#B`{+#j1)y?nc7hcnGqBku9Wo|SSPgG9yiqo|d@ z{~n5#qPdi)3YqF6|DcIo8sqpVT2lqTz92VkPWS(@;r6D0Kh_KA3de)?+>G21{s7 zju`hb^*lnBdC~)FI6z5A`do;h;J|%#mZcCa*3%71yy=IQJ7~JobHPv{mkR+CeF%w~ zyd5W?d!@CDcEuUK&e=o~@sU3a(j**=WY=fA9n+%I`R5eaP@{qv&@F|G8Z}}Lr$wA4 zuxz@n_BwT0PaSo4>W!b9As;WNd>gC3hqF7I+=jU@ojPcqF*%vT{hot*{eH!y7f=I2 zz7tzLEt1@+Lf6Mnb9twm!6N4+{~0OFUJ0!KI-_$(UG-!Z5Feq?^b{ujwnVsp%lgaM4?L4XW#fcefA&vCm4v!1xCFOtciR zr0M!774Okp&g&sPD|Gb@=Fybl|W%yfPK%`CVEuN?>`C69vdJS;uB&>co@YamKFRCW{^-w%sND z3#u4>yI%mbC|Yp(yetQC*O|w3Jwjrf4Z1962%C>meSEX4by5Z`&xd~cJkG%z(QhU?FXiatJ|EGu;S=&GQ4W!!|h0W)On{xk78N` z8Q~~KrP?X-=DcEK&B@K_P|i^3m*HFbCl>xQx6(M@59PKzBj4qtQk{)p5rvYg*t0DM z`TRY**dvl}G9gnllBgk41B|y%5>KwQ5M9twE$F#I^OLU_A|_4+UhQg=Zy6iH>qJ#m z_NB}Bv4_RSwbty2fusEd)fn`p@t3Nbq%XRw_#}L_j<`89ST+OhV)g_3Vg2P`b*>3& zo**W|RA6i-l5yUiX(^>UeH6Xaa>I+{4^4TlKtXP1 zCR9DHLnXt4f#O!fU*+wq4O*}-e>c3L>3frc435piy&bSbVt%4T;0tbzu>J%9`f?mX zdNNWx9Z?bIDGJY6ZLI1dlZX##H*6B(Zl<h1$@L+ccDTTxA?Bfi|Oi`#sV$4;*A z&M6qrlFQj%SIt}yW;EYxjKu6w6407ROrKylS_{?iovNNJp2aSqtZ!EN5P$hq$HR2t zaLc}TQ~dfhJ@03>?}pF5lzcJwM&U{gR(#5~%B*%I)99ZhU8X0j;xHRn5oBlt?w!eg zy102UJ!QQ&_hIr!lGhkE^7Y!?gAK|^?!(@_XK8f)Vc&$I0rH!8O*(VVE4d(VSD_Vy zHC=N$yDkOgjYz(KY0)?);B|xdATOz#e-&~6%q`e0k5P|NmnW*LaYK@KRe*_o6%?l-E=qwh z4Nx+UR_)p8g2bwY94n@7_kzlLJtT$~dmis{R9Qz6qaB&19(RktFsd9YQnRzjA!5k@ zagzXFt9Dgcmhwy+{!!#$cYnhxF%z<816xZ#dYjFFT%m(@k9Lw-C*^$@nH#B=$zfqjPc7)o30L z;z)Yff-YMK(74d)VszMz+*+`ge`qdL>c~bDgAUtoP8ayHH~E?lTs4OZDNl zx3{|(Cf#_tVUuq`U|{2(fD|7HY>J&ZM=xEzY?#(caVUDqZ*r4OzZQY8l7g%f)$B}s zx->&JpYzaMY9KKo#d^vY72KM&#im|MF>FuM5vP;R*m;v7|524(N8%a=3yH)F_9-PH zbn3VI(OEC>YzJ+w7H6HRvkh1xi_xqTz<&T48I(S5Q5ki%>L@_G0e3(&;4bu-yic0u z(#;tQC;2Wfh>&NoQ9a80PgEH^5DTfLip zEvOJubRUpOUqzZ_YGnc48}^7jjLi7YaTyw7s@R~L+Q`%KIU@IbQo=xhy{avkqA<^SD0dQu4BY&DLoR);5m*ubp28M<+n_H%%& zBfE{%$x)$%(~Hr5#3he3N-WUw!BnNOsbQN&79Y)JOvMDfvA@m@KQv%pb`;3FB;BC~ zY~esqD331Dwt9bL1^9F|JG44iQBDC&NLdB~p3arTP(gdrj2b2SShH9zfPTuWk?g~J)Dzu-Ef+}aSqa*2@3 za;#45l20$x8HjZJ!{Z4l_Qxv^p#oSZbZ#2_T+f4W5%S7)WHmeeOzk%uC`h6?Q%=ye ziPja{3?URU{bm=*IgxM8`Qd|TGtUq2XUonOCNHnFty3(Q9k|2{!cq?Ew=BoT9D@@z*TF~mb zU#&CmYKBbUNamJGO5nOL)y2W)Z^3|w<5sKjW>$)mrWuW2`%F)Sam{~bg(s3yq#0J6 zg`x{!{dZdd!R)4mBBTQ?UmQj&Yi@sf>Igpn*uC+JdxGZE?h_w55(lUjQnXLlAYYfK|qGy43}z!Tkg*Vp7zAfM0_od~@5tH%^wwxZh0C?;l)XSa4Kh0sy`P^t!j zX~UHWYRvcGy}(!leO2Fu!b0_qcdZ3LQfN5eTcoRwx-usCyme!E|0eq7*akypyg#L3 zrPpu^I8SoRn8! zO)8H4nK%$8{d3&=PiaR`jMm}~kl(SAU9Uov4d>C+rL42V8GyhZ(u05mApzDuCtJc; zBc=5gdXogMX8+O$eeJn1zV6$rDBa9GuRm}9B2&rCR2a@r&1n#{H2f~_%`Z6KnsI5| zUeP?cP=(hU0t^14L~CiB|#-~1EF zd_GrO4b@k6IDaxjEbMo)Q`JJNB$aFLD0?wkTRV;7uGV0drYJ5XDS=r+zHc~5o}Bh zt$kkmx;S^&r%+b3#@D1^G~=P39nXPb47}2l^EeNlm;-ww&)I8_Ec08Ewqr9|o4yWK z)l9DC16iFb6iTYk>7~^cY<6l_H%Qo~lDZ$}n!03_y%UA>GVD?`>>doXjK(=DbIo!b z*N^!`_5~4d61f+3xWpQ?F?cF3=<3LGcWgfSB0}0oT<}k$gLVO#IRa$$j$5vyU>wUL z42Cz*_|71!1_)MNE!BN)v>RegLP;V*0nud#wSYBUk|Q!E&6e=82!d(#P~^z90<&>f z!A6x<+fUy$2_F~7yqUfiQ5)Ze5uU#uGX`YwdJxe^MHeCVXaUS#*6sy{*c3-d+}$s1 z>0S|>@+s<+`1tza)ecHTw@w+Cl3Tp=yi(lz5YOi|=6%pW1tCjd=qhNEkt&wWA1Zi% zP$>6Iu0!s2GoHCV|De(0o3+i~quy~5t6K;v!2y_$P1|kV8o}}DxYI&2 zc_e7I%Y?^*QeswrJI43CR;d?med~g4p}L#Dd26n_R~Lzy)#wL# zRAp0vWAu?2iP-blG{&i>!oQW#Bb1#n>seDV{{4jG${9i?HTcG+dcN7}IZ5!xVR&_p zvvT^n<%sNoPqL@&J>-=9l;>OID8)T;&!o1%yPqhi8b}KAVXEYZ+n#oS`&HMYenpzt z#)tsod9k+`zr`C<14>QNWs~acr3XmN^4P9$W0Zrywd zQ+%~h!eWs_^+6U*`|pNiIwcZ>D<9(7k~bFr8Y8o;h} ze7-S5XsL(X@6MP?X-i&`E*apu*mAQ#-vr2Bzc*p%*>8xTBUR$D)BXj$=qdi&o40C;Y-$i=Lrna=b zDLCfc8@{|jq>2G7X%O~H;hL%8d^^Ncj}#ZtzoqfQTDeaFOlLSBNykbv0fcIlN(0>~ zvo$-<_@*Ycc|io#KbUn#O7Q9?)f5%ygb@qd>MKn$b+`6MEk7f^Ho0&Yw9K2|8=RHd zC~kRA9^L)+SG3J0|HL`0gkDT8Un;#Op2foN3B4{cy-6+}fL7%tW6U;VwVv$LDAv=ao1h(POsB-4y=EuWbYAq@^}5gI`dD`v^{1>Ol21iO}W6y6e9}V_u>5d zI!8Z$9Z!Y5Din^X8oKpi3+v9kh5bA5@=OTMs;Ov^4PQ2|C35s@cO0VN37wd9I6C|I zb;ToOQrUBL4@!5Y>HK5QW`p14%9JNwY(`tGE5{10B-FGgq^qwZg5nP=d{`HO!k@j> z(U+V4o?1pan;!sB!dB+ARSZ~rDfK_dH~{lI8N&v>YIY#xEu_B02VVg~QjV$t@HbU; z|H$7dH6wWaPdCQ|Ll4rHP40W=?zQ-kH(u+o8+*cTBo3VX4Ef-6Sb87+&w) zORstS!ZvKNFgx`Wgi~R>E|&iTBlfF%pDHcVPj~o!i%c9yU-mp{G^+FwVNWKY-O?8y0GL*g8jIxy6LwQD&I^c7ZGhpD z;ctJnoqe$5l%yB(t|`PB>&Ru@Wg9RF7T5OJ+D&eVy89Xv_wCyOZO6}7DFrB>DtOr4 zMbaJ}2>LyGSJ{VK3=l+==?Cinf|#*Mx8NsFoK}ms>4^6P*{ddq`~f;& z3l|}iC#?7LNr0}12;@df9sf<=G`mC+uTU`?1=r=jg@@oG6b$B-wLQ~E88nRrzdm|f z6ZxTXb!L^)qZg-{U9 zMn+3jBHlDY@RLT-9VdHm@$LoymeYR+QstF0ddVEgmxZAFc0jejxap=*2jE!$WARlo zYmN%NqM*=XxC|Z0i5FlR0T)tSCa*g7UJ{08f|LdkCA)bFY=q? ziSd*te`2UqrP3-ua|uMt<+1XlTdoH=iiuOyhLNL`f}R%o)&SgKB=L?3@=XB(Jq{c? zSmL9s@EF!pKDJ*Cagz8iz+>t1QJxArJI$y#XKZ(V8CushJd!?wv?IhSVH{id#A_wb zhE_HBH$2scOIc57&hYtZo8&7JMKaK==*KsTx6ioX|2d=!JA;D~?ezxpJGo)Db?D?C z?nA8xF&^uC(u%4Se&URN2ZAD!JRwi93m68#hKjTm%awnH3C z0i$W`k4;#ZppdWbyr!#~wOge>W6KcXy?3PtKvE``ZT}tBjpi^G@ha01G1Mn#Y?kYr zqor%1qtvZWGlaq_+Q6bJZmyiiy0qlOmK!sCjSDB;fQ->b0GR>x8k!Z!x9F;foG=p( zussY1Y6{fBYCKFEm2({;HNqQp*6;QuDD-%O3>8i8NvnHqg}HWjmyB59y7UYD7(Ydw zt35r_?_>zPpCJ-th{V*muel%q_6Z$}uGxuV5~Ii>x}_r#E{58xTgV!y+t-rNgaVR1#BQ`<&{Iq;^d>&MbfPAn!?EB!)Xzl!^ib! zWvv$uH~pnoR9o=6g5$>4K_SGz)NsJ#0KxR*uk+RkTmYMgR@{0yJKPXenCl;upHVN`o5vXFUGB*?jNx4 z_%EwXAb*IGXIB2LP43;A#gGIlXf21?2`WIp53P>E`^|NA#TpW(sweC zB$K)K{PvHTB=_FAcOH3kKYljf&nKVE`JLZ+%*_3r^Lw2$gfYS&T)RY)79`+QfMxdr_JXRX_N-smq28q9C(2W58 z2*9)=F|de$r_`G89T=a$O<@8sEp$Zyz!Cre(^oMG-${9}=_CSOnHi&MJbn2>k?VvI zF95X6ruscFwSZ}jfT-0aFI3#6({zXh5~3Q!*8q4JK+M{ZfW|=$+td|nms#l92-HS` zQt}f3rdz5GfN)2gVl9QLNfk8>Cu$I10V52c#8k<6D849mH6D3HFQl(2Ogz+@raG`3 zQ7;I5X{@XL*AqHMD7+EC+Dx%O?p;_u7DL@141bRaU$Cj_V8X45&V|nvq$1Weczpe0 z4H0ev@P13}5h>+`{XLs{02={OSFBwIq5l!UGy_?JaF?$?b>~Sc)Y5?w8?m{h0As^)SC{ZW{OqmP&` ziFdZ$G9E8fysQ1yxOd^XApF=!C!o<6Y`lJ~{Z3(gAfVbA%KqSl1h*a%05bfJ38!T|tVH>P6UZFjPy3jp*PYAOiwsj{l(IVJ(oT-px-oMIrw z7~e5ZR?ovZ6le8zU>P&EuWAvrcfN%EFb_tGV5U+6TG3O`H*PL z>K|ON9pbIMY91IrF;n@O_%z?&_rK8g?)Ya3BzjtZFCjh%;IM%_L%8(wMg&ByX{ux3 zN&|x~m_C{4YMD5n_3(tc+inK1*+6m=GdCynceK7Fh(7Fb^I6f3)@~qUnD~PY=P+TVf%l+C z!Cmpmm_-K1wr=cGYnpBV>@N2z>Tr?b1MKu6m8}ZT6_OXTF6If6 z2`Tuk#Qemolj13>Yw!dH=?uwOz^SB4da_@p`gZoe`s^D3inG^z6=%4FH?vX_LmSk! zfulVD%-$aF^m*ve`nRBcu)P7}(wO-v9=SoxWhXf}Qqk#s~euaBZx+t)tj} z0YGT@3ILB?9RWV|?TXXlz>b`KNPGom0se8UW7BJ6nNOjJb+*1>^Z}f;wtJlrwd;Y#- z85UA&!u4QY2My-|L_`b*fh0lQALci=qtyfc14?({Py2l`^#>wh>4Nd7g$ zGOcvpnO{?;E&6BSDW3)50$75j4qwHh9XuGlr~mb*TL7dAeIG1PZSfs43=gCv+ydaE zh4wd$va04efVS2`HucSR8~_M}8WyXe@D?zB22dUBRfY-x&8%R%Kh*eL zEUPcE;mM1#RpGh*P~^A7cn080VD-)bfKWhg&m31^cBEv%(icikT6H%Nn3mm}8M1fJ zZ=k0Ri&8&CsUM1Ox?h<#>w^FQ|Mb&JPON?0J$LyHqHKCvH~@eINO6Nafn}n6?!9Hp zFWZdy^SuSXA)|rLNCIdvRuk}%vYN)CBj1W!zY&0O9umNProJhU9YV@K8toSXejWc) z`=MgzGG2(lz`kk#@1KtRSh(HW3YprlpV$)x7X`FVe#iYr&K=%|GXu6hS)C?NA{1l z^+1sObeDH#ZRw2DepNE_w3Q=*!3@?@02{_60pY3mSXpKD-N=(d1OUEJ(<~U$G?=hW zmwX$5A>&&RjG{bRzEIOFX8yE={Bt4bKP%&k1B! zg`Eb{S$2E`e%n}P!et=(4FNs_`Bh@v!h$Z85Y-a$A^^8Bu*=e0jnw(sSp8@QT8&-!gppaL!BXw z&w*Apn0mUT~zK@fi- z-nr=CZBh-Mz_R+=Re#{iU|tVE$&7hc!sGsT0H4Y=wer>LD)m(IQ2@q7wqwk6ZKAt% zGqUZAdrN1Wv;>tWf5%lm=O0ON10cAjzpuJvFuGyronx27!Jy{UkjZsP-oUx;*-I+{ zfb|rH_HW;kI=rh1$-cjias}@3E(zYRR4%*-03Z}alm)H_@F(QRPLMDQ2ot&iKw)v0 zPO=ezJ5wc~XkP_zg;QU{P|s;RK@QP6la5RLvv_CAUjTqY;q5%?;#Ji>>Zx(yRhzB| zOs5!>WMysbgs0OC70qwq0%(4giQA)emjIiBpM@+&FT3 zE>LGJx9)c4NeXj-P$B99fWGdQT?{;EW~R#np@zk|`U3z?(ibf|x=zHmjdp_`6F-E# z9a7cowb&o;Hpcht_NJ9eNo?ytCz(%={H2 zO$fTpi9b*qIR&6m4<92b8J894Y8*)J?)jp2@QufZZ3QavDK2$(1JX`|EU=z5d7y+R z?e!gccl+z9V|%U!aL}&Z(6Ixr>-~Kr?pP8oudJLO%r!YON6-&PJ6cBsPZ4w{0Ha%6 z(IIZm)t^6j)pR0UV;TeWu<^4ZrXF)W+u9+465aIuGjPm~?|48k!Srx~^qI$d5^Ro? z!lk6Lwo_>lqBQ8&@vhd**>qKH%*ny?nLO|->SGLOdpw zo8H9Gdq&}6L5^NJ7nxNXX70|K3elY~tB`J2I{24`Ky73i7}pwU1pL$*&q{dbzX})^ zWIhvEpY3r=r;!-$GT^ycJjV$xFui7u9PrZPJqa>3mYRUn1xqQwDFMj@)`gUs zAm#u|4U}=yQi*yZ(yFaFpNOT)7*Est_NRETw-wfCx(%dE9R^E)GCl|<)rpq{uP~bL ztfcJ^BjTMce^En`M*zOuz!@%E6~5h?Bq$|c53?4a5oEUNx@EnFQ1~SPt1=Y7zxM4V zu?}RZtVDR78N+G!GyK7;rpLPPcTlq@4=$m|Zk8q?z68KGkJkkJwGF4SH0>Z+UzR=0 zYC%zWazhfq@HUL5*PLc$Q8bx$t&Tzk@S)M+1>h2(jI+lnW5dpmDZGJ70P_I6Wva03 zAll(K-_{3XKg;=EhNnYhaTquzkytZ$SpjSp8W#E;B*hvGzdb^!&~PJw8-{r*&TOdA zfX|!f!Vg6|H#w<v*TI5lt-fQ{|tgo4%!f;pRZ zXO>GngG9%s?I1jAAvdcluN<*cwGwe1%z=<+kXT9r>%^=Rq)w2! z0Lx4{Qzvy{EQ6(hlq3gx-7u$O000D;RVkKw5yAl!a$2B#MffD$aFm+z?AVg%Ho<(U zq4P4<-8TTZxw5LcJVn}dra1!-Mmt*Hbm~tby+;}X{mk}fLt@j^r_`f~qZnPNJ`$*n zoRX{W$$$$06xLfL#1}ED5nNNgBD|rmr|m80rsId-rG1`(S)&yCRcz}r1|x@NbDnlNJfMrhN6f@1(S$&w@=mq2^NboC$@|u1#;Mp>QuoUELW&-4erkazm1W z;Z24!=p|C;I{;jc;|Ku&(XLH9L2T4@OEluUF^LD5)KKG62EJh+F$?{U!h7H##yY>c z1`mKU43xiES=>5IS--?tDKLz*mq$Z7_?SdGLf(&bs(TFSW1PK6`|*zeID~X8;{br{ zOD{-hrni)3SLB3c7ZkR;B3n)nO4jq*d|hiRfarCj1Rd0xHgI zr~vSD7|shYxV^}O{Br`ekuhB^;2#v{!;mk;K-QUuzMi&SFh+0*7iR=(EU1JzMqm%? z4gyj585W{K4G^P%4EFaY@?OK4IXv_`BQ3MC8|}1|OHx4#9D_}Eh!pmANdMN2`vH8{ zKq5%_^NN~=6AKxqSOkKRivWI5acTdfhQgl(){7AV$)p(EYHGo_)fWt3HIb8;7P=y! zl`4|Mk>IN|<6RdOQ56fhKI zA5jN@xL4i@;|@*mTkmN}mV#90OT{;7^Xx$9<_9j1+Ynpsyx;0|BAoG67qoT`iq< zyv(@QSSH}>{_60lcDzn-pK!|6C3BWtQ8Ig3sPv@Td);%FKaH|7hn>&5gD3?Hh(a#2 zIwrDJEPw6WOJI&y6zEbJ?D>$uC+HF(re$vEtIY{h>ztbP<41He*id$>)5>wczwkqO5YnODE$856ZrMx4<#^&a6@u zrDA0n^w~Rz^$si-XGhR;A>P@t1)x7!$U~Y^nC)Q%lgUso_XpSJ%oa{`xBrd>T>`*( zRf3m+Z>cWrO|_=+^Tzm0E(B^Lr}%>5yV#|@3Cc?2f6h0M`p9rLhwVe&_xydwaN*}Z ztR78#55P4Tm2Lv#G7w#+4(#bgDEt(Nb`W@<&hntbBqX7Nh~5uIRf@GbhihN<63Z*% z-J7f*b}@Jz={vj>03Hx;Q1|#hhfriQ5pM2y; zblm@r8f;tx&=ZFDGK)cvs-f_IL)RW-J@tl6C6dtP0w_wAaC<*YZtqNEBwh|iqsjd{ zZqfVRbw2o`P|CE9QsDwZBy>~cUr}s4U}QXEduS|fdZ2mjo*yV!n+xd)&T z*D;ek8^9Ymnjq29`jju!w2B}%!+d4v95DR=y8GYMP-8bh&k$obzy|<2A{8nZ%(GeI zYyhjJ_bzhZ2NCt_cU8|*UYUEs+L$vAMWFWC4v5&aOu7aWhG0>CH0>El)c0!%Rd zCWFnUDFL(yK{rP`HrYmRi@=irbpd=oMt6lJ0N*;2vq$W-Z+8fVK3eD@oOpNpbG}g1 z#RU1EFzfJVLvfu=iVI4Kv+XVWkM7qF?)aj9bXU&IW09p|DFt7^0C$O19W(DB_4Nw& zjJm2>6mEYQk?3lBS`9UJ0#pNFNY~tV=51sHNz`N9;|Q_}=DX4y#&3BARn65tMZ1-N z>n#hg$jqxIooJKEkdc`qfzCmJy~MaK*4>tMtT*r3;=|pHo9tqC@4Zy=J zHY$cASW6D+$M!v{ z9o&BY(1D!~BTvMh$YmA?QMi$wrD9o(cHchVQW#v@#P<1Dg+Op7m4whXuyx~}JnZ&; zxJ{REfE`j!Y_WA?f1L=A3k!1qNHRp z4nk%tE+UdnLeWRgrh`DjtYQgC)kf*$B0N}DM7vv`P($GkUFR|F9zgmo0$(WnHUTGQ zNW+GD!3n4t$H2PEs^*(Cm;P}k{wNsdA$K0#y9E56!H>nKsgGwgB}a@aQhq|+HEJT6 zTui_$(&2|79W0M261)3OjP0Jt%i79<;c8(%vUzvBt1T;EOkH{9Ok&q5EHRgbh>H~O zioMwO2F9BCn+tUfo*<;44 zus=nmBa3w3^Y>+?{EP+OP}6)x^7}~0N@h_a1;-WQK7x3ySk2*Y<}Tbq1!M|Sm`LQc zL8dT;|IaAuS8k^;g)zd7!c-?KkK#bjSGL+Rg}k;!J^>j%XKnZA&eXee@VE-7&)3&I z0N%+l$ICIAnah^HavR$4pER_Q82zv$(ym--C|lb|hyGMcnQeHj{V~$sC8svhGOn@C zV?!&XTw*hpw0Ui9=d!7HYLlycyxJ7nha>U{2x0UAE`{+j0`u$P`qDhSPVIBmJ&4=_ zvJ~l~u=yS^4e2u&`5szpS%b}wJ~|*~X_v6{TZo|Q7|p8A*Vg!_V>C6~iN{)pOuM9| z&gr$WERjz0cj9-V4SyF_J~+^f|8*O_eC@6CuxsOpd;+4OfC~fukMRwd&y`GJ3KNDY zBYg@}$bl&VnZguuU`jxyFohiWn{o5HUHUOyDvu18sqp-(DNJF^_+Ng<<~XIen3(_o N002ovPDHLkV1o1n77PFY diff --git a/docs/source/_images/openmc_logo.png b/docs/source/_images/openmc_logo.png new file mode 100644 index 0000000000000000000000000000000000000000..73d4765387f00dbd0506f0edd38bffec474ac304 GIT binary patch literal 15452 zcmX|I1z6Ne*I$$d=~!Z=Q$Ubfx|Z%xX;`|uk2frT{S$dm~dySzD;NshTO&#MAMosUQF#uqz7@)9_e4 z$Wl}z`+3<^N0`0gj!380EY2i_RkVM45Wn>Q1M0l?8CXjvN1v3QGtTHRg zud*;yV+teh#WQIufPDhz0qy(D?%ER)_pU|<|!`rO1# zgsxzxnbkM<(z{i*QZ#66rtLV4&uo*$XI;U?HUmt8NUT$v!$KGYg&xloS{d^YG(uWA z=|L$h1$DtvvFDzQ17_jV*AGeURv8zI2taaylfBQ4qXhBODCk~Xt+Z_CCp%4Ca@rS; zl%f(4N<8x!W+3cV`U)j|C%$-CX0@r-b}@?+^;}|7&!gh7JpRQl)*eb zUgj`H1ZY~ZTjKw;c+(PctY!5DJf(u0?L_NzZ>?qK{$tsal}d$=D&gZsAX89sTr>X5 z9<=Ayn8{KnNeD1d%86ziTz{RdeEV$EIL=lG+A%QZpbXxsOF(rYim! zA0<;$Hp;Pl@iLQqrpY+;nMBD#4)(_51q(8F{!CwxCn~>iQtCOD|krT9f(|}y%-3N_bYZM`DfvRJN-;JwSotV-KIzk ziSIx-tKa)@{C^36uc!U+)t?=`(`YdwHO?rd z9x0^KQHZ?kB*I5Vq^TnPIbTI>1KZD6T#oVx?IT4;a2DX82(#g zAaevn7Wj010{txlfh=d;xZHlXt64!=;5r|h8Q#7+J9FNmfV1ZYc|3W{9v!qi^@J%N zD6Q!q4A0q&C4P#P3Puyz347dl6HC!q6&Y$k&KOF9Mo5L+d_B|v+UGl&%hD$_`o~j5 zg)O!^D{))YJyIKl{!HD`Zsq#~Vc_%QHWY9(c3NlNgcuI6bB#seP>GjaB-UClGvY7< z0FBB24x@%i5j+^zMy3 zU)O5+uh^VA5q_S>V{@BPiIm5@`?kfL1kGIvu@E#hFd*4Y(4FrUYIhhmMX?bihTHb_ zvh&9jj(GT77|%_3&o9sw+(|WRRHAthX#3`Bnh+}fu)LdO9ol&wXESY(5XT{p z)G2O)AMX;6TbKR}K=TLF3LlAi_#W5xN5fvEViPMH9_05*_( zeqgW-=3P|ZTiwBLlqvu~XJ}*VK~2u9mBpt;Ogd$gO6JH*6{Rl3z#gYd^Bwm(;+ zP7%PHx(#AsewSe}L^Qvj$oC^ua;>2@sZX~Yq}E>WnJ23B+xUKu9lIus6tvpjv6*%^ zMgx#-KU;$6|ozaV=g2dDbKEKYA*kUKNzH_R>_PJj#D;ay05SLbr6Uib93QEg`8rY!4}CyA-ZNhKoJ zk+!YiWUR*|>8Z5;Tbg%SM@!7tE|96g%!yK-DxQ%Er5smSFg>0Sd??GNWbZVyc5I79 z1~layhXleHC%BwCAKynBvw?&H*r7Ow3-_HeE-|%lL^P~rHDXS_H{t|sWl)@b&V#RR zkOftSHuV@TC70?Nmiw@9Jb@@FYnfYBI!n%>`yo_(w>%HMs~yO6aQ3gc7)BWI)9fO? zPbhpHR;!N{_>f{foM<1@2xxO45A78V&Xo(v`%b9dEUr8AqQB~d;t3oK8WJob>|lw2u|vsiDYW0=WzpR5UB$3b zx$K@4>eDE}%@qGM%=M5&5W}six*ixW^8#jeQ)GJW9bZ7HLFl2m&WTJbT=|em%73J*-6gwQ z)iM%Qp{NnG8`vN8D2pe`Bvh2;M6$=DN?%lfz2>4IC19nB>L~?9Gzcs>S}+U@Z_BPq zF}W#!H}KBV^}HUWO0CD%mUBltK&4QK#B6oNsPaI=1L9F%d1pzUFVhiz$n6QzT3IrtKXwpfK@-yjgN}IOBKO>DACgGMom^Bm=cb6LsbEh|7QcBRxPO}AMjo- zbHAbVm$NyF4ECuhwz(u40+gD=<%R6-4pDIknyc<^Q`Xy@2XbOexT68hm*_r#b(gYy zP4+&#ZmQW3Wih2yTR-~q6wV$`=}Y~`wAGxJKWKqEbp=W|vB{~PC`5DE2%ijwX|YqY zL1hP33JX>iC|Ga<|7>!rihK0W(?&uRKN$o0gHQp+^1Ov0a9>(X6RICxmP<*Fc*R^| zyM7?5bj$4+Unrqm`0o7J&F}mBdx0t%xG`qs-nFhMEiRe2Rs3e}s`w)KVy&D<7t_@R z&beFEfu#>=lbq+W*z*QRP-@jepHe@*x1b2U9ACL!71TU`OQld9FUtR^Huuq*OZ1B+1kaSqN+K$i(>Ek9+Q`;8HW9vuJdoYc z^={c6|E)PBStNZ}K{I+O|eyUhx zn)->8t6GB)os~I|b6u-_$)Kp?_K`wa`;yF*9C=789ET$wK|AX*Fn$M7b&3y9VGVrt z(elQ{fJQ|wMIRA~3G{WWF7vuI=f=f(%SRhkzVyP2If_t6cmMY<7{{ES@#E&0j{I<* z%E%{Q^VbQfR|~C&lE()@cU2b4^$t!&$2PVPS7zh^Z_aA1>595zo%hPB&&Yp^c<|)l zh#8r1M_cr{Cpy}`3p8SuV`-bCQrM(X(qLO2chJmE@SL$B&e)#}%d3K$z`@l0U!4_Q*ql@!w8|60_zmtJFxH@58+#KzGcQu;301PNK z1(bTf7t0rq)nJm$#z;krxTEJf0To8uP0RgWDyl?7NnPndGX=r735>w8pL|c)^?hzt z_1k4+rtm0O&^Y;L!%B-*7*{cZAR%)f6GQ{DLTtyyYB#ky{;|8&!ttmCDYgGY&ob6D z%bjyl!cu6TX1mDJW$r69d?o@kpIRYS_W&*SzFOp%7qMsBC}maX_|3A%{F)g+2F=GG zpHciT^U|hVW;S>+c{%3GOE6%0?y23d0{(Ah^~MjY#m3pq?0Lswzzd>$z1J5DB=OAE z)6O4(E4&nDnRAOfK>9ZdK)Hg1{9)4rz9WD38jX+s$o(ICs@~3*youJbSu)gGstafJ zLLxnuAT9>&r%lgjZi{@n>U$D)aq*gWsS}1dIMJOYO9+T-R;%0*w_Dv<%Wrg3hc2_H z)W2~2d_kbRwQ##Dw~I+&7QHhclT~q^!RlHc`|z+_sxIKho>gb*0@*z%QMCQadbzNA zRWxea%VKwrd4AE>fQJC~F;GM$AgO!1E*hn7D-_R`7Iw=Bmicfz1bupL|HV(ao4vV{_K)@HgKIg)`^fhUhmG;?D_-6RJ$@c&kl$$X zD#!9FAt&>b)#&2gI1UBt6F=%(NAK;_Q|4KxfqjVCvb5RiM)!GPD&t-v)sLJfRds2} zjL!E}D>Pgt8P@tfOX0Sod28qv^}&D6p^lx{&bV}DKC+4=PfJ45hLwSkSWMoV#}%@- z*fyh=#R+v=FX_lg4<*iuxfxV?8}im_z_Jb5J?Jk&z=J#>C8|c?o z#tzNGPg-%VXYyTVcMbcNny)63A1LqEdAi%)6;0gjMKPXsRk0Aw=LD(rig?=!2!Lom zJaUt1t<*ztjp z4)9CaQ7$`(yr&0)P?PfaZHb-?Y4_8EQHg4ueVbz)kkxTnQq?G;@$r7@LEn9y);eZ3M&@5L=I~ZQG0iOx6@(PnG)G*`Yg18((HqC0YlR2yHzp~>qFQT z$?(VB17cd+&u`Td8!#;SBOz!hIjJ5Bi%QoG$^tkh#&_m{k^<{xbQW_AQb^$o&GY}h-f>BHta!EU4hUkC|<=+DC;1T=0h;~c+hcxHxXT6 zq&+JBAxFftaJ|V+onB|gCv#Tj*@Xtg<+y-420=L!vA?~G{Pc1P3$hb=Gk0y>Xfl3m zxOBgc(}pr6O}(+7xi4LvokSwZ$7J~nl|51~UN~yG4Zwjsx`H!tRCz#O-kXMBw}_x z)`ar+uRqcbC1N04`G^Y=Ci1h1L@euL_xO2sIp8s(^c8XOdN9+Se`#_<+NK9;5%v3z zk=ZdayQl)$8Tb8*zKQj<7d5o?LR{}95VoeQ2Y&>7k6L;ZugLcr9N+!}MjVZ--%^ny z5dF#ZcI*CLcQ%nsqNzSXlfNG|=A%o@uB_dlA-TOI99f_Q~a5Ep9d z25d>B;Yh>d${H>#+;P+bodEcZa9m`ORK!e8vH^$F7Aw!utJ!wdbT6h8s`9sUe>>Iz zj-4_;1vHW^^S<*fb6cx3c^d_S6!6xmirqyTEL^5oEoG2}ijdu=*xHice+6M))yYREw}hCyq>G1u%Us7@(`5Tp zQVJRpFFJ{wS5!F^@}bOR-J}-NJwZu#KfmLPN|Q)?g!b(6o9%Ru(Onww*!&sC`w;m` zn-5I(R;ye7XhEEI&b4b7CHVj`|0gvKF20Qd$Q!NRX9q41iE1j~GWx#Q`q*$z;czh# zNIs4uu$W#$sPXMB$)#DQRu~`?Us9^0TTZV3SA+CYq&AKs1%8g7yxR+FBNH6@^*Rcs z!de%CuHgX6zC2zjK|) z5KXp8nkJ?PBRart<~NZba+QcGa%r_BJLa@Rum z-tzaov5XG z%8PZAVnrL;@E~(Jpw_ZdURUWe=Wdu!9TJYe?qW-J8(lP1;%{Uge`UVBJ)H~Y_{D6v-Uhf zIYAL5@+M?#9{87h%(u^&B@jZwY$ci{zVK4+FqzV?kwryjx(8Y+!|*WjAzy0r{@98o zRu?Hk4JiT^@5!e3sEUt$G#6gH7j0g1>>F(+*ZX1g%kS#N61A${K!0~=ukrw8NcU%- z6UH$T`T`L2g-y*W#Nz-j9Knp^iZI<^9tAvkLsfDrwWqK-pSGcj1XB$wVM6njEm_x3tEdy>j_D*D;?zoL86*DRrGecDk;0fqaPwzHl*^B*+xcaUU{* ziZsCz#1GTO#Zt0Asey|%JyLM!FvUJ%WhB}*9Z#;b(Me$2V0n^oARj^BbaO1*swH}? zQU1IN->_GlnGVXt6?{*V!}yfYLH@ec^#TXDLPR2aCs6%py{rDx)MC;gZ0R-W9C z?}s1Pt+%u53I+;ak%K}@hGy>iMD5G2V8tG=bsJ&)m+E@#!`Ck{$n1(sTx!Z#e?gYS zt%GDdo^lTbKivrX!Or~5Ty(_`J{AS6I^t1jpO$Rw+l1ruHDEnV_uTD0wNDP{wl%Dn z5f>HbU6z1W$DHNXaBk;iPw;}#T0YH^_ct;Dy}|(#c%Js(_<(~9inHXJlT@FNML57b zVHx4d&2zn@Ya-Efa>K(_;rMv9ug|!AbqYSJN%b9%g{oJWjD2>jWmMXvMtlk^bOE=- z8!IO0prg3qwW48|g&~DGfj6SUfhY!s9oCP6Ax;#tJup_+(z&82 z+w*;2OSqa#1CvA;;TD(UI#^koE^9#iDBMxe5enU7fuC z6+=C0YlVwn6mT^&c*Zs0<2aib^NYq?`NCb?nZ3PFG*u?30!Wr^fYV>uLKv5)XB7&- zbeuKIzNNR&h)Pm_66=2Qm3du0gXL*z*EGvX+IJh?Woet|IO#?Qa)e9CrPJklS1fFcRhZWR}^>eU(_RkZZ=U8b7*O{d+gy3!QLwv9F2sI1*G+h(VduPU!wnGaFquWTR9Z zO4tJ*52|svC8N>8M}k*zL+Oxe9%A^e%lKP-Ha67AT1yhAetf1oO^gX7iRymuhkZAN z18=U=L39W3t@_V#bMsU!@G_NAKFPepnY*Un+Wq}6asK{aSA1&w`_tPb(DjS$HXaXU zhysLqA~A-u&HhsNU~KeEi6Y7;lvqC=44|vqDw&+4+?rmWaPTzagX_+^>w{JZN`HWX4RVaOt(iwZ z4sCQjw$LqJyarF7B-2PlU-$z%4*?ki%>&)b3hLS&39SRGVM zAoxpgjd};YLrX|;t_$v)j+QYUcQn&Ik=u=brHlO2`QI8*6|iil^PTD;W4M zpDKyp5o(qqYBw%D>cUBK7Bh<5T%Ud2S0?f)zS_T!m!8H4{_OF>cwc6{?o=vXE8L9> zR_ANeC2CBM+ssM!D00&fU7nR+(@b@EM;-%-W=#{HQC4h&umP0xAyi~x1XvlD2#qDUUO;xOiCd$R zUh7fd?PeA>PgqnO7Q?8_*Lp{uU*<;-9%AYlAn<3U?KqO71q{4`9H5d-`=-8 zsR#(5WFT?k<-G1QPDVC}y3C`Iz8h~Wd-etCire6^IBB%3C zSR%H>dyxHSVTqHD6{{GRhlu&w)KfY?;{5A@lA0%`j(Ksq&5?W=66>Xuo`#!lukkFU z+qx6*GbI^JO~*tYxh-L%^~RN+7L?6zj=M|ZDo2UQkxy^(L&S84 zY&$PtWHg;d>)psuf%s@mp3WCMvX)33d6q2l6u)@{5x z0#cpU{+<%%s*!}>kLw0w$usFe6boC@gm1T*=z>_jZbw}q0}QDogBqEc$iG5=>d&|D z9do$H?~H)ylEp14bL;TRd>FZ3q!%$!Ioce!$E58QfmG!u@XE6A{(P`$pt!l)m)R#S zY*+`J8Xc6860%}5TNZnI4UK#7^34HvqtKx|mb+!`433c* ztFOoCg17^l9zqs{&p+p3D@`5m4$L1fyA`^$MxgV|_R;r#UmPuG>+ZZJm9!Gorlb{X zL{!=Ze(tK_A;RXW7Sh>TJ`j<~LbK0=EjIOzoAj&8d%mb&i;1;0uJ0UI0&MY8kbC?~+`$CJU`Ht*I$QaL5y zWbifRWiNjh8xuJcD4^37eA0Wz@m`&2x-36S9x_fAOrD{?q=XHT>`F2Qniw63u-}V) zo2%~T*%~G_Vo@mNJ*AbyG^ykAux1kYG^y{&d%7*Y5O;c(1KDG8h-G>pAmMoyK(psh|lz;SZH)b|5$x{ zLH~tKfM`!pN4JE6C@Zm~TA1GViD(j>%spzHE&OU;4)B}M`#pK(dkp_ooIMJHNaOF( ze)pNpYgIK*8~tfde5K$C?_R2#Ba5?CP8|+?uaTA(=rym|uhS8+c}Gz6ftpyKksey! zuwghhK5-k%@Q>PKGX-cMQ4ipRs2BYS-X|mT-(!)FtqrpND&^zS)OJvGtmz=*GI$|p zQds2hPrh6>-{D=a86oFoO|;@hlI}VZNSe*-c+V4J8`N|aCks_u&<%qLQ-s8u8!vEw ziSURGD|7~QT9;-puQeV^!>h{mAI^>20W2r%pNxy&(MvE?Ydf&(_$Mi8-gZ;)TI0E zoT-9_TiQvoM4jkb!2DyP$3S_sb)riXV80n@?Mi=2ETt&Iql!)k9Wb1idiE6GH9<)l z;$(SER7cEEpB;P+BDS7av5v%V<9k|7{YWEV$0aBditw%X#J;EC!<(2t=bBRFmjZ_x z@q~}OApVwJvNA!6eTQ7*oGwo4jE*F?;jiegI0`tMxF3)Jzqz+P^#9od9bRK4<2q6g z*?}fsObbzjS1dKV+B)W&PaFE1F;WZt>c@QOi!%19iO`zGNExKQW>pjTRZOmQ!)lOa ztNJ#1eqm6+Tu9NP82E+o3-@ets^q8~3g86Uz7W1AoEV*-d)yelwtfAr^YESVhF_%M z+3K$?6eDXGyb;Jjf;e-n@ZChQ^CAz+N}OJa6gdgJO_0hY()SO_`o-B&Cd z?@AV>b4)T!oN+@JT&2)@!$GX>4i0Qvw+)Y1#|{Hm2K4Xr@8ECWG<;za$2UrTL9V0c z-v`osL>P<<#v*>r1?1_Z*Wa&(tk)HH)52D~|76qjp*PC=SGUUSf}?F?U2P6JTDYi} zdJvx0AF}24_ljs`TIak-sPbufh#0&uGUIFNsx*EW4H9B}M)Dj7P2lBKzl-S(o(*Dk z1@949A%#m&UF~Hj?GfzNKyU?U9BmXc$M0&Ashg5J4GURM!M#Vq%n%Ju>pn@PT^m3s zxrt}hm@$;ma;)3o(L3Y`v>!Mybg{=Kzw)&MGHpr)f=lIe%1WwqFWxJRKCi zTd+#R7G1}3>Sct`6mh{yf;b$n>Ds#UMX7KPIW5U_SYw979oLuf{3ZreTiHz`P_aL) z7R6KKyhtG(G<)GF+&b~5!5j;ZOj@;uj>SO9tT#8O$Ep0nq|w1mTd;D_wGZhnF$`bX znRsJi;`lY0&lQhrur%^zR7pKZq$%W2m_zaJ-S!!CLI5Q#HM-}jI6Bl3Kav&~B|qY%ts6j)$)r2Ei zT}t*8*hDBvr4%S79WHU|*Q`~$Hq&Y9<3yN zr*=O(Ob#X)`;a>(y>=2BQG^V4f=&!MTvReqM{FO{Pv!@^Jj!Fj@9Q7|1DdNLIA_x& zagy%2AKCV|sTiu2Gj2tE;VxV}aVvC7t`)1@iZy*&#`RaA1j`3n^lcw`;+?>(F9x8K zU*^1;Xyz&TWDGpM-}RJC)*$dzrLkQ~a(FTO0~rq7lA1(zDoJ+k9LJ20G$cwi4JaVa z^B;^AE2+gUg+w2khRxPAN^#!nuTW^QckEvzLq*2N+Vxx5#%`HQ;;F04smIzOKf)c@ z(V0|AY4}|%r}s#ZwAE=Nb|=HZ9IiVMp@zH}NU+IQcEHtIZ8^ zTfD$D&KbE*m)*W(;XEg=Z@YZx+D!C&ahkC7V$Oxa$I#F{ie{HR{$)g>uW}L{4@@GO z1>K>x5l+kspgj*8lv`G6CCr8;Qh?}}^&XXiI_EkecRWI1UMak@J{=_q63(ymswkB-dJz(O3{$D?8;_(G_CtFk~yLCRq4@7 z{gCH7y&D)uF*Yfa+@Yqm*Afq2c5cELh?Mv{361d|@A3M(eC_p@lV2(X=_f|Ad!gIY zOrd9_Ww3%lm8R-D3OcKPFKZ2WFdO&6XPD1kp!Bw$(D~Cbo`9e&w)p4sW$w?j)O8Q5 z#8I@V6ha1qmSCWFbR|PsJyQxCzRg6NuE9^}li}<$k{OI!v>=#6T@?HRU^5&A zb{w%t-I^Sm`OV!?(zIrv`>ydcy0@`I<#w3gHVIkGV|=h{fn%!TWl0zPySH#C^musq zO;l-}pdl$~jxjY9fzq;P;UXn`p%8&|*}9&7;HwR-0`2VNn~xIfY+YB$m_8~InSPm- z_%Fl2&{t`(2>FOCzgNT~Vvk%QhTP`z1aMT*W*~n@A;Hr6l%T5k4$?K=9%hd-vk^po z;H9MZxP$IpO8wx2^cVFB_sz#NM;=pQYP5_GbL!V1**S8cwsaqsBQ+OqI36I5W~H4C zq}uYqrW{{{0;AwB4-x%zwUF_!jyXZuZ_}USC_btJGzUM=_Nw$|BXQn5_bK-9#Sh;o1dkiGSgK0S&tL9Y!l zZil<;eF}9rdmb5CpPhHT@a#XC3^cY3v{$4?v$FLG3ahFG#WN#RYA87x*|G~@6fDQC zdPWXT1LLPt;s^7PLlRh<(z& zI%y!jHbW8PwWmt~Vh|S`8f9sZEQNvWv+{pW(mWoyM1-&ploJLlEO> z9o9kHqaO1 z%|Bu$4-NZ$5#$WxBM_Um0EoNz<*e#KQ%+6MI)WFGZ?0rXkZ9(iUbiBwURWi$7PJTqDB15W--u!C=hiUN^#Ds~> z?QY@VGY0DyfKH>H8X^>(D{sVz-hn(t5f=gvnyz5}@$)_`3q!~tMiLu2K`UnWEBFDt zipEOg%JH54A7ZhSl>0Z$^rI3u z#|vHf`9?)3PVwj$l);_8;Xr-sZK~)v^Wy7sP6>)(6&7G z*xG^Yr!O+ED>C!3TTXt-+y+Xu-v4Qim-tk+{~|sH03h=hJ>7eJM@YL@+`g&TEhTmX z2W1C5Z3v%z=b4LFYUsoCVJG#032Xv==eC^=<_KW}cwq(d%kZX>s)|tGIy2s-5C34y zYg5e5s{PuFU1BRa^P0r$o1O07eeD=L2?>ITUL(*T0gbuq7kpjp=C0#A+(OhYL5_OE zQHWSB9c2LvkWfJHP4i20r0o_UeOHWsY;o=y9e-sau&szt(mj3VUJ_%smK}+w<8#T+ zD%1MT{vB3QC0S7=sBDowwwWhNlMDbzmy8tnK&q+_hivbTwSC~ocYP{ygo-QiUAUZ) z#a;|%O1%Tan*uBm)!t^ujIMn9*LUFO#XuEjOS3Ab?qZfnVT@Rj5hCAuV~&W{oBlpT zblf_6a5=_Ybk|r0frs$)K4sp8*la7Br7i$)m>>VLb~Od=y=BXuJCnprl{0RsT&e0f`=0rabVC=Le0}$-fJQxYiOg#TDj8)63O~wU%Nvbz)ezJL z9OlfQg-8*4KFCYLFuzn!-|%s5-xaf|Zk9*_UiF}f-?t8P6UX{yqiSR+i>?Jk@}mx^ z6O}R?4FX}zY6^#=1o@*SVTMUX14kmh_tz7PDL2sunEuAE{}h2$jhhrz*V||g=@Se1 za=A-E)6mQ|BRf)FD?vd$h@{|1tj!ij!QLZ!)=}=*FAd<*tJX^4QZdI8BSdRO&auth zyc-?2=u#J7SQH%3i>eq&F+88ef$c+Kb}-7YL}jKCOac|>^(ILKyHXMn5*^{ow2edvftI zXb_elR8}~@FuL8oBAe%&ay9Z=-o2-OhXB}o^P#mKiDyAtzfQmu`$qU-dJ*^&MuIk0n?y_h#BDB5U z?Kj|&XzRE8-tTRaIBT+4BAgB{nnC*IZ$CH&Rw8xaeQhHvL&=~@vbDtcC;tu(A4P6U zHK@$1V^HYE$MzK*Z0tDy(!sDQ>pgyft7>G}(9HU8HGfyH8hcS4BWfiG*bxpRxBxLc zWcV}wV*y+Sl@JDR0QhwJwF};wx4bmMr0f55ql%6VM!V={hU2nwoKW<)zOWAXj|nF{ zyC`F~6zo_QBDz$(d1`ahJITsH<;qN#PQr<>5I^|JE2s z6&m6#X&^>k?=c5oN-i0c9*T*a%~$pOQT5TGl>4_1fJm3s3(K?-k^$Os*~ z4A#eU=!Lh;AKj;M+ld*x`Dkrr_GF_-o(I8yYq!X)d!_8{L;ot!99CwRTe*rnI?&U5c%+u>*K|=datkyBTy62cK_5Mv?Klcj5n>= zP67CzP2h8DNEyP8pT_cgmDyGa!uyF! zqb=TNlp}r5Tgdh^NBY9{=l?2}q%2`iR)YNa1wuHa;}?TsJX`%XhJ~756z%_)Ifc9O z=ehu%V$^&3>)~WA2*RAmrhnRxLYzm__ufz0_lhU>N*Pkfj6?gZ`e0({ravnuCz27T z?h*nIgBCvv9)b%eWLRN(|4kJ@#<|l|0A~Y3p>w}Pq=659nNVpsM?H$`(V~$1FUJ2m z9?){${8z%9-(ALWm`6LZ@d(G@;%f5>sp=!J)$rT|nJ{-uvDfefLM%PrrJx8(s=(o$()NeZglFzJ_!}NcF4Ef; zBGv%P=ZP4HfQs)qUUZp@iU+u%I{vG;dHhczm%q&Y=|(@c6cx1c*4yWas)F75xaa?Q zs2~c3%DaDw=n99VWL4$=749!cWkmf??UzrdtW!DeyH%&~smw9;tocT~2eO!$_V5}O z6Xz1mu}k{=c^LqQn!|ZQh6;pAr(=LE<{axWY0rs>J4A&}wBd~?6?)ALbdadj0auMx zrBu~tWk+y%=`!Bq(-4@q-}{(>|5D0Q1HxfFn@ n;Y{RxR;H$UCQ;Hp8TcarI9P3v3rtRY{`n{?sRXGMHw^qgjvZG= literal 0 HcmV?d00001 diff --git a/docs/source/_static/theme_overrides.css b/docs/source/_static/theme_overrides.css index bee03f4150..dea941814d 100644 --- a/docs/source/_static/theme_overrides.css +++ b/docs/source/_static/theme_overrides.css @@ -16,3 +16,7 @@ .wy-table, .rst-content table.docutils, .rst-content table.field-list { margin-bottom: 0px; } + +.wy-side-nav-search { + background-color: #343131; +} diff --git a/docs/source/conf.py b/docs/source/conf.py index 1baea2b03c..75e621bc17 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -129,7 +129,7 @@ if not on_rtd: html_theme = 'sphinx_rtd_theme' html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] -html_logo = '_images/openmc200px.png' +html_logo = '_images/openmc_logo.png' # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". From 08219585f9328968e3ed49e42b89bb2eb3d4796a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 20 Aug 2016 16:50:14 -0500 Subject: [PATCH 06/23] Fix two documentation issues --- .../pythonapi/examples/mdgxs-part-i.rst | 13 ++++++++++++ .../pythonapi/examples/mdgxs-part-ii.rst | 13 ++++++++++++ openmc/mesh.py | 21 ++++++++++++------- 3 files changed, 39 insertions(+), 8 deletions(-) create mode 100644 docs/source/pythonapi/examples/mdgxs-part-i.rst create mode 100644 docs/source/pythonapi/examples/mdgxs-part-ii.rst diff --git a/docs/source/pythonapi/examples/mdgxs-part-i.rst b/docs/source/pythonapi/examples/mdgxs-part-i.rst new file mode 100644 index 0000000000..953dcf4700 --- /dev/null +++ b/docs/source/pythonapi/examples/mdgxs-part-i.rst @@ -0,0 +1,13 @@ +.. _notebook_mdgxs_part_i: + +========================== +MDGXS Part I: Introduction +========================== + +.. only:: html + + .. notebook:: mdgxs-part-i.ipynb + +.. only:: latex + + IPython notebooks must be viewed in the online HTML documentation. diff --git a/docs/source/pythonapi/examples/mdgxs-part-ii.rst b/docs/source/pythonapi/examples/mdgxs-part-ii.rst new file mode 100644 index 0000000000..a42eb766b7 --- /dev/null +++ b/docs/source/pythonapi/examples/mdgxs-part-ii.rst @@ -0,0 +1,13 @@ +.. _notebook_mdgxs_part_ii: + +================================ +MDGXS Part II: Advanced Features +================================ + +.. only:: html + + .. notebook:: mdgxs-part-ii.ipynb + +.. only:: latex + + IPython notebooks must be viewed in the online HTML documentation. diff --git a/openmc/mesh.py b/openmc/mesh.py index 58b9c7c0e5..7d7b483f73 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -187,15 +187,20 @@ class Mesh(object): of the mesh. For example the following code: - for mesh_index in mymesh.cell_generator(): - print mesh_index - will produce the following output for a 3-D 2x2x2 mesh in mymesh: - [1, 1, 1] - [1, 1, 2] - [1, 2, 1] - [1, 2, 2] - ... + .. code-block:: python + + for mesh_index in mymesh.cell_generator(): + print mesh_index + + will produce the following output for a 3-D 2x2x2 mesh in mymesh:: + + [1, 1, 1] + [1, 1, 2] + [1, 2, 1] + [1, 2, 2] + ... + """ From aaa80eb98926a271b59df6b992a90c48df88d637 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 23 Aug 2016 06:26:20 -0500 Subject: [PATCH 07/23] 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 9e23f11d88..29950a0bc2 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 08/23] 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 29950a0bc2..dc007aab61 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 09/23] 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 0000000000..a7352b79ab --- /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 10/23] 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 dc007aab61..ee5aeeaf1c 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 11/23] 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 33c955e3be..d0d67b6878 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 0df6042f6a..aed053263e 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 12/23] 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 d0d67b6878..83be99a33a 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 13/23] 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 83be99a33a..b75ffd88fb 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 14/23] 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 b75ffd88fb..e8ed04e3a0 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 15/23] 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 4bbc96cf86..2ab823e6ed 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 73b51da5e8..4dacabf288 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 16/23] 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 4dacabf288..cc5c0d44b3 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 17/23] 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 e8ed04e3a0..76293051b0 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 18/23] 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 76293051b0..0c372aeefd 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 19/23] 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 0c372aeefd..5525ea559c 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 20/23] 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 de7475308c..d0bc07987c 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 21/23] 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 d0bc07987c..dd5c2056c6 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 22/23] 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 dd5c2056c6..543caf088a 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 23/23] 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 543caf088a..1d2e9098d3 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':