From f3cfc181f47d959158d5e871209c45b9aa810f33 Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 16 Aug 2016 20:14:34 -0500 Subject: [PATCH] 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')