From 108d8c835b58df505f02b124a1f0e9d13f4eb0d0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 31 May 2016 12:53:12 -0500 Subject: [PATCH 1/7] Add colors and filename argument to Universe.plot --- openmc/universe.py | 47 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 39 insertions(+), 8 deletions(-) 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. From d16e3fca8976ccdc11d1eaf5583992e38bf96a2c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 1 Jun 2016 01:40:17 -0500 Subject: [PATCH 2/7] Remove many __deepcopy__ implementations --- openmc/arithmetic.py | 112 ------------------------------------------- openmc/filter.py | 21 -------- openmc/material.py | 25 ---------- openmc/mesh.py | 23 --------- openmc/tallies.py | 47 +----------------- openmc/trigger.py | 19 -------- 6 files changed, 1 insertion(+), 246 deletions(-) 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 52560a193..e06958490 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..f29a42117 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -154,31 +154,6 @@ class Material(object): 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 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/tallies.py b/openmc/tallies.py index 2313a6173..1d895022c 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 @@ -2875,7 +2830,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 From 07ef67fee019b3c2eb580fc651e6bbbf0bf2e0cb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 26 May 2016 07:46:07 -0500 Subject: [PATCH 3/7] Initial TRISO modeling capabilities --- docs/source/pythonapi/index.rst | 27 +++++ openmc/model/__init__.py | 1 + openmc/model/triso.py | 172 ++++++++++++++++++++++++++++++ setup.py | 3 +- tests/test_triso/inputs_true.dat | 1 + tests/test_triso/results_true.dat | 2 + tests/test_triso/test_triso.py | 107 +++++++++++++++++++ 7 files changed, 312 insertions(+), 1 deletion(-) create mode 100644 openmc/model/__init__.py create mode 100644 openmc/model/triso.py create mode 100644 tests/test_triso/inputs_true.dat create mode 100644 tests/test_triso/results_true.dat create mode 100644 tests/test_triso/test_triso.py diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index bf35e7587..54cb93594 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -296,6 +296,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/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..dc53e6ec0 --- /dev/null +++ b/openmc/model/triso.py @@ -0,0 +1,172 @@ +import copy +from collections import Iterable +from numbers import Real + +import numpy as np + +import openmc +import openmc.checkvalue as cv + +class TRISO(object): + """Tristructural-isotopic (TRISO) micro fuel particle + + Parameters + ---------- + materials : Iterable of openmc.Material + Material to be assigned to each layer of the TRISO particle starting + with the innermost and proceeding outwards + radii : Iterable of float + Outer radii in cm of each layer of the TRISO particle in ascending order + center : Iterable of float + Cartesian coordinates of the center of the TRISO particle in cm + + Attributes + ---------- + cells : list of opemc.Cell + Each layer of the TRISO particle + center : numpy.ndarray + Cartesian coordinates of the center of the TRISO particle in cm + outside : openmc.Region + Region of space outside of the TRISO particle + bounding_box : tuple of numpy.ndarray + Lower-left and upper-right coordinates of an axis-aligned bounding box + for the TRISO particle + + """ + + def __init__(self, materials, radii, center=(0., 0., 0.)): + surfaces = [openmc.Sphere(R=r) for r in radii] + cells = [] + for i, m in enumerate(materials): + c = openmc.Cell(fill=m) + if i == 0: + c.region = -surfaces[i] + else: + c.region = +surfaces[i-1] & -surfaces[i] + cells.append(c) + self._cells = cells + self._surfaces = surfaces + self.center = np.asarray(center) + + @property + def bounding_box(self): + return self.cells[-1].region.bounding_box + + @property + def cells(self): + return self._cells + + @property + def center(self): + return self._center + + @property + def outside(self): + return ~self.cells[-1].region.nodes[-1] + + @center.setter + def center(self, center): + cv.check_type('TRISO center', center, Iterable, Real) + for s in self._surfaces: + s.x0, s.y0, s.z0 = 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.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) + for c, c_copy in zip(t.cells, t_copy.cells): + c_copy.id = None + c_copy.fill = c.fill + for s in t_copy._surfaces: + s.id = None + triso_locations[idx].append(t_copy) + + # 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.outside 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_cells(t.cells) + 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/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..40f312850 --- /dev/null +++ b/tests/test_triso/inputs_true.dat @@ -0,0 +1 @@ +6c6fecedf3db0b91b69b7b7b57b3a52c42947253bea4ee3889d6bbd6e74935cc4e599c1c743eb589a1045868412f3f88c5fef7cf10e180bdc4bd182970c9ee15 \ No newline at end of file diff --git a/tests/test_triso/results_true.dat b/tests/test_triso/results_true.dat new file mode 100644 index 000000000..8f968f65a --- /dev/null +++ b/tests/test_triso/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.685303E+00 1.121936E-01 diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py new file mode 100644 index 000000000..f7f4daa7b --- /dev/null +++ b/tests/test_triso/test_triso.py @@ -0,0 +1,107 @@ +#!/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 + materials = [fuel, porous_carbon, ipyc, sic, opyc] + radii = np.array([212.5, 312.5, 347.5, 382.5, 422.5])*1e-4 + trisos = [] + random.seed(1) + for i in range(100): + # Randomly sample location + x = random.uniform(-0.5, 0.5) + y = random.uniform(-0.5, 0.5) + z = random.uniform(-0.5, 0.5) + t = openmc.model.TRISO(materials, radii, (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*radii[-1]: + 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) + lattice = openmc.model.create_triso_lattice( + trisos, ll, (ur - ll)/shape, 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 = 50 + 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() From 55e1ebff2d06c4c70599ab6c3dbe07cb57a5b3ae Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 7 Jun 2016 17:11:59 -0500 Subject: [PATCH 4/7] Use multiple universes for TRISOs --- openmc/model/triso.py | 52 +++++++++++++------------------ tests/test_triso/inputs_true.dat | 2 +- tests/test_triso/plots.xml | 6 ++++ tests/test_triso/results_true.dat | 2 +- tests/test_triso/test_triso.py | 28 +++++++++++------ 5 files changed, 48 insertions(+), 42 deletions(-) create mode 100644 tests/test_triso/plots.xml diff --git a/openmc/model/triso.py b/openmc/model/triso.py index dc53e6ec0..ca519ed6f 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -12,18 +12,17 @@ class TRISO(object): Parameters ---------- - materials : Iterable of openmc.Material - Material to be assigned to each layer of the TRISO particle starting - with the innermost and proceeding outwards - radii : Iterable of float - Outer radii in cm of each layer of the TRISO particle in ascending order + outer_radius : int + Outer radius of TRISO particle + inner_univ : 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 ---------- - cells : list of opemc.Cell - Each layer of the TRISO particle + cell : opemc.Cell + Cell which contains the TRISO universe center : numpy.ndarray Cartesian coordinates of the center of the TRISO particle in cm outside : openmc.Region @@ -34,27 +33,18 @@ class TRISO(object): """ - def __init__(self, materials, radii, center=(0., 0., 0.)): - surfaces = [openmc.Sphere(R=r) for r in radii] - cells = [] - for i, m in enumerate(materials): - c = openmc.Cell(fill=m) - if i == 0: - c.region = -surfaces[i] - else: - c.region = +surfaces[i-1] & -surfaces[i] - cells.append(c) - self._cells = cells - self._surfaces = surfaces + def __init__(self, outer_radius, inner_univ, center=(0., 0., 0.)): + self._surface = openmc.Sphere(R=outer_radius) + self._cell = openmc.Cell(fill=inner_univ, region=-self._surface) self.center = np.asarray(center) @property def bounding_box(self): - return self.cells[-1].region.bounding_box + return self.cell.region.bounding_box @property - def cells(self): - return self._cells + def cell(self): + return self._cell @property def center(self): @@ -62,13 +52,15 @@ class TRISO(object): @property def outside(self): - return ~self.cells[-1].region.nodes[-1] + return +self._surface @center.setter def center(self, center): cv.check_type('TRISO center', center, Iterable, Real) - for s in self._surfaces: - s.x0, s.y0, s.z0 = center + self._surface.x0 = center[0] + self._surface.y0 = center[1] + self._surface.z0 = center[2] + self.cell.translation = center self._center = center def classify(self, lattice): @@ -136,11 +128,9 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): # Create copy of TRISO particle with materials preserved and # different cell/surface IDs t_copy = copy.deepcopy(t) - for c, c_copy in zip(t.cells, t_copy.cells): - c_copy.id = None - c_copy.fill = c.fill - for s in t_copy._surfaces: - s.id = None + t_copy.cell.id = None + t_copy.cell.fill = t.cell.fill + t_copy._surface.id = None triso_locations[idx].append(t_copy) # Create universes @@ -155,7 +145,7 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): u = openmc.Universe() u.add_cell(background_cell) for t in triso_list: - u.add_cells(t.cells) + u.add_cell(t.cell) iz, iy, ix = idx t.center = lattice.get_local_coordinates(t.center, (ix, iy, iz)) diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index 40f312850..3e54e5e6f 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -1 +1 @@ -6c6fecedf3db0b91b69b7b7b57b3a52c42947253bea4ee3889d6bbd6e74935cc4e599c1c743eb589a1045868412f3f88c5fef7cf10e180bdc4bd182970c9ee15 \ No newline at end of file +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 index 8f968f65a..ea7da21ed 100644 --- a/tests/test_triso/results_true.dat +++ b/tests/test_triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.685303E+00 1.121936E-01 +1.662675E+00 1.475968E-02 diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py index f7f4daa7b..d1ac4e5cd 100644 --- a/tests/test_triso/test_triso.py +++ b/tests/test_triso/test_triso.py @@ -50,22 +50,31 @@ class TRISOTestHarness(PyAPITestHarness): graphite.add_s_alpha_beta('Graph', '71t') # Create TRISO particles - materials = [fuel, porous_carbon, ipyc, sic, opyc] - radii = np.array([212.5, 312.5, 347.5, 382.5, 422.5])*1e-4 + 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 - x = random.uniform(-0.5, 0.5) - y = random.uniform(-0.5, 0.5) - z = random.uniform(-0.5, 0.5) - t = openmc.model.TRISO(materials, radii, (x, y, z)) + 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*radii[-1]: + if distance <= 2*outer_radius: break else: trisos.append(t) @@ -82,8 +91,9 @@ class TRISOTestHarness(PyAPITestHarness): # Create lattice ll, ur = box.region.bounding_box shape = (3, 3, 3) + pitch = (ur - ll) / shape lattice = openmc.model.create_triso_lattice( - trisos, ll, (ur - ll)/shape, shape, graphite) + trisos, ll, pitch, shape, graphite) box.fill = lattice root = openmc.Universe(0, cells=[box]) @@ -93,7 +103,7 @@ class TRISOTestHarness(PyAPITestHarness): settings = openmc.Settings() settings.batches = 5 settings.inactive = 0 - settings.particles = 50 + settings.particles = 100 settings.source = openmc.Source(space=openmc.stats.Point()) settings.export_to_xml() From 69f2bc10a8a628f5feb545b393150e2b5eff01a3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 Jun 2016 20:08:54 -0500 Subject: [PATCH 5/7] Warn user if TRISO particle is placed outside of lattice --- openmc/model/triso.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index ca519ed6f..6f0d145b2 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -1,6 +1,7 @@ import copy from collections import Iterable from numbers import Real +import warnings import numpy as np @@ -132,6 +133,9 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): t_copy.cell.fill = t.cell.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) From 933fdd3c818dfaefc7745af4a435308850ba3811 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 Jun 2016 22:17:46 -0500 Subject: [PATCH 6/7] Allow multiple temperatures to be used in a material --- openmc/material.py | 95 +++++++++++++++++++++------------------------- 1 file changed, 44 insertions(+), 51 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index f29a42117..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,10 +140,8 @@ 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 @@ -301,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 @@ -319,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 @@ -439,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 @@ -454,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 @@ -488,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 @@ -506,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) @@ -523,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': @@ -550,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': @@ -569,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 @@ -577,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 @@ -625,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: @@ -820,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") From 14aaaa6bb8f556b9f6775ba92f2797bcbf5f2a72 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Jun 2016 19:47:05 -0500 Subject: [PATCH 7/7] Make TRISO a special kind of openmc.Cell --- openmc/model/triso.py | 49 +++++++++++++++++-------------------------- 1 file changed, 19 insertions(+), 30 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 6f0d145b2..89e0d8aa7 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -8,60 +8,49 @@ import numpy as np import openmc import openmc.checkvalue as cv -class TRISO(object): +class TRISO(openmc.Cell): """Tristructural-isotopic (TRISO) micro fuel particle Parameters ---------- - outer_radius : int + outer_radius : float Outer radius of TRISO particle - inner_univ : openmc.Universe + 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 ---------- - cell : opemc.Cell - Cell which contains the TRISO universe + 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 - outside : openmc.Region - Region of space outside of the TRISO particle - bounding_box : tuple of numpy.ndarray - Lower-left and upper-right coordinates of an axis-aligned bounding box - for the TRISO particle + fill : openmc.Universe + Universe that contains the TRISO layers + region : openmc.Region + Region of space within the TRISO particle """ - def __init__(self, outer_radius, inner_univ, center=(0., 0., 0.)): + def __init__(self, outer_radius, fill, center=(0., 0., 0.)): self._surface = openmc.Sphere(R=outer_radius) - self._cell = openmc.Cell(fill=inner_univ, region=-self._surface) + super(TRISO, self).__init__(fill=fill, region=-self._surface) self.center = np.asarray(center) - @property - def bounding_box(self): - return self.cell.region.bounding_box - - @property - def cell(self): - return self._cell - @property def center(self): return self._center - @property - def outside(self): - return +self._surface - @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.cell.translation = center + self.translation = center self._center = center def classify(self, lattice): @@ -80,7 +69,7 @@ class TRISO(object): """ - ll, ur = self.bounding_box + 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) @@ -129,8 +118,8 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): # Create copy of TRISO particle with materials preserved and # different cell/surface IDs t_copy = copy.deepcopy(t) - t_copy.cell.id = None - t_copy.cell.fill = t.cell.fill + t_copy.id = None + t_copy.fill = t.fill t_copy._surface.id = None triso_locations[idx].append(t_copy) else: @@ -141,7 +130,7 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): 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.outside for t in triso_list]) + 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) @@ -149,7 +138,7 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): u = openmc.Universe() u.add_cell(background_cell) for t in triso_list: - u.add_cell(t.cell) + u.add_cell(t) iz, iy, ix = idx t.center = lattice.get_local_coordinates(t.center, (ix, iy, iz))