Merge pull request #667 from paulromano/various-pyapi-changes

Various changes in Python API
This commit is contained in:
Sterling Harper 2016-06-09 10:21:32 -04:00
commit 47e992b889
15 changed files with 395 additions and 306 deletions

View file

@ -299,6 +299,33 @@ Multi-group Cross Section Libraries
openmc.mgxs.Library
-------------------------------------
:mod:`openmc.model` -- Model Building
-------------------------------------
TRISO Fuel Modeling
-------------------
Classes
+++++++
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.model.TRISO
Functions
+++++++++
.. autosummary::
:toctree: generated
:nosignatures:
openmc.model.create_triso_lattice
.. _Jupyter: https://jupyter.org/
.. _NumPy: http://www.numpy.org/
.. _Codecademy: https://www.codecademy.com/tracks/python

View file

@ -67,24 +67,6 @@ class CrossScore(object):
def __ne__(self, other):
return not self == other
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._left_score = self.left_score
clone._right_score = self.right_score
clone._binary_op = self.binary_op
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
string = '({0} {1} {2})'.format(self.left_score,
self.binary_op, self.right_score)
@ -169,28 +151,9 @@ class CrossNuclide(object):
def __ne__(self, other):
return not self == other
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._left_nuclide = self.left_nuclide
clone._right_nuclide = self.right_nuclide
clone._binary_op = self.binary_op
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
return self.name
@property
def left_nuclide(self):
return self._left_nuclide
@ -325,27 +288,6 @@ class CrossFilter(object):
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', filter_bins)
return string
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._left_filter = self.left_filter
clone._right_filter = self.right_filter
clone._binary_op = self.binary_op
clone._type = self.type
clone._bins = self._bins
clone._stride = self.stride
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def left_filter(self):
return self._left_filter
@ -532,23 +474,6 @@ class AggregateScore(object):
def __ne__(self, other):
return not self == other
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._scores = self.scores
clone._aggregate_op = self.aggregate_op
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
string = ', '.join(map(str, self.scores))
string = '{0}({1})'.format(self.aggregate_op, string)
@ -622,23 +547,6 @@ class AggregateNuclide(object):
def __ne__(self, other):
return not self == other
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._nuclides = self.nuclides
clone._aggregate_op = self._aggregate_op
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
# Append each nuclide in the aggregate to the string
@ -757,26 +665,6 @@ class AggregateFilter(object):
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
return string
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._type = self.type
clone._aggregate_filter = self.aggregate_filter
clone._aggregate_op = self.aggregate_op
clone._bins = self._bins
clone._stride = self.stride
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def aggregate_filter(self):
return self._aggregate_filter

View file

@ -104,27 +104,6 @@ class Filter(object):
def __hash__(self):
return hash(repr(self))
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._type = self.type
clone._bins = copy.deepcopy(self.bins, memo)
clone._num_bins = self.num_bins
clone._mesh = copy.deepcopy(self.mesh, memo)
clone._stride = self.stride
clone._distribcell_paths = copy.deepcopy(self.distribcell_paths)
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
string = 'Filter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)

View file

@ -50,14 +50,14 @@ class Material(object):
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/cm3',
'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only
applies in the case of a multi-group calculation.
elements : collections.OrderedDict
Dictionary whose keys are element names and values are 3-tuples
consisting of an :class:`openmc.Element` instance, the percent density,
and the percent type (atom or weight fraction).
nuclides : collections.OrderedDict
Dictionary whose keys are nuclide names and values are 3-tuples
consisting of an :class:`openmc.Nuclide` instance, the percent density,
and the percent type (atom or weight fraction).
elements : list of tuple
List in which each item is a 3-tuple consisting of an
:class:`openmc.Element` instance, the percent density, and the percent
type ('ao' or 'wo').
nuclides : list of tuple
List in which each item is a 3-tuple consisting of an
:class:`openmc.Nuclide` instance, the percent density, and the percent
type ('ao' or 'wo').
"""
@ -68,19 +68,15 @@ class Material(object):
self._density = None
self._density_units = ''
# An ordered dictionary of Nuclides (order affects OpenMC results)
# Keys - Nuclide names
# Values - tuple (nuclide, percent, percent type)
self._nuclides = OrderedDict()
# A list of tuples (nuclide, percent, percent type)
self._nuclides = []
# The single instance of Macroscopic data present in this material
# (only one is allowed, hence this is different than _nuclides, etc)
self._macroscopic = None
# An ordered dictionary of Elements (order affects OpenMC results)
# Keys - Element names
# Values - tuple (element, percent, percent type)
self._elements = OrderedDict()
# A list of tuples (element, percent, percent type)
self._elements = []
# If specified, a list of tuples of (table name, xs identifier)
self._sab = []
@ -134,10 +130,8 @@ class Material(object):
string += '{0: <16}\n'.format('\tNuclides')
for nuclide in self._nuclides:
percent = self._nuclides[nuclide][1]
percent_type = self._nuclides[nuclide][2]
string += '{0: <16}'.format('\t{0}'.format(nuclide))
for nuclide, percent, percent_type in self._nuclides:
string += '{0: <16}'.format('\t{0.name}.{0.xs}'.format(nuclide))
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
if self._macroscopic is not None:
@ -146,39 +140,12 @@ class Material(object):
string += '{0: <16}\n'.format('\tElements')
for element in self._elements:
percent = self._elements[element][1]
percent_type = self._elements[element][2]
string += '{0: <16}'.format('\t{0}'.format(element))
for element, percent, percent_type in self._elements:
string += '{0: <16}'.format('\t{0.name}.{0.xs}'.format(element))
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
return string
def __deepcopy__(self, memo):
existing = memo.get(id(self))
if existing is None:
# If this is the first time we have tried to copy this object, create a copy
clone = type(self).__new__(type(self))
clone._id = self._id
clone._name = self._name
clone._density = self._density
clone._density_units = self._density_units
clone._nuclides = deepcopy(self._nuclides, memo)
clone._macroscopic = self._macroscopic
clone._elements = deepcopy(self._elements, memo)
clone._sab = deepcopy(self._sab, memo)
clone._convert_to_distrib_comps = self._convert_to_distrib_comps
clone._distrib_otf_file = self._distrib_otf_file
memo[id(self)] = clone
return clone
else:
# If this object has been copied before, return the first copy made
return existing
@property
def id(self):
return self._id
@ -326,7 +293,7 @@ class Material(object):
else:
nuclide = openmc.Nuclide(nuclide)
self._nuclides[nuclide._name] = (nuclide, percent, percent_type)
self._nuclides.append((nuclide, percent, percent_type))
def remove_nuclide(self, nuclide):
"""Remove a nuclide from the material
@ -344,8 +311,9 @@ class Material(object):
raise ValueError(msg)
# If the Material contains the Nuclide, delete it
if nuclide._name in self._nuclides:
del self._nuclides[nuclide._name]
for nuc in self._nuclides:
if nuclide == nuc:
self._nuclides.remove(nuc)
def add_macroscopic(self, macroscopic):
"""Add a macroscopic to the material. This will also set the
@ -464,10 +432,9 @@ class Material(object):
raise NotImplementedError('Expanding natural element based on '
'weight percent is not yet supported.')
for isotope, abundance in element.expand():
self._nuclides[isotope.name] = (
isotope, percent*abundance, percent_type)
self._nuclides.append((isotope, percent*abundance, percent_type))
else:
self._elements[element.name] = (element, percent, percent_type)
self._elements.append((element, percent, percent_type))
def remove_element(self, element):
"""Remove a natural element from the material
@ -479,9 +446,15 @@ class Material(object):
"""
# If the Material contains the Element, delete it
if element._name in self._elements:
del self._elements[element._name]
if not isinstance(nuclide, openmc.Element):
msg = 'Unable to remove "{0}" in Material ID="{1}" ' \
'since it is not an Element'.format(self.id, element)
raise ValueError(msg)
# If the Material contains the Nuclide, delete it
for elm in self._elements:
if element == elm:
self._nuclides.remove(elm)
def add_s_alpha_beta(self, name, xs):
r"""Add an :math:`S(\alpha,\beta)` table to the material
@ -513,10 +486,10 @@ class Material(object):
self._sab.append((name, xs))
def make_isotropic_in_lab(self):
for nuclide_name in self._nuclides:
self._nuclides[nuclide_name][0].scattering = 'iso-in-lab'
for element_name in self._elements:
self._elements[element_name][0].scattering = 'iso-in-lab'
for nuclide, percent, percent_type in self._nuclides:
nuclide.scattering = 'iso-in-lab'
for element, percent, percent_type in self._elements:
element.scattering = 'iso-in-lab'
def get_all_nuclides(self):
"""Returns all nuclides in the material
@ -531,15 +504,10 @@ class Material(object):
nuclides = OrderedDict()
for nuclide_name, nuclide_tuple in self._nuclides.items():
nuclide = nuclide_tuple[0]
density = nuclide_tuple[1]
nuclides[nuclide._name] = (nuclide, density)
for element_name, element_tuple in self._elements.items():
element = element_tuple[0]
density = element_tuple[1]
for nuclide, density, density_type in self._nuclides:
nuclides[nuclide.name] = (nuclide, density)
for element, density, density_type in self._elements:
# Expand natural element into isotopes
for isotope, abundance in element.expand():
nuclides[isotope.name] = (isotope, density*abundance)
@ -548,7 +516,7 @@ class Material(object):
def _get_nuclide_xml(self, nuclide, distrib=False):
xml_element = ET.Element("nuclide")
xml_element.set("name", nuclide[0]._name)
xml_element.set("name", nuclide[0].name)
if not distrib:
if nuclide[2] == 'ao':
@ -575,7 +543,7 @@ class Material(object):
def _get_element_xml(self, element, distrib=False):
xml_element = ET.Element("element")
xml_element.set("name", str(element[0]._name))
xml_element.set("name", str(element[0].name))
if not distrib:
if element[2] == 'ao':
@ -594,7 +562,7 @@ class Material(object):
def _get_nuclides_xml(self, nuclides, distrib=False):
xml_elements = []
for nuclide in nuclides.values():
for nuclide in nuclides:
xml_elements.append(self._get_nuclide_xml(nuclide, distrib))
return xml_elements
@ -602,7 +570,7 @@ class Material(object):
def _get_elements_xml(self, elements, distrib=False):
xml_elements = []
for element in elements.values():
for element in elements:
xml_elements.append(self._get_element_xml(element, distrib))
return xml_elements
@ -650,7 +618,7 @@ class Material(object):
subelement = ET.SubElement(element, "compositions")
comps = []
allnucs = self._nuclides.values() + self._elements.values()
allnucs = self._nuclides + self._elements
dist_per_type = allnucs[0][2]
for nuc, per, typ in allnucs:
if not typ == dist_per_type:
@ -845,4 +813,4 @@ class Materials(cv.CheckedList):
# Write the XML Tree to the materials.xml file
tree = ET.ElementTree(self._materials_file)
tree.write("materials.xml", xml_declaration=True,
encoding='utf-8', method="xml")
encoding='utf-8', method="xml")

View file

@ -1,5 +1,4 @@
from collections import Iterable
import copy
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
@ -83,28 +82,6 @@ class Mesh(object):
else:
return True
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._id = self._id
clone._name = self._name
clone._type = self._type
clone._dimension = copy.deepcopy(self._dimension, memo)
clone._lower_left = copy.deepcopy(self._lower_left, memo)
clone._upper_right = copy.deepcopy(self._upper_right, memo)
clone._width = copy.deepcopy(self._width, memo)
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def id(self):
return self._id

1
openmc/model/__init__.py Normal file
View file

@ -0,0 +1 @@
from .triso import *

155
openmc/model/triso.py Normal file
View file

@ -0,0 +1,155 @@
import copy
from collections import Iterable
from numbers import Real
import warnings
import numpy as np
import openmc
import openmc.checkvalue as cv
class TRISO(openmc.Cell):
"""Tristructural-isotopic (TRISO) micro fuel particle
Parameters
----------
outer_radius : float
Outer radius of TRISO particle
fill : openmc.Universe
Universe which contains all layers of the TRISO particle
center : Iterable of float
Cartesian coordinates of the center of the TRISO particle in cm
Attributes
----------
id : int
Unique identifier for the TRISO cell
name : str
Name of the TRISO cell
center : numpy.ndarray
Cartesian coordinates of the center of the TRISO particle in cm
fill : openmc.Universe
Universe that contains the TRISO layers
region : openmc.Region
Region of space within the TRISO particle
"""
def __init__(self, outer_radius, fill, center=(0., 0., 0.)):
self._surface = openmc.Sphere(R=outer_radius)
super(TRISO, self).__init__(fill=fill, region=-self._surface)
self.center = np.asarray(center)
@property
def center(self):
return self._center
@center.setter
def center(self, center):
cv.check_type('TRISO center', center, Iterable, Real)
self._surface.x0 = center[0]
self._surface.y0 = center[1]
self._surface.z0 = center[2]
self.translation = center
self._center = center
def classify(self, lattice):
"""Determine lattice element indices which might contain the TRISO particle.
Parameters
----------
lattice : openmc.RectLattice
Lattice to check
Returns
-------
list of tuple
(z,y,x) lattice element indices which might contain the TRISO
particle.
"""
ll, ur = self.region.bounding_box
if lattice.ndim == 2:
(i_min, j_min), p = lattice.find_element(ll)
(i_max, j_max), p = lattice.find_element(ur)
return list(np.broadcast(*np.ogrid[
j_min:j_max+1, i_min:i_max+1]))
else:
(i_min, j_min, k_min), p = lattice.find_element(ll)
(i_max, j_max, k_max), p = lattice.find_element(ur)
return list(np.broadcast(*np.ogrid[
k_min:k_max+1, j_min:j_max+1, i_min:i_max+1]))
def create_triso_lattice(trisos, lower_left, pitch, shape, background):
"""Create a lattice containing TRISO particles for optimized tracking.
Parameters
----------
trisos : list of openmc.model.TRISO
List of TRISO particles to put in lattice
lower_left : Iterable of float
Lower-left Cartesian coordinates of the lattice
pitch : Iterable of float
Pitch of the lattice elements in the x-, y-, and z-directions
shape : Iterable of float
Number of lattice elements in the x-, y-, and z-directions
background : openmc.Material
A background material that is used anywhere within the lattice but
outside a TRISO particle
Returns
-------
lattice : openmc.RectLattice
A lattice containing the TRISO particles
"""
lattice = openmc.RectLattice()
lattice.lower_left = lower_left
lattice.pitch = pitch
indices = list(np.broadcast(*np.ogrid[:shape[2], :shape[1], :shape[0]]))
triso_locations = {idx: [] for idx in indices}
for t in trisos:
for idx in t.classify(lattice):
if idx in sorted(triso_locations):
# Create copy of TRISO particle with materials preserved and
# different cell/surface IDs
t_copy = copy.deepcopy(t)
t_copy.id = None
t_copy.fill = t.fill
t_copy._surface.id = None
triso_locations[idx].append(t_copy)
else:
warnings.warn('TRISO particle is partially or completely '
'outside of the lattice.')
# Create universes
universes = np.empty(shape[::-1], dtype=openmc.Universe)
for idx, triso_list in sorted(triso_locations.items()):
if len(triso_list) > 0:
outside_trisos = openmc.Intersection(*[~t.region for t in triso_list])
background_cell = openmc.Cell(fill=background, region=outside_trisos)
else:
background_cell = openmc.Cell(fill=background)
u = openmc.Universe()
u.add_cell(background_cell)
for t in triso_list:
u.add_cell(t)
iz, iy, ix = idx
t.center = lattice.get_local_coordinates(t.center, (ix, iy, iz))
if len(shape) == 2:
universes[-1 - idx[0], idx[1]] = u
else:
universes[idx[0], -1 - idx[1], idx[2]] = u
lattice.universes = universes
# Set outer universe
background_cell = openmc.Cell(fill=background)
lattice.outer = openmc.Universe(cells=[background_cell])
return lattice

View file

@ -129,51 +129,6 @@ class Tally(object):
self._sp_filename = None
self._results_read = False
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone.id = self.id
clone.name = self.name
clone.estimator = self.estimator
clone.num_realizations = self.num_realizations
clone._sum = copy.deepcopy(self._sum, memo)
clone._sum_sq = copy.deepcopy(self._sum_sq, memo)
clone._mean = copy.deepcopy(self._mean, memo)
clone._std_dev = copy.deepcopy(self._std_dev, memo)
clone._with_summary = self.with_summary
clone._with_batch_statistics = self.with_batch_statistics
clone._derived = self.derived
clone._sparse = self.sparse
clone._sp_filename = self._sp_filename
clone._results_read = self._results_read
clone._filters = []
for self_filter in self.filters:
clone.filters.append(copy.deepcopy(self_filter, memo))
clone._nuclides = []
for nuclide in self.nuclides:
clone.nuclides.append(copy.deepcopy(nuclide, memo))
clone._scores = []
for score in self.scores:
clone.scores.append(score)
clone._triggers = []
for trigger in self.triggers:
clone.triggers.append(trigger)
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __eq__(self, other):
if not isinstance(other, Tally):
return False
@ -2879,7 +2834,7 @@ class Tally(object):
return other * self**-1
def __pos__(self):
def __abs__(self):
"""The absolute value of this tally.
Returns

View file

@ -39,25 +39,6 @@ class Trigger(object):
self.threshold = threshold
self._scores = []
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._trigger_type = self._trigger_type
clone._threshold = self._threshold
clone.scores = self.scores
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __eq__(self, other):
if str(self) == str(other):
return True

View file

@ -155,7 +155,7 @@ class Universe(object):
return []
def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200),
basis='xy', color_by='cell', seed=None):
basis='xy', color_by='cell', colors=None, filename=None, seed=None):
"""Display a slice plot of the universe.
Parameters
@ -170,6 +170,22 @@ class Universe(object):
The basis directions for the plot
color_by : {'cell', 'material'}
Indicate whether the plot should be colored by cell or by material
colors : dict
Assigns colors to specific materials or cells. Keys are instances of
:class:`Cell` or :class:`Material` and values are RGB 3-tuples or RGBA
4-tuples. Red, green, blue, and alpha should all be floats in the
range [0.0, 1.0], for example:
.. code-block:: python
# Make water blue
water = openmc.Cell(fill=h2o)
universe.plot(..., colors={water: (0., 0., 1.))
filename : str or None
Filename to save plot to. If no filename is given, the plot will be
displayed using the currently enabled matplotlib backend.
seed : hashable object or None
Hashable object which is used to seed the random number generator
used to select colors. If None, the generator is seeded from the
@ -182,6 +198,15 @@ class Universe(object):
if seed is not None:
random.seed(seed)
if colors is None:
# Create default dictionary if none supplied
colors = {}
else:
# Convert to RGBA if necessary
for obj, rgb in colors.items():
if len(rgb) == 3:
colors[obj] = rgb + (1.0,)
if basis == 'xy':
x_min = center[0] - 0.5*width[0]
x_max = center[0] + 0.5*width[0]
@ -206,7 +231,7 @@ class Universe(object):
y_coords = np.linspace(y_max, y_min, pixels[1], endpoint=False) - \
0.5*(y_max - y_min)/pixels[1]
colors = {}
# Search for locations and assign colors
img = np.zeros(pixels + (4,)) # Use RGBA form
for i, x in enumerate(x_coords):
for j, y in enumerate(y_coords):
@ -220,21 +245,27 @@ class Universe(object):
if len(path) > 0:
try:
if color_by == 'cell':
uid = path[-1].id
obj = path[-1]
elif color_by == 'material':
if path[-1].fill_type == 'material':
uid = path[-1].fill.id
obj = path[-1].fill
else:
continue
except AttributeError:
continue
if uid not in colors:
colors[uid] = (random.random(), random.random(),
if obj not in colors:
colors[obj] = (random.random(), random.random(),
random.random(), 1.0)
img[j,i,:] = colors[uid]
img[j,i,:] = colors[obj]
# Display image
plt.imshow(img, extent=(x_min, x_max, y_min, y_max))
plt.show()
# Show or save the plot
if filename is None:
plt.show()
else:
plt.savefig(filename)
def add_cell(self, cell):
"""Add a cell to the universe.

View file

@ -11,7 +11,8 @@ except ImportError:
kwargs = {'name': 'openmc',
'version': '0.7.1',
'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.stats'],
'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.model',
'openmc.stats'],
'scripts': glob.glob('scripts/openmc-*'),
# Metadata

View file

@ -0,0 +1 @@
2dcfd1a17cba671874e60192a7355deb57e2e51467a474fd168c8b51e454a977edb34df07ae11625c0a43906112152c75113e442a9a8f240a4c9d1a11ee4771d

View file

@ -0,0 +1,6 @@
<?xml version="1.0"?>
<plots>
<plot id="1" type="slice" basis="xy" color="material"
origin="0.0 0.0 0.0" width="1.0 1.0" pixels="400 400">
</plot>
</plots>

View file

@ -0,0 +1,2 @@
k-combined:
1.662675E+00 1.475968E-02

View file

@ -0,0 +1,117 @@
#!/usr/bin/env python
import os
import sys
import glob
import random
from math import sqrt
import numpy as np
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
import openmc
import openmc.model
class TRISOTestHarness(PyAPITestHarness):
def _build_inputs(self):
# 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)
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')
ipyc = openmc.Material()
ipyc.set_density('g/cm3', 1.90)
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('C-Nat', 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')
graphite = openmc.Material()
graphite.set_density('g/cm3', 1.1995)
graphite.add_nuclide('C-Nat', 1.0)
graphite.add_s_alpha_beta('Graph', '71t')
# Create TRISO particles
spheres = [openmc.Sphere(R=r*1e-4)
for r in [212.5, 312.5, 347.5, 382.5]]
c1 = openmc.Cell(fill=fuel, region=-spheres[0])
c2 = openmc.Cell(fill=porous_carbon, region=+spheres[0] & -spheres[1])
c3 = openmc.Cell(fill=ipyc, region=+spheres[1] & -spheres[2])
c4 = openmc.Cell(fill=sic, region=+spheres[2] & -spheres[3])
c5 = openmc.Cell(fill=opyc, region=+spheres[3])
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)
# Define box to contain lattice
min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')
max_x = openmc.XPlane(x0=0.5, boundary_type='reflective')
min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective')
max_y = openmc.YPlane(y0=0.5, boundary_type='reflective')
min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective')
max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective')
box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z)
# Create lattice
ll, ur = box.region.bounding_box
shape = (3, 3, 3)
pitch = (ur - ll) / shape
lattice = openmc.model.create_triso_lattice(
trisos, ll, pitch, shape, graphite)
box.fill = lattice
root = openmc.Universe(0, cells=[box])
geom = openmc.Geometry(root)
geom.export_to_xml()
settings = openmc.Settings()
settings.batches = 5
settings.inactive = 0
settings.particles = 100
settings.source = openmc.Source(space=openmc.stats.Point())
settings.export_to_xml()
mats = openmc.Materials([fuel, porous_carbon, ipyc, sic, opyc, graphite])
mats.default_xs = '71c'
mats.export_to_xml()
if __name__ == '__main__':
harness = TRISOTestHarness('statepoint.5.h5')
harness.main()