diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 4b1d8cca7..09b3d5135 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -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 diff --git a/openmc/arithmetic.py b/openmc/arithmetic.py index 5869cad80..f9ef58db8 100644 --- a/openmc/arithmetic.py +++ b/openmc/arithmetic.py @@ -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 diff --git a/openmc/filter.py b/openmc/filter.py index b34400cbb..72fb3b14e 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -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) diff --git a/openmc/material.py b/openmc/material.py index 6b2b07f2d..66c2320ea 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -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") diff --git a/openmc/mesh.py b/openmc/mesh.py index 8bad6c537..0c88d4a68 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -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 diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py new file mode 100644 index 000000000..ffb1f4282 --- /dev/null +++ b/openmc/model/__init__.py @@ -0,0 +1 @@ +from .triso import * diff --git a/openmc/model/triso.py b/openmc/model/triso.py new file mode 100644 index 000000000..89e0d8aa7 --- /dev/null +++ b/openmc/model/triso.py @@ -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 diff --git a/openmc/tallies.py b/openmc/tallies.py index 38aa0bf12..af19549a7 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -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 diff --git a/openmc/trigger.py b/openmc/trigger.py index b8383bd27..537af1c8d 100644 --- a/openmc/trigger.py +++ b/openmc/trigger.py @@ -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 diff --git a/openmc/universe.py b/openmc/universe.py index 7f32d68a9..84a5acdca 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -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. diff --git a/setup.py b/setup.py index 770f280ad..99e465c9d 100644 --- a/setup.py +++ b/setup.py @@ -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 diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat new file mode 100644 index 000000000..3e54e5e6f --- /dev/null +++ b/tests/test_triso/inputs_true.dat @@ -0,0 +1 @@ +2dcfd1a17cba671874e60192a7355deb57e2e51467a474fd168c8b51e454a977edb34df07ae11625c0a43906112152c75113e442a9a8f240a4c9d1a11ee4771d \ No newline at end of file diff --git a/tests/test_triso/plots.xml b/tests/test_triso/plots.xml new file mode 100644 index 000000000..60ae7d9d8 --- /dev/null +++ b/tests/test_triso/plots.xml @@ -0,0 +1,6 @@ + + + + + diff --git a/tests/test_triso/results_true.dat b/tests/test_triso/results_true.dat new file mode 100644 index 000000000..ea7da21ed --- /dev/null +++ b/tests/test_triso/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.662675E+00 1.475968E-02 diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py new file mode 100644 index 000000000..d1ac4e5cd --- /dev/null +++ b/tests/test_triso/test_triso.py @@ -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()