crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge remote-tracking branch 'upstream/develop' into mg_docs

Browse source
This commit is contained in:
Adam Nelson 2016-06-15 19:36:49 -04:00
commit 5651505918
28 changed files with 1218 additions and 1088 deletions
Download patch file Download diff file Expand all files Collapse all files

180
docs/source/pythonapi/examples/mgxs-part-iv.ipynb
View file

File diff suppressed because one or more lines are too long

27
docs/source/pythonapi/index.rst
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

8
docs/source/usersguide/mgxs_library.rst
View file

@ -173,9 +173,9 @@ attributes/sub-elements required to describe the meta-data:
provided via the ``scatt_type`` element above, is represented and thus used
during the scattering process. Specifically, the options are to either
convert the Legendre expansion to a tabular representation or leave it as
a set of Legendre coefficients. Converting to a tabular representation will
cost memory but can allow for a decrease in runtime compared to leaving as a
set of Legendre coefficients. This element has the following
a set of Legendre coefficients. Converting to a tabular representation
will cost memory but can allow for a decrease in runtime compared to
leaving as a set of Legendre coefficients. This element has the following
attributes/sub-elements:
:enable:
@ -183,7 +183,7 @@ attributes/sub-elements required to describe the meta-data:
tabular format should be performed or not. A value of "true" means
the conversion should be performed, "false" means it should not.
*Default*: "false"
*Default*: "true"
:num_points:
If the conversion is to take place the number of tabular points is

112
openmc/arithmetic.py
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

2
openmc/element.py
View file

@ -126,7 +126,7 @@ class Element(object):
"""
isotopes = []
for isotope, abundance in natural_abundance.items():
for isotope, abundance in sorted(natural_abundance.items()):
if isotope.startswith(self.name + '-'):
nuc = openmc.Nuclide(isotope, self.xs)
isotopes.append((nuc, abundance))

21
openmc/filter.py
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)

120
openmc/material.py
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")

23
openmc/mesh.py
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

65
openmc/mgxs/library.py
View file

@ -753,7 +753,8 @@ class Library(object):
return pickle.load(open(full_filename, 'rb'))
def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
xs_id='1m', order=None):
xs_id='1m', order=None, tabular_legendre=None,
tabular_points=33):
"""Generates an openmc.XSdata object describing a multi-group cross section
data set for eventual combination in to an openmc.MGXSLibrary object
(i.e., the library).
@ -773,9 +774,22 @@ class Library(object):
nuclide this will be set to 'macro' regardless.
xs_ids : str
Cross section set identifier. Defaults to '1m'.
order : Scattering order for this data entry. Default is None,
order : int
Scattering order for this data entry. Default is None,
which will set the XSdata object to use the order of the
Library.
tabular_legendre : None or bool
Flag to denote whether or not the Legendre expansion of the
scattering angular distribution is to be converted to a tabular
representation by OpenMC. A value of `True` means that it is to be
converted while a value of `False` means that it will not be.
Defaults to `None` which leaves the default behavior of OpenMC in
place (the distribution is converted to a tabular representation).
tabular_points : int
This parameter is not used unless the ``tabular_legendre``
parameter is set to `True`. In this case, this parameter sets the
number of equally-spaced points in the domain of [-1,1] to be used
in building the tabular distribution. Default is `33`.
Returns
-------
@ -804,6 +818,10 @@ class Library(object):
if order is not None:
cv.check_greater_than('order', order, 0, equality=True)
cv.check_less_than('order', order, 10, equality=True)
cv.check_type('tabular_legendre', tabular_legendre,
(type(None), bool))
if tabular_points is not None:
cv.check_greater_than('tabular_points', tabular_points, 1)
# Make sure statepoint has been loaded
if self._sp_filename is None:
@ -830,6 +848,11 @@ class Library(object):
# the provided order or the Library's order.
xsdata.order = min(order, self.legendre_order)
# Set the tabular_legendre option if needed
if tabular_legendre is not None:
xsdata.tabular_legendre = {'enable': tabular_legendre,
'num_points': tabular_points}
if nuclide is not 'total':
xsdata.zaid = self._nuclides[nuclide][0]
xsdata.awr = self._nuclides[nuclide][1]
@ -906,7 +929,8 @@ class Library(object):
return xsdata
def create_mg_library(self, xs_type='macro', xsdata_names=None,
xs_ids=None):
xs_ids=None, tabular_legendre=None,
tabular_points=33):
"""Creates an openmc.MGXSLibrary object to contain the MGXS data for the
Multi-Group mode of OpenMC.
@ -923,6 +947,18 @@ class Library(object):
Cross section set identifier (i.e., '71c') for all
data sets (if only str) or for each individual one
(if iterable of str). Defaults to '1m'.
tabular_legendre : None or bool
Flag to denote whether or not the Legendre expansion of the
scattering angular distribution is to be converted to a tabular
representation by OpenMC. A value of `True` means that it is to be
converted while a value of `False` means that it will not be.
Defaults to `None` which leaves the default behavior of OpenMC in
place (the distribution is converted to a tabular representation).
tabular_points : int
This parameter is not used unless the ``tabular_legendre``
parameter is set to `True`. In this case, this parameter sets the
number of equally-spaced points in the domain of [-1,1] to be used
in building the tabular distribution. Default is `33`.
Returns
-------
@ -983,13 +1019,16 @@ class Library(object):
xsdata_name += '_' + nuclide
xsdata = self.get_xsdata(domain, xsdata_name, nuclide=nuclide,
xs_type=xs_type, xs_id=xs_ids[i])
xs_type=xs_type, xs_id=xs_ids[i],
tabular_legendre=tabular_legendre,
tabular_points=tabular_points)
mgxs_file.add_xsdata(xsdata)
return mgxs_file
def create_mg_mode(self, xsdata_names=None, xs_ids=None):
def create_mg_mode(self, xsdata_names=None, xs_ids=None,
tabular_legendre=None, tabular_points=33):
"""Creates an openmc.MGXSLibrary object to contain the MGXS data for the
Multi-Group mode of OpenMC as well as the associated openmc.Materials
and openmc.Geometry objects. The created Geometry is the same as that
@ -1007,6 +1046,18 @@ class Library(object):
Cross section set identifier (i.e., '71c') for all
data sets (if only str) or for each individual one
(if iterable of str). Defaults to '1m'.
tabular_legendre : None or bool
Flag to denote whether or not the Legendre expansion of the
scattering angular distribution is to be converted to a tabular
representation by OpenMC. A value of `True` means that it is to be
converted while a value of `False` means that it will not be.
Defaults to `None` which leaves the default behavior of OpenMC in
place (the distribution is converted to a tabular representation).
tabular_points : int
This parameter is not used unless the ``tabular_legendre``
parameter is set to `True`. In this case, this parameter sets the
number of equally-spaced points in the domain of [-1,1] to be used
in building the tabular distribution. Default is `33`.
Returns
-------
@ -1071,7 +1122,9 @@ class Library(object):
# Create XSdata and Macroscopic for this domain
xsdata = self.get_xsdata(domain, xsdata_name, nuclide='total',
xs_type=xs_type, xs_id=xs_ids[i])
xs_type=xs_type, xs_id=xs_ids[i],
tabular_legendre=tabular_legendre,
tabular_points=tabular_points)
mgxs_file.add_xsdata(xsdata)
macroscopic = openmc.Macroscopic(name=xsdata_name, xs=xs_ids[i])

2
openmc/mgxs_library.py
View file

@ -445,7 +445,7 @@ class XSdata(object):
raise ValueError(msg)
if 'num_points' in tabular_legendre:
num_points = tabular_legendre['num_points']
check_value('num_points', num_points, Integral)
check_type('num_points', num_points, Integral)
check_greater_than('num_points', num_points, 0)
else:
if not enable:

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

47
openmc/tallies.py
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

19
openmc/trigger.py
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

47
openmc/universe.py
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.

3
setup.py
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

10
src/mgxs_header.F90
View file

@ -387,11 +387,13 @@ module mgxs_header
! Get scattering treatment information
! Tabular_legendre tells us if we are to treat the provided
! Legendre polynomials as tabular data (if enable is true) or leaving
! them as Legendres (if enable is false, or the default)
! Legendre polynomials as tabular data (if enable is true, default)
! or leaving them as Legendres (if enable is false)
! Set the default (leave as Legendre polynomials)
enable_leg_mu = .false.
! Set the default (Convert to tabular format w/ 33 points)
enable_leg_mu = .true.
legendre_mu_points = 33
! Get the user-provided values
if (check_for_node(node_xsdata, "tabular_legendre")) then
call get_node_ptr(node_xsdata, "tabular_legendre", node_legendre_mu)
if (check_for_node(node_legendre_mu, "enable")) then

22
src/scattdata_header.F90
View file

@ -398,29 +398,18 @@ contains
end do
! Re-normalize fmu for numerical integration issues and in case
! the negative fix-up introduced un-normalized data
! the negative fix-up introduced un-normalized data while
! accruing the CDF
norm = ZERO
do imu = 2, order
norm = norm + HALF * this % dmu * &
(this % fmu(gin) % data(imu - 1, gout) + &
this % fmu(gin) % data(imu, gout))
this % dist(gin) % data(imu, gout) = norm
end do
if (norm > ZERO) then
this % fmu(gin) % data(:, gout) = &
this % fmu(gin) % data(:, gout) / norm
end if
! Now create CDF from fmu with trapezoidal rule
this % dist(gin) % data(1, gout) = ZERO
do imu = 2, order
this % dist(gin) % data(imu, gout) = &
this % dist(gin) % data(imu - 1, gout) + &
HALF * this % dmu * (this % fmu(gin) % data(imu - 1, gout) + &
this % fmu(gin) % data(imu, gout))
end do
! Ensure we normalize to 1 still
norm = this % dist(gin) % data(order, gout)
if (norm > ZERO) then
this % dist(gin) % data(:, gout) = &
this % dist(gin) % data(:, gout) / norm
end if
@ -628,11 +617,10 @@ contains
p1 = this % fmu(gin) % data(k + 1, gout)
mu1 = this % mu(k + 1)
frac = (p1 - p0) / (mu1 - mu0)
if (frac == ZERO) then
if (p0 == p1) then
mu = mu0 + (xi - c_k) / p0
else
frac = (p1 - p0) / (mu1 - mu0)
mu = mu0 + &
(sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac
end if

2
src/tally.F90
View file

@ -107,8 +107,6 @@ contains
real(8) :: atom_density_ ! atom/b-cm
real(8) :: f ! interpolation factor
real(8) :: score ! analog tally score
real(8) :: macro_total ! material macro total xs
real(8) :: macro_scatt ! material macro scatt xs
real(8) :: E ! particle energy
i = 0

2
tests/test_mg_basic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.033731E+00 4.974463E-02
1.047136E+00 2.765964E-02

2
tests/test_mg_max_order/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.055274E+00 1.715904E-02
1.102093E+00 4.962190E-02

2
tests/test_mg_nuclide/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.033731E+00 4.974463E-02
1.047136E+00 2.765964E-02

1306
tests/test_mg_tallies/results_true.dat
View file

File diff suppressed because it is too large Load diff

2
tests/test_mgxs_library_ce_to_mg/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.094839E+00 1.203524E-02
1.140804E+00 2.937150E-02

1
tests/test_triso/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
2dcfd1a17cba671874e60192a7355deb57e2e51467a474fd168c8b51e454a977edb34df07ae11625c0a43906112152c75113e442a9a8f240a4c9d1a11ee4771d

6
tests/test_triso/plots.xml Normal file
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>

2
tests/test_triso/results_true.dat Normal file
View file

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

117
tests/test_triso/test_triso.py Normal file
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()