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

Fixed merge conflicts with develop

Browse source
This commit is contained in:
Will Boyd 2015-07-26 16:09:15 -07:00
commit 4e6aca6028
331 changed files with 6506 additions and 222815 deletions
Download patch file Download diff file Expand all files Collapse all files

138
openmc/checkvalue.py
View file

@ -1,13 +1,137 @@
import numpy as np
def check_type(name, value, expected_type, expected_iter_type=None):
"""Ensure that an object is of an expected type. Optionally, if the object is
iterable, check that each element is of a particular type.
Parameters
----------
name : str
Description of value being checked
value : object
Object to check type of
expected_type : type
type to check object against
expected_iter_type : type or None, optional
Expected type of each element in value, assuming it is iterable. If
None, no check will be performed.
"""
if not isinstance(value, expected_type):
msg = 'Unable to set {0} to {1} which is not of type {2}'.format(
name, value, expected_type.__name__)
raise ValueError(msg)
if expected_iter_type:
for item in value:
if not isinstance(item, expected_iter_type):
msg = 'Unable to set {0} to {1} since each item must be ' \
'of type {2}'.format(name, value,
expected_iter_type.__name__)
raise ValueError(msg)
def is_integer(val):
return isinstance(val, (int, np.int32, np.int64))
def check_length(name, value, length_min, length_max=None):
"""Ensure that a sized object has length within a given range.
Parameters
----------
name : str
Description of value being checked
value : collections.Sized
Object to check length of
length_min : int
Minimum length of object
length_max : int or None, optional
Maximum length of object. If None, it is assumed object must be of
length length_min.
"""
if length_max is None:
if len(value) != length_min:
msg = 'Unable to set {0} to {1} since it must be of ' \
'length {2}'.format(name, value, length_min)
raise ValueError(msg)
elif not length_min <= len(value) <= length_max:
if length_min == length_max:
msg = 'Unable to set {0} to {1} since it must be of ' \
'length {2}'.format(name, value, length_min)
else:
msg = 'Unable to set {0} to {1} since it must have length ' \
'between {2} and {3}'.format(name, value, length_min,
length_max)
raise ValueError(msg)
def is_float(val):
return isinstance(val, (float, np.float32, np.float64))
def check_value(name, value, accepted_values):
"""Ensure that an object's value is contained in a set of acceptable values.
Parameters
----------
name : str
Description of value being checked
value : collections.Iterable
Object to check
accepted_values : collections.Container
Container of acceptable values
def is_string(val):
return isinstance(val, (str, np.str))
"""
if value not in accepted_values:
msg = 'Unable to set {0} to {1} since it is not in {2}'.format(
name, value, accepted_values)
raise ValueError(msg)
def check_less_than(name, value, maximum, equality=False):
"""Ensure that an object's value is less than a given value.
Parameters
----------
name : str
Description of the value being checked
value : object
Object to check
maximum : object
Maximum value to check against
equality : bool, optional
Whether equality is allowed. Defaluts to False.
"""
if equality:
if value > maximum:
msg = 'Unable to set {0} to {1} since it is greater than ' \
'{2}'.format(name, value, maximum)
raise ValueError(msg)
else:
if value >= maximum:
msg = 'Unable to set {0} to {1} since it is greater than ' \
'or equal to {2}'.format(name, value, maximum)
raise ValueError(msg)
def check_greater_than(name, value, minimum, equality=False):
"""Ensure that an object's value is less than a given value.
Parameters
----------
name : str
Description of the value being checked
value : object
Object to check
minimum : object
Minimum value to check against
equality : bool, optional
Whether equality is allowed. Defaluts to False.
"""
if equality:
if value < minimum:
msg = 'Unable to set {0} to {1} since it is less than ' \
'{2}'.format(name, value, minimum)
raise ValueError(msg)
else:
if value <= minimum:
msg = 'Unable to set {0} to {1} since it is less than ' \
'or equal to {2}'.format(name, value, minimum)
raise ValueError(msg)

746
openmc/cmfd.py
View file

@ -1,15 +1,77 @@
"""This module can be used to specify parameters used for coarse mesh finite
difference (CMFD) acceleration in OpenMC. CMFD was first proposed by [Smith]_
and is widely used in accelerating neutron transport problems.
References
----------
.. [Smith] K. Smith, "Nodal method storage reduction by non-linear
iteration", *Trans. Am. Nucl. Soc.*, **44**, 265 (1983).
"""
from collections import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than, check_less_than)
if sys.version_info[0] >= 3:
basestring = str
class CMFDMesh(object):
"""A structured Cartesian mesh used for Coarse Mesh Finite Difference (CMFD)
acceleration.
Attributes
----------
lower_left : Iterable of float
The lower-left corner of the structured mesh. If only two coordinates are
given, it is assumed that the mesh is an x-y mesh.
upper_right : Iterable of float
The upper-right corner of the structrued mesh. If only two coordinates
are given, it is assumed that the mesh is an x-y mesh.
dimension : Iterable of int
The number of mesh cells in each direction.
width : Iterable of float
The width of mesh cells in each direction.
energy : Iterable of float
Energy bins in MeV, listed in ascending order (e.g. [0.0, 0.625e-7,
20.0]) for CMFD tallies and acceleration. If no energy bins are listed,
OpenMC automatically assumes a one energy group calculation over the
entire energy range.
albedo : Iterable of float
Surface ratio of incoming to outgoing partial currents on global
boundary conditions. They are listed in the following order: -x +x -y +y
-z +z.
map : Iterable of int
An optional acceleration map can be specified to overlay on the coarse
mesh spatial grid. If this option is used, a ``1`` is used for a
non-accelerated region and a ``2`` is used for an accelerated region.
For a simple 4x4 coarse mesh with a 2x2 fuel lattice surrounded by
reflector, the map is:
::
[1, 1, 1, 1,
1, 2, 2, 1,
1, 2, 2, 1,
1, 1, 1, 1]
Therefore a 2x2 system of equations is solved rather than a 4x4. This is
extremely important to use in reflectors as neutrons will not contribute
to any tallies far away from fission source neutron regions. A ``2``
must be used to identify any fission source region.
"""
def __init__(self):
self._lower_left = None
self._upper_right = None
self._dimension = None
@ -18,625 +80,441 @@ class CMFDMesh(object):
self._albedo = None
self._map = None
@property
def lower_left(self):
return self._lower_left
@property
def upper_right(self):
return self._upper_right
@property
def dimension(self):
return self._dimension
@property
def width(self):
return self._width
@property
def energy(self):
return self._energy
@property
def albedo(self):
return self._albedo
@property
def map(self):
return self._mape
return self._map
@lower_left.setter
def lower_left(self, lower_left):
if not isinstance(lower_left, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with lower_left {0} which is ' \
'not a Python list, tuple or NumPy array'.format(lower_left)
raise ValueError(msg)
elif len(lower_left) != 2 and len(lower_left) != 3:
msg = 'Unable to set CMFD Mesh with lower_left {0} since it ' \
'must include 2 or 3 dimensions'.format(lower_left)
raise ValueError(msg)
for coord in lower_left:
if not is_integer(coord) and not is_float(coord):
msg = 'Unable to set CMFD Mesh with lower_left {0} which is ' \
'not an integer or a floating point value'.format(coord)
raise ValueError(msg)
check_type('CMFD mesh lower_left', lower_left, Iterable, Real)
check_length('CMFD mesh lower_left', lower_left, 2, 3)
self._lower_left = lower_left
@upper_right.setter
def upper_right(self, upper_right):
if not isinstance(upper_right, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with upper_right {0} which is ' \
'not a Python list, tuple or NumPy array'.format(upper_right)
raise ValueError(msg)
if len(upper_right) != 2 and len(upper_right) != 3:
msg = 'Unable to set CMFD Mesh with upper_right {0} since it ' \
'must include 2 or 3 dimensions'.format(upper_right)
raise ValueError(msg)
for coord in upper_right:
if not is_integer(coord) and not is_float(coord):
msg = 'Unable to set CMFD Mesh with upper_right {0} which ' \
'is not an integer or floating point value'.format(coord)
raise ValueError(msg)
check_type('CMFD mesh upper_right', upper_right, Iterable, Real)
check_length('CMFD mesh upper_right', upper_right, 2, 3)
self._upper_right = upper_right
@dimension.setter
def dimension(self, dimension):
if not isinstance(dimension, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with dimension {0} which is ' \
'not a Python list, tuple or NumPy array'.format(dimension)
raise ValueError(msg)
elif len(dimension) != 2 and len(dimension) != 3:
msg = 'Unable to set CMFD Mesh with dimension {0} since it ' \
'must include 2 or 3 dimensions'.format(dimension)
raise ValueError(msg)
for dim in dimension:
if not is_integer(dim):
msg = 'Unable to set CMFD Mesh with dimension {0} which ' \
'is a non-integer'.format(dim)
raise ValueError(msg)
check_type('CMFD mesh dimension', dimension, Iterable, Integral)
check_length('CMFD mesh dimension', dimension, 2, 3)
self._dimension = dimension
@width.setter
def width(self, width):
if not width is None:
if not isinstance(width, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with width {0} which ' \
'is not a Python list, tuple or NumPy array'.format(width)
raise ValueError(msg)
if len(width) != 2 and len(width) != 3:
msg = 'Unable to set CMFD Mesh with width {0} since it must ' \
'include 2 or 3 dimensions'.format(width)
raise ValueError(msg)
for dim in width:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to set CMFD Mesh with width {0} which is ' \
'not an integer or floating point value'.format(width)
raise ValueError(msg)
check_type('CMFD mesh width', width, Iterable, Real)
check_length('CMFD mesh width', width, 2, 3)
self._width = width
@energy.setter
def energy(self, energy):
if not isinstance(energy, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh energy to {0} which is not ' \
'a Python tuple/list or NumPy array'.format(energy)
raise ValueError(msg)
check_type('CMFD mesh energy', energy, Iterable, Real)
for e in energy:
if not is_integer(e) and not is_float(e):
msg = 'Unable to set CMFD Mesh energy to {0} which is not ' \
'an integer or floating point value'.format(e)
raise ValueError(msg)
elif e < 0:
msg = 'Unable to set CMFD Mesh energy to {0} which is ' \
'is a negative integer'.format(e)
raise ValueError(msg)
check_greater_than('CMFD mesh energy', e, 0, True)
self._energy = energy
@albedo.setter
def albedo(self, albedo):
if not isinstance(albedo, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \
'a Python tuple/list or NumPy array'.format(albedo)
raise ValueError(msg)
if not len(albedo) == 6:
msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \
'length 6 for +/-x,y,z'.format(albedo)
raise ValueError(msg)
check_type('CMFD mesh albedo', albedo, Iterable, Real)
check_length('CMFD mesh albedo', albedo, 6)
for a in albedo:
if not is_integer(a) and not is_float(a):
msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \
'an integer or floating point value'.format(a)
raise ValueError(msg)
elif a < 0 or a > 1:
msg = 'Unable to set CMFD Mesh albedo to {0} which is ' \
'is not in [0,1]'.format(a)
raise ValueError(msg)
check_greater_than('CMFD mesh albedo', a, 0, True)
check_less_than('CMFD mesh albedo', a, 1, True)
self._albedo = albedo
@map.setter
def map(self, map):
if not isinstance(map, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh map to {0} which is not ' \
'a Python tuple/list or NumPy array'.format(map)
raise ValueError(msg)
for m in map:
if m != 1 and m != 2:
msg = 'Unable to set CMFD Mesh map to {0} which is ' \
'is not 1 or 2'.format(m)
raise ValueError(msg)
self._map = map
def get_mesh_xml(self):
def map(self, meshmap):
check_type('CMFD mesh map', meshmap, Iterable, Integral)
for m in meshmap:
check_value('CMFD mesh map', m, [1, 2])
self._map = meshmap
def _get_xml_element(self):
element = ET.Element("mesh")
if len(self._lower_left) == 2:
subelement = ET.SubElement(element, "lower_left")
subelement.text = '{0} {1}'.format(self._lower_left[0],
self._lower_left[1])
else:
subelement = ET.SubElement(element, "lower_left")
subelement.text = '{0} {1} {2}'.format(self._lower_left[0],
self._lower_left[1],
self._lower_left[2])
subelement = ET.SubElement(element, "lower_left")
subelement.text = ' '.join(map(str, self._lower_left))
if not self._upper_right is None:
if len(self._upper_right) == 2:
subelement = ET.SubElement(element, "upper_right")
subelement.text = '{0} {1}'.format(self._upper_right[0],
self._upper_right[1])
else:
subelement = ET.SubElement(element, "upper_right")
subelement.text = '{0} {1} {2}'.format(self._upper_right[0],
self._upper_right[1],
self._upper_right[2])
if self.upper_right is not None:
subelement = ET.SubElement(element, "upper_right")
subelement.text = ' '.join(map(str, self.upper_right))
if len(self._dimension) == 2:
subelement = ET.SubElement(element, "dimension")
subelement.text = '{0} {1}'.format(self._dimension[0],
self._dimension[1])
else:
subelement = ET.SubElement(element, "dimension")
subelement.text = '{0} {1} {2}'.format(self._dimension[0],
self._dimension[1],
self._dimension[2])
subelement = ET.SubElement(element, "dimension")
subelement.text = ' '.join(map(str, self.dimension))
if not self._width is None:
if len(self._width) == 2:
subelement = ET.SubElement(element, "width")
subelement.text = '{0} {1}'.format(self._width[0],
self._width[1])
else:
subelement = ET.SubElement(element, "width")
subelement.text = '{0} {1} {2}'.format(self._width[0],
self._width[1],
self._width[2])
if not self._energy is None:
if self.width is not None:
subelement = ET.SubElement(element, "width")
subelement.text = ' '.join(map(str, self.width))
if self.energy is not None:
subelement = ET.SubElement(element, "energy")
subelement.text = ' '.join(map(str, self.energy))
energy = ''
for e in self._energy:
energy += '{0} '.format(e)
subelement.set("energy", energy.rstrip(' '))
if not self._albedo is None:
if self.albedo is not None:
subelement = ET.SubElement(element, "albedo")
subelement.text = ' '.join(map(str, self.albedo))
albedo = ''
for a in self._albedo:
albedo += '{0} '.format(a)
subelement.set("albedo", albedo.rstrip(' '))
if not self._map is None:
if self.map is not None:
subelement = ET.SubElement(element, "map")
map = ''
for m in self._map:
map += '{0} '.format(m)
subelement.set("map", map.rstrip(' '))
subelement.text = ' '.join(map(str, self.map))
return element
class CMFDFile(object):
"""Parameters that control the use of coarse-mesh finite difference acceleration
in OpenMC. This corresponds directly to the cmfd.xml input file.
Attributes
----------
begin : int
Batch number at which CMFD calculations should begin
dhat_reset : bool
Indicate whether :math:`\widehat{D}` nonlinear CMFD parameters should be
reset to zero before solving CMFD eigenproblem.
display : {'balance', 'dominance', 'entropy', 'source'}
Set one additional CMFD output column. Options are:
* "balance" - prints the RMS [%] of the resdiual from the neutron balance
equation on CMFD tallies.
* "dominance" - prints the estimated dominance ratio from the CMFD
iterations.
* "entropy" - prints the *entropy* of the CMFD predicted fission source.
* "source" - prints the RMS [%] between the OpenMC fission source and
CMFD fission source.
downscatter : bool
Indicate whether an effective downscatter cross section should be used
when using 2-group CMFD.
feedback : bool
Indicate or not the CMFD diffusion result is used to adjust the weight
of fission source neutrons on the next OpenMC batch. Defaults to False.
gauss_seidel_tolerance : Iterable of float
Two parameters specifying the absolute inner tolerance and the relative
inner tolerance for Gauss-Seidel iterations when performing CMFD.
ktol : float
Tolerance on the eigenvalue when performing CMFD power iteration
cmfd_mesh : CMFDMesh
Structured mesh to be used for acceleration
norm : float
Normalization factor applied to the CMFD fission source distribution
power_monitor : bool
View convergence of power iteration during CMFD acceleration
run_adjoint : bool
Perform adjoint calculation on the last batch
shift : float
Optional Wielandt shift parameter for accelerating power iterations. By
default, it is very large so there is effectively no impact.
spectral : float
Optional spectral radius that can be used to accelerate the convergence
of Gauss-Seidel iterations during CMFD power iteration.
stol : float
Tolerance on the fission source when performing CMFD power iteration
tally_reset : list of int
List of batch numbers at which CMFD tallies should be reset
write_matrices : bool
Write sparse matrices that are used during CMFD acceleration (loss,
production) to file
"""
def __init__(self):
self._active_flush = None
self._begin = None
self._dhat_reset = None
self._display = None
self._downscatter = None
self._feedback = None
self._inactive = None
self._inactive_flush = None
self._gauss_seidel_tolerance = None
self._ktol = None
self._cmfd_mesh = None
self._norm = None
self._num_flushes = None
self._power_monitor = None
self._run_adjoint = None
self._shift = None
self._spectral = None
self._stol = None
self._tally_reset = None
self._write_matrices = None
self._cmfd_file = ET.Element("cmfd")
self._cmfd_mesh_element = None
@property
def active_flush(self):
return self._active_flush
@property
def begin(self):
return self._begin
@property
def dhat_reset(self):
return self._dhat_reset
@property
def display(self):
return self._display
@property
def downscatter(self):
return self._downscatter
@property
def feedback(self):
return self._feedback
@property
def gauss_seidel_tolerance(self):
return self._gauss_seidel_tolerance
@property
def inactive(self):
return self._inactive
@property
def inactive_flush(self):
return self._inactive_flush
def ktol(self):
return self._ktol
@property
def cmfd_mesh(self):
return self._cmfd_mesh
@property
def norm(self):
return self._norm
@property
def num_flushes(self):
return self._num_flushes
@property
def power_monitor(self):
return self._power_monitor
@property
def run_adjoint(self):
return self._run_adjoint
@property
def shift(self):
return self._shift
@property
def solver(self):
return self._solver
def spectral(self):
return self._spectral
@property
def stol(self):
return self._stol
@property
def tally_reset(self):
return self._tally_reset
@property
def write_matrices(self):
return self._write_matrices
@active_flush.setter
def active_flush(self, active_flush):
if not is_integer(active_flush):
msg = 'Unable to set CMFD active flush batch to a non-integer ' \
'value {0}'.format(active_flush)
raise ValueError(msg)
if active_flush < 0:
msg = 'Unable to set CMFD active flush batch to a negative ' \
'value {0}'.format(active_flush)
raise ValueError(msg)
self._active_flush = active_flush
@begin.setter
def begin(self, begin):
if not is_integer(begin):
msg = 'Unable to set CMFD begin batch to a non-integer ' \
'value {0}'.format(begin)
raise ValueError(msg)
if begin <= 0:
msg = 'Unable to set CMFD begin batch batch to a negative ' \
'value {0}'.format(begin)
raise ValueError(msg)
check_type('CMFD begin batch', begin, Integral)
check_greater_than('CMFD begin batch', begin, 0)
self._begin = begin
@dhat_reset.setter
def dhat_reset(self, dhat_reset):
check_type('CMFD Dhat reset', dhat_reset, bool)
self._dhat_reset = dhat_reset
@display.setter
def display(self, display):
if not is_string(display):
msg = 'Unable to set CMFD display to a non-string ' \
'value'.format(display)
raise ValueError(msg)
if display not in ['balance', 'dominance', 'entropy', 'source']:
msg = 'Unable to set CMFD display to {0} which is ' \
'not an accepted value'.format(display)
raise ValueError(msg)
check_type('CMFD display', display, basestring)
check_value('CMFD display', display,
['balance', 'dominance', 'entropy', 'source'])
self._display = display
@downscatter.setter
def downscatter(self, downscatter):
check_type('CMFD downscatter', downscatter, bool)
self._downscatter = downscatter
@feedback.setter
def feedback(self, feedback):
if not isinstance(feedback, bool):
msg = 'Unable to set CMFD feedback to {0} which is ' \
'a non-boolean value'.format(feedback)
raise ValueError(msg)
check_type('CMFD feedback', feedback, bool)
self._feedback = feedback
@gauss_seidel_tolerance.setter
def gauss_seidel_tolerance(self, gauss_seidel_tolerance):
check_type('CMFD Gauss-Seidel tolerance', gauss_seidel_tolerance,
Iterable, Real)
check_length('Gauss-Seidel tolerance', gauss_seidel_tolerance, 2)
self._gauss_seidel_tolerance = gauss_seidel_tolerance
@inactive.setter
def inactive(self, inactive):
if not isinstance(inactive, bool):
msg = 'Unable to set CMFD inactive batch to {0} which is ' \
' a non-boolean value'.format(inactive)
raise ValueError(msg)
self._inactive = inactive
@inactive_flush.setter
def inactive_flush(self, inactive_flush):
if not is_integer(inactive_flush):
msg = 'Unable to set CMFD inactive flush batch to {0} which is ' \
'a non-integer value'.format(inactive_flush)
raise ValueError(msg)
if inactive_flush <= 0:
msg = 'Unable to set CMFD inactive flush batch to {0} which is ' \
'a negative value {0}'.format(inactive_flush)
raise ValueError(msg)
self._inactive_flush = inactive_flush
@ktol.setter
def ktol(self, ktol):
check_type('CMFD eigenvalue tolerance', ktol, Real)
self._ktol = ktol
@cmfd_mesh.setter
def cmfd_mesh(self, mesh):
if not isinstance(mesh, CMFDMesh):
msg = 'Unable to set CMFD mesh to {0} which is not a ' \
'CMFDMesh object'.format(mesh)
raise ValueError(msg)
check_type('CMFD mesh', mesh, CMFDMesh)
self._mesh = mesh
@norm.setter
def norm(self, norm):
if not is_integer(norm) and not is_float(norm):
msg = 'Unable to set the CMFD norm to {0} which is not ' \
'an integer or floating point value'.format(norm)
raise ValueError(msg)
check_type('CMFD norm', norm, Real)
self._norm = norm
@num_flushes.setter
def num_flushes(self, num_flushes):
if not is_integer(num_flushes):
msg = 'Unable to set the CMFD number of flushes to {0} ' \
'which is not an integer value'.format(num_flushes)
raise ValueError(msg)
if num_flushes < 0:
msg = 'Unable to set CMFD number of flushes to a negative ' \
'value {0}'.format(num_flushes)
raise ValueError(msg)
self._num_flushes = num_flushes
@power_monitor.setter
def power_monitor(self, power_monitor):
if not isinstance(power_monitor, bool):
msg = 'Unable to set CMFD power monitor to {0} which is a ' \
'non-boolean value'.format(power_monitor)
raise ValueError(msg)
check_type('CMFD power monitor', power_monitor, bool)
self._power_monitor = power_monitor
@run_adjoint.setter
def run_adjoint(self, run_adjoint):
if not isinstance(run_adjoint, bool):
msg = 'Unable to set CMFD run adjoint to {0} which is a ' \
'non-boolean value'.format(run_adjoint)
raise ValueError(msg)
check_type('CMFD run adjoint', run_adjoint, bool)
self._run_adjoint = run_adjoint
@shift.setter
def shift(self, shift):
check_type('CMFD Wielandt shift', shift, Real)
self._shift = shift
@spectral.setter
def spectral(self, spectral):
check_type('CMFD spectral radius', spectral, Real)
self._spectral = spectral
@stol.setter
def stol(self, stol):
check_type('CMFD fission source tolerance', stol, Real)
self._stol = stol
@tally_reset.setter
def tally_reset(self, tally_reset):
check_type('tally reset batches', tally_reset, Iterable, Integral)
self._tally_reset = tally_reset
@write_matrices.setter
def write_matrices(self, write_matrices):
if not isinstance(write_matrices, bool):
msg = 'Unable to set CMFD write matrices to {0} which is a ' \
'non-boolean value'.format(write_matrices)
raise ValueError(msg)
check_type('CMFD write matrices', write_matrices, bool)
self._write_matrices = write_matrices
def create_active_flush_subelement(self):
if not self._active_flush is None:
element = ET.SubElement(self._cmfd_file, "active_flush")
element.text = '{0}'.format(str(self._active_flush))
def create_begin_subelement(self):
if not self._begin is None:
def _create_begin_subelement(self):
if self._begin is not None:
element = ET.SubElement(self._cmfd_file, "begin")
element.text = '{0}'.format(str(self._begin))
element.text = str(self._begin)
def _create_dhat_reset_subelement(self):
if self._dhat_reset is not None:
element = ET.SubElement(self._cmfd_file, "dhat_reset")
element.text = str(self._dhat_reset).lower()
def create_display_subelement(self):
if not self._display is None:
def _create_display_subelement(self):
if self._display is not None:
element = ET.SubElement(self._cmfd_file, "display")
element.text = '{0}'.format(str(self._display))
element.text = str(self._display)
def _create_downscatter_subelement(self):
if self._downscatter is not None:
element = ET.SubElement(self._cmfd_file, "downscatter")
element.text = str(self._downscatter).lower()
def create_feedback_subelement(self):
if not self._feedback is None:
def _create_feedback_subelement(self):
if self._feedback is not None:
element = ET.SubElement(self._cmfd_file, "feeback")
element.text = '{0}'.format(str(self._feedback).lower())
element.text = str(self._feedback).lower()
def _create_gauss_seidel_tolerance_subelement(self):
if self._gauss_seidel_tolerance is not None:
element = ET.SubElement(self._cmfd_file, "gauss_seidel_tolerance")
element.text = ' '.join(map(str, self._gauss_seidel_tolerance))
def create_inactive_subelement(self):
def _create_ktol_subelement(self):
if self._ktol is not None:
element = ET.SubElement(self._ktol, "ktol")
element.text = str(self._ktol)
if not self._inactive is None:
element = ET.SubElement(self._cmfd_file, "inactive")
element.text = '{0}'.format(str(self._inactive).lower())
def create_inactive_flush_subelement(self):
if not self._inactive_flush is None:
element = ET.SubElement(self._cmfd_file, "inactive_flush")
element.text = '{0}'.format(str(self._inactive_flush))
def create_mesh_subelement(self):
if not self._mesh is None:
xml_element = self._mesh.get_mesh_xml()
def _create_mesh_subelement(self):
if self._mesh is not None:
xml_element = self._mesh._get_xml_element()
self._cmfd_file.append(xml_element)
def create_norm_subelement(self):
if not self._num_flushes is None:
def _create_norm_subelement(self):
if self._norm is not None:
element = ET.SubElement(self._cmfd_file, "norm")
element.text = '{0}'.format(str(self._norm))
element.text = str(self._norm)
def create_num_flushes_subelement(self):
if not self._num_flushes is None:
element = ET.SubElement(self._cmfd_file, "num_flushes")
element.text = '{0}'.format(str(self._num_flushes))
def create_power_monitor_subelement(self):
if not self._power_monitor is None:
def _create_power_monitor_subelement(self):
if self._power_monitor is not None:
element = ET.SubElement(self._cmfd_file, "power_monitor")
element.text = '{0}'.format(str(self._power_monitor).lower())
element.text = str(self._power_monitor).lower()
def create_run_adjoint_subelement(self):
if not self._run_adjoint is None:
def _create_run_adjoint_subelement(self):
if self._run_adjoint is not None:
element = ET.SubElement(self._cmfd_file, "run_adjoint")
element.text = '{0}'.format(str(self._run_adjoint).lower())
element.text = str(self._run_adjoint).lower()
def _create_shift_subelement(self):
if self._shift is not None:
element = ET.SubElement(self._shift, "shift")
element.text = str(self._shift)
def create_write_matrices_subelement(self):
def _create_spectral_subelement(self):
if self._spectral is not None:
element = ET.SubElement(self._spectral, "spectral")
element.text = str(self._spectral)
if not self._write_matrices is None:
def _create_stol_subelement(self):
if self._stol is not None:
element = ET.SubElement(self._stol, "stol")
element.text = str(self._stol)
def _create_tally_reset_subelement(self):
if self._tally_reset is not None:
element = ET.SubElement(self._tally_reset, "tally_reset")
element.text = ' '.join(map(str, self._tally_reset))
def _create_write_matrices_subelement(self):
if self._write_matrices is not None:
element = ET.SubElement(self._cmfd_file, "write_matrices")
element.text = '{0}'.format(str(self._write_matrices).lower())
element.text = str(self._write_matrices).lower()
def export_to_xml(self):
"""Create a cmfd.xml file using the class data that can be used for an OpenMC
simulation.
self.create_active_flush_subelement()
self.create_begin_subelement()
self.create_display_subelement()
self.create_feedback_subelement()
self.create_inactive_subelement()
self.create_inactive_flush_subelement()
self.create_mesh_subelement()
self.create_norm_subelement()
self.create_num_flushes_subelement()
self.create_power_monitor_subelement()
self.create_run_adjoint_subelement()
self.create_write_matrices_subelement()
"""
self._create_begin_subelement()
self._create_dhat_reset_subelement()
self._create_display_subelement()
self._create_downscatter_subelement()
self._create_feedback_subelement()
self._create_gauss_seidel_tolerance_subelement()
self._create_ktol_subelement()
self._create_mesh_subelement()
self._create_norm_subelement()
self._create_power_monitor_subelement()
self._create_run_adjoint_subelement()
self._create_shift_subelement()
self._create_spectral_subelement()
self._create_stol_subelement()
self._create_tally_reset_subelement()
self._create_write_matrices_subelement()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._cmfd_file)

65
openmc/element.py
View file

@ -1,79 +1,78 @@
from openmc.checkvalue import *
import sys
from openmc.checkvalue import check_type
if sys.version_info[0] >= 3:
basestring = str
class Element(object):
"""A natural element used in a material via <element>. Internally, OpenMC will
expand the natural element into isotopes based on the known natural
abundances.
Parameters
----------
name : str
Chemical symbol of the element, e.g. Pu
xs : str
Cross section identifier, e.g. 71c
Attributes
----------
name : str
Chemical symbol of the element, e.g. Pu
xs : str
Cross section identifier, e.g. 71c
"""
def __init__(self, name='', xs=None):
# Initialize class attributes
self._name = ''
self._xs = None
# Set the Material class attributes
# Set class attributes
self.name = name
if not xs is None:
if xs is not None:
self.xs = xs
def __eq__(self, element2):
# Check type
if not isinstance(element2, Element):
return False
# Check name
# Check name and xs
if self._name != element2._name:
return False
# Check xs
elif self._xs != element2._xs:
return False
else:
return True
def __hash__(self):
hashable = []
hashable.append(self._name)
hashable.append(self._xs)
return hash(tuple(hashable))
return hash((self._name, self._xs))
@property
def xs(self):
return self._xs
@property
def name(self):
return self._name
@xs.setter
def xs(self, xs):
if not is_string(xs):
msg = 'Unable to set cross-section identifier xs for Element ' \
'with a non-string value {0}'.format(xs)
raise ValueError(msg)
check_type('cross section identifier', xs, basestring)
self._xs = xs
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Element with a non-string ' \
'value {0}'.format(name)
raise ValueError(msg)
check_type('name', name, basestring)
self._name = name
def __repr__(self):
string = 'Element - {0}\n'.format(self._name)
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
return string

95
openmc/executor.py
View file

@ -1,83 +1,128 @@
from __future__ import print_function
import subprocess
from numbers import Integral
import os
import sys
from openmc.checkvalue import *
from openmc.checkvalue import check_type
if sys.version_info[0] >= 3:
basestring = str
class Executor(object):
"""Control execution of OpenMC
Attributes
----------
working_directory : str
Path to working directory to run in
"""
def __init__(self):
self._working_directory = '.'
def _run_openmc(self, command, output):
# Launch a subprocess to run OpenMC
p = subprocess.Popen(command, shell=True,
p = subprocess.Popen(command, shell=True,
cwd=self._working_directory,
stdout=subprocess.PIPE)
# Capture and re-print OpenMC output in real-time
while (True and output):
line = p.stdout.readline()
print(line),
print(line, end='')
# If OpenMC is finished, break loop
if line == '' and p.poll() != None:
if not line and p.poll() != None:
break
# Return the returncode (integer, zero if no problems encountered)
return p.returncode
@property
def working_directory(self):
return self._working_directory
@working_directory.setter
def working_directory(self, working_directory):
if not is_string(working_directory):
msg = 'Unable to set Executor\'s working directory to {0} ' \
'since it is not a string'.format(working_directory)
raise ValueError(msg)
elif not os.path.isdir(working_directory):
check_type("Executor's working directory", working_directory,
basestring)
if not os.path.isdir(working_directory):
msg = 'Unable to set Executor\'s working directory to {0} ' \
'which does not exist'.format(working_directory)
raise ValueError(msg)
self._working_directory = working_directory
def plot_geometry(self, output=True, openmc_exec='openmc'):
"""Run OpenMC in plotting mode"""
def plot_geometry(self, output=True):
self._run_openmc('openmc -p', output)
return self._run_openmc(openmc_exec + ' -p', output)
def run_simulation(self, particles=None, threads=None,
geometry_debug=False, restart_file=None,
tracks=False, mpi_procs=1, output=True):
tracks=False, mpi_procs=1, output=True,
openmc_exec='openmc', mpi_exec=None):
"""Run an OpenMC simulation.
Parameters
----------
particles : int
Number of particles to simulate per generation
threads : int
Number of OpenMP threads
geometry_debug : bool
Turn on geometry debugging during simulation
restart_file : str
Path to restart file to use
tracks : bool
Write tracks for all particles
mpi_procs : int
Number of MPI processes
output : bool
Capture OpenMC output from standard out
openmc_exec : str
Path to OpenMC executable
"""
post_args = ' '
pre_args = ''
if is_integer(particles) and particles > 0:
if isinstance(particles, Integral) and particles > 0:
post_args += '-n {0} '.format(particles)
if is_integer(threads) and threads > 0:
if isinstance(threads, Integral) and threads > 0:
post_args += '-s {0} '.format(threads)
if geometry_debug:
post_args += '-g '
if is_string(restart_file):
if isinstance(restart_file, basestring):
post_args += '-r {0} '.format(restart_file)
if tracks:
post_args += '-t'
if is_integer(mpi_procs) and mpi_procs > 1:
pre_args += 'mpirun -n {0} '.format(mpi_procs)
if isinstance(mpi_procs, Integral) and mpi_procs > 1:
np_present = True
else:
np_present = False
command = pre_args + 'openmc ' + post_args
if mpi_exec is not None and isinstance(mpi_exec, basestring):
mpi_exec_present = True
else:
mpi_exec_present = False
self._run_openmc(command, output)
if np_present or mpi_exec_present:
if mpi_exec_present:
pre_args += mpi_exec + ' '
else:
pre_args += 'mpirun '
pre_args += '-n {0} '.format(mpi_procs)
command = pre_args + openmc_exec + ' ' + post_args
return self._run_openmc(command, output)

160
openmc/filter.py
View file

@ -1,15 +1,38 @@
from collections import Iterable
import copy
from numbers import Real, Integral
import numpy as np
from openmc import Mesh
from openmc.checkvalue import *
from openmc.constants import *
from openmc.checkvalue import check_type
class Filter(object):
"""A filter used to constrain a tally to a specific criterion, e.g. only tally
events when the particle is in a certain cell and energy range.
Parameters
----------
type : str
The type of the tally filter. Acceptable values are "universe",
"material", "cell", "cellborn", "surface", "mesh", "energy",
"energyout", and "distribcell".
bins : int or Iterable of int or Iterable of float
The bins for the filter. This takes on different meaning for different
filters.
Attributes
----------
type : str
The type of the tally filter.
bins : int or Iterable of int or Iterable of float
The bins for the filter
"""
# Initialize Filter class attributes
def __init__(self, type=None, bins=None):
self.type = type
self._num_bins = 0
self.bins = bins
@ -17,9 +40,7 @@ class Filter(object):
self._offset = -1
self._stride = None
def __eq__(self, filter2):
# Check type
if self._type != filter2._type:
return False
@ -35,21 +56,14 @@ class Filter(object):
else:
return True
def __hash__(self):
hashable = []
hashable.append(self._type)
hashable.append(self._bins)
return hash(tuple(hashable))
return hash((self._type, self._bins))
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)
@ -66,64 +80,52 @@ class Filter(object):
else:
return existing
@property
def type(self):
return self._type
@property
def bins(self):
return self._bins
@property
def num_bins(self):
return self._num_bins
@property
def mesh(self):
return self._mesh
@property
def offset(self):
return self._offset
@property
def stride(self):
return self._stride
@type.setter
def type(self, type):
if type is None:
self._type = type
elif not type in FILTER_TYPES.values():
elif type not in FILTER_TYPES.values():
msg = 'Unable to set Filter type to "{0}" since it is not one ' \
'of the supported types'.format(type)
raise ValueError(msg)
self._type = type
@bins.setter
def bins(self, bins):
if bins is None:
self.num_bins = 0
elif self._type is None:
msg = 'Unable to set bins for Filter to "{0}" since ' \
'the Filter type has not yet been set'.format(bins)
raise ValueError(msg)
# If the bin edge is a single value, it is a Cell, Material, etc. ID
if not isinstance(bins, (tuple, list, np.ndarray)):
if not isinstance(bins, Iterable):
bins = [bins]
# If the bins are in a collection, convert it to a list
@ -132,30 +134,23 @@ class Filter(object):
if self._type in ['cell', 'cellborn', 'surface', 'material',
'universe', 'distribcell']:
for edge in bins:
if not is_integer(edge):
if not isinstance(edge, Integral):
msg = 'Unable to add bin "{0}" to a {1} Filter since ' \
'it is a non-integer'.format(edge, self._type)
'it is not an integer'.format(edge, self._type)
raise ValueError(msg)
elif edge < 0:
msg = 'Unable to add bin "{0}" to a {1} Filter since ' \
'it is a negative integer'.format(edge, self._type)
'it is negative'.format(edge, self._type)
raise ValueError(msg)
elif self._type in ['energy', 'energyout']:
for edge in bins:
if not is_integer(edge) and not is_float(edge):
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, self._type)
raise ValueError(msg)
elif edge < 0.:
msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \
'since it is a negative value'.format(edge, self._type)
@ -163,27 +158,22 @@ class Filter(object):
# Check that bin edges are monotonically increasing
for index in range(len(bins)):
if index > 0 and bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a {1} Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, self._type)
raise ValueError(msg)
# mesh filters
elif self._type == 'mesh':
if not len(bins) == 1:
msg = 'Unable to add bins "{0}" to a mesh Filter since ' \
'only a single mesh can be used per tally'.format(bins)
raise ValueError(msg)
elif not is_integer(bins[0]):
elif not isinstance(bins[0], Integral):
msg = 'Unable to add bin "{0}" to mesh Filter since it ' \
'is a non-integer'.format(bins[0])
raise ValueError(msg)
elif bins[0] < 0:
msg = 'Unable to add bin "{0}" to mesh Filter since it ' \
'is a negative integer'.format(bins[0])
@ -192,12 +182,10 @@ class Filter(object):
# If all error checks passed, add bin edges
self._bins = bins
# FIXME
@num_bins.setter
def num_bins(self, num_bins):
if not is_integer(num_bins) or num_bins < 0:
if not isinstance(num_bins, Integral) or num_bins < 0:
msg = 'Unable to set the number of bins "{0}" for a {1} Filter ' \
'since it is not a positive ' \
'integer'.format(num_bins, self._type)
@ -205,39 +193,22 @@ class Filter(object):
self._num_bins = num_bins
@mesh.setter
def mesh(self, mesh):
if not isinstance(mesh, Mesh):
msg = 'Unable to set Mesh to "{0}" for Filter since it is not a ' \
'Mesh object'.format(mesh)
raise ValueError(msg)
check_type('filter mesh', mesh, Mesh)
self._mesh = mesh
self.type = 'mesh'
self.bins = self._mesh._id
@offset.setter
def offset(self, offset):
if not is_integer(offset):
msg = 'Unable to set offset "{0}" for a {1} Filter since it is a ' \
'non-integer value'.format(offset, self._type)
raise ValueError(msg)
check_type('filter offset', offset, Integral)
self._offset = offset
@stride.setter
def stride(self, stride):
if not is_integer(stride):
msg = 'Unable to set stride "{0}" for a {1} Filter since it is a ' \
'non-integer value'.format(stride, self._type)
raise ValueError(msg)
check_type('filter stride', stride, Integral)
if stride < 0:
msg = 'Unable to set stride "{0}" for a {1} Filter since it is a ' \
'negative value'.format(stride, self._type)
@ -245,8 +216,20 @@ class Filter(object):
self._stride = stride
def can_merge(self, filter):
"""Determine if filter can be merged with another.
Parameters
----------
filter : Filter
Filter to compare with
Returns
-------
bool
Whether the filter can be merged
"""
if not isinstance(filter, Filter):
return False
@ -270,8 +253,20 @@ class Filter(object):
else:
return True
def merge(self, filter):
"""Merge this filter with another.
Parameters
----------
filter : Filter
Filter to merge with
Returns
-------
merged_filter : Filter
Filter resulting from the merge
"""
if not self.can_merge(filter):
msg = 'Unable to merge {0} with {1} filters'.format(self._type, filter._type)
@ -287,33 +282,30 @@ class Filter(object):
return merged_filter
def get_bin_index(self, filter_bin):
"""Returns the index in the Filter for some bin.
Parameters
----------
filter_bin : int, tuple
The bin is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. The bin is an integer for
the cell instance ID for 'distribcell' Filters. The bin is
a 2-tuple of floats for 'energy' and 'energyout' filters
corresponding to the energy boundaries of the bin of interest.
The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding to
the mesh cell of interest.
filter_bin : int or tuple
The bin is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. The bin is an integer for the
cell instance ID for 'distribcell' Filters. The bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest.
Returns
-------
filter_index : int
The index in the Tally data array for this filter bin.
"""
# FIXME: This does not work for distribcells!!!
try:
# Filter bins for a mesh are an (x,y,z) tuple
if self.type == 'mesh':
# Convert (x,y,z) to a single bin -- this is similar to
# subroutine mesh_indices_to_bin in openmc/src/mesh.F90.
if (len(self.mesh.dimension) == 3):
@ -350,9 +342,7 @@ class Filter(object):
return filter_index
def __repr__(self):
string = 'Filter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self._bins)

125
openmc/geometry.py
View file

@ -2,7 +2,7 @@ from xml.etree import ElementTree as ET
import openmc
from openmc.clean_xml import *
from openmc.checkvalue import check_type
def reset_auto_ids():
openmc.reset_auto_material_id()
@ -12,28 +12,28 @@ def reset_auto_ids():
class Geometry(object):
"""Geometry representing a collection of surfaces, cells, and universes.
Attributes
----------
root_universe : openmc.universe.Universe
Root universe which contains all others
"""
def __init__(self):
# Initialize Geometry class attributes
self._root_universe = None
self._offsets = {}
@property
def root_universe(self):
return self._root_universe
@root_universe.setter
def root_universe(self, root_universe):
if not isinstance(root_universe, openmc.Universe):
msg = 'Unable to add root Universe {0} to Geometry since ' \
'it is not a Universe'.format(root_universe)
raise ValueError(msg)
elif root_universe._id != 0:
check_type('root universe', root_universe, openmc.Universe)
if root_universe._id != 0:
msg = 'Unable to add root Universe {0} to Geometry since ' \
'it has ID={1} instead of ' \
'ID=0'.format(root_universe, root_universe._id)
@ -41,24 +41,26 @@ class Geometry(object):
self._root_universe = root_universe
def get_offset(self, path, filter_offset):
"""
Returns the corresponding location in the results array for a given
path and filter number.
"""Returns the corresponding location in the results array for a given path and
filter number. This is primarily intended to post-processing result when
a distribcell filter is used.
Parameters
----------
path : list
A list of IDs that form the path to the target. It should begin
with 0 for the base universe, and should cover every universe,
cell, and lattice passed through. For the case of the lattice,
a tuple should be provided to indicate which coordinates in the
lattice should be entered. This should be in the
form: (lat_id, i_x, i_y, i_z)
A list of IDs that form the path to the target. It should begin with
0 for the base universe, and should cover every universe, cell, and
lattice passed through. For the case of the lattice, a tuple should
be provided to indicate which coordinates in the lattice should be
entered. This should be in the form: (lat_id, i_x, i_y, i_z)
filter_offset : int
An integer that specifies which offset map the filter is using
An integer that specifies which offset map the filter is using
Returns
-------
offset : int
Location in the results array for the path and filter
"""
@ -74,16 +76,39 @@ class Geometry(object):
# Return the final offset
return offset
def get_all_cells(self):
"""Return all cells defined
Returns
-------
list of openmc.universe.Cell
Cells in the geometry
"""
return self._root_universe.get_all_cells()
def get_all_universes(self):
"""Return all universes defined
Returns
-------
list of openmc.universe.Universe
Universes in the geometry
"""
return self._root_universe.get_all_universes()
def get_all_nuclides(self):
"""Return all nuclides assigned to a material in the geometry
Returns
-------
list of openmc.nuclide.Nuclide
Nuclides in the geometry
"""
nuclides = {}
materials = self.get_all_materials()
@ -93,8 +118,15 @@ class Geometry(object):
return nuclides
def get_all_materials(self):
"""Return all materials assigned to a cell
Returns
-------
list of openmc.material.Material
Materials in the geometry
"""
material_cells = self.get_all_material_cells()
materials = set()
@ -104,9 +136,7 @@ class Geometry(object):
return list(materials)
def get_all_material_cells(self):
all_cells = self.get_all_cells()
material_cells = set()
@ -116,16 +146,21 @@ class Geometry(object):
return list(material_cells)
def get_all_material_universes(self):
"""Return all universes composed of at least one non-fill cell
Returns
-------
list of openmc.universe.Universe
Universes with non-fill cells
"""
all_universes = self.get_all_universes()
material_universes = set()
for universe_id, universe in all_universes.items():
cells = universe._cells
for cell_id, cell in cells.items():
if cell._type == 'normal':
material_universes.add(universe)
@ -133,33 +168,35 @@ class Geometry(object):
return list(material_universes)
class GeometryFile(object):
"""Geometry file used for an OpenMC simulation. Corresponds directly to the
geometry.xml input file.
Attributes
----------
geometry : Geometry
The geometry to be used
"""
def __init__(self):
# Initialize GeometryFile class attributes
self._geometry = None
self._geometry_file = ET.Element("geometry")
@property
def geometry(self):
return self._geometry
@geometry.setter
def geometry(self, geometry):
if not isinstance(geometry, Geometry):
msg = 'Unable to set the Geometry to {0} for the GeometryFile ' \
'since it is not a Geometry object'.format(geometry)
raise ValueError(msg)
check_type('the geometry', geometry, Geometry)
self._geometry = geometry
def export_to_xml(self):
"""Create a geometry.xml file that can be used for a simulation.
"""
root_universe = self._geometry._root_universe
root_universe.create_xml_subelement(self._geometry_file)
@ -168,7 +205,7 @@ class GeometryFile(object):
sort_xml_elements(self._geometry_file)
clean_xml_indentation(self._geometry_file)
# Write the XML Tree to the materials.xml file
# Write the XML Tree to the geometry.xml file
tree = ET.ElementTree(self._geometry_file)
tree.write("geometry.xml", xml_declaration=True,
encoding='utf-8', method="xml")

305
openmc/material.py
View file

@ -1,10 +1,14 @@
from collections import MappingView
from collections import Iterable
from copy import deepcopy
from numbers import Real, Integral
import warnings
from xml.etree import ElementTree as ET
import sys
if sys.version_info[0] >= 3:
basestring = str
import openmc
from openmc.checkvalue import *
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.clean_xml import *
@ -14,6 +18,7 @@ MATERIAL_IDS = []
# A static variable for auto-generated Material IDs
AUTO_MATERIAL_ID = 10000
def reset_auto_material_id():
global AUTO_MATERIAL_ID, MATERIAL_IDS
AUTO_MATERIAL_ID = 10000
@ -23,20 +28,36 @@ def reset_auto_material_id():
# Units for density supported by OpenMC
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'at/b-cm', 'at/cm3', 'sum']
# ENDF temperatures
ENDF_TEMPS = np.array([300, 600, 700, 900, 1200, 1500])
# ENDF ZAIDs
ENDF_ZAIDS = np.array(['70c', '71c', '72c', '73c', '74c'])
# Constant for density when not needed
NO_DENSITY = 99999.
class Material(object):
"""A material composed of a collection of nuclides/elements that can be assigned
to a region of space.
Parameters
----------
material_id : int, optional
Unique identifier for the material. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the material. If not specified, the name will be the empty
string.
Attributes
----------
id : int
Unique identifier for the material
density : float
Density of the material (units defined separately)
density_units : str
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/cm3',
'atom/b-cm', 'atom/cm3', or 'sum'.
"""
def __init__(self, material_id=None, name=''):
# Initialize class attributes
self.id = material_id
self.name = name
@ -62,95 +83,75 @@ class Material(object):
# If specified, this file will be used instead of composition values
self._distrib_otf_file = None
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def density(self):
return self._density
@property
def density_units(self):
return self._density_units
@property
def convert_to_distrib_comps(self):
return self._convert_to_distrib_comps
@property
def distrib_otf_file(self):
return self._distrib_otf_file
@id.setter
def id(self, material_id):
global AUTO_MATERIAL_ID, MATERIAL_IDS
# If the Material already has an ID, remove it from global list
if hasattr(self, '_id') and not self._id is None:
if hasattr(self, '_id') and self._id is not None:
MATERIAL_IDS.remove(self._id)
if material_id is None:
self._id = AUTO_MATERIAL_ID
MATERIAL_IDS.append(AUTO_MATERIAL_ID)
AUTO_MATERIAL_ID += 1
# Check that the ID is an integer and wasn't already used
elif not is_integer(material_id):
msg = 'Unable to set a non-integer Material ' \
'ID {0}'.format(material_id)
raise ValueError(msg)
elif material_id in MATERIAL_IDS:
msg = 'Unable to set Material ID to {0} since a Material with ' \
'this ID was already initialized'.format(material_id)
raise ValueError(msg)
elif material_id < 0:
msg = 'Unable to set Material ID to {0} since it must be a ' \
'non-negative integer'.format(material_id)
raise ValueError(msg)
else:
check_type('material ID', material_id, Integral)
if material_id in MATERIAL_IDS:
msg = 'Unable to set Material ID to {0} since a Material with ' \
'this ID was already initialized'.format(material_id)
raise ValueError(msg)
check_greater_than('material ID', material_id, 0)
self._id = material_id
MATERIAL_IDS.append(material_id)
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Material ID={0} with a non-string ' \
'value {1}'.format(self._id, name)
raise ValueError(msg)
else:
self._name = name
check_type('name for Material ID={0}'.format(self._id),
name, basestring)
self._name = name
def set_density(self, units, density=NO_DENSITY):
"""Set the density of the material
if not is_float(density):
msg = 'Unable to set the density for Material ID={0} to a ' \
'non-floating point value {1}'.format(self._id, density)
raise ValueError(msg)
Parameters
----------
units : str
Physical units of density
density : float, optional
Value of the density. Must be specified unless units is given as
'sum'.
elif not units in DENSITY_UNITS:
msg = 'Unable to set the density for Material ID={0} with ' \
'units {1}'.format(self._id, units)
raise ValueError(msg)
"""
check_type('the density for Material ID={0}'.format(self._id),
density, Real)
check_value('density units', units, DENSITY_UNITS)
if density == NO_DENSITY and units is not 'sum':
msg = 'Unable to set the density Material ID={0} ' \
@ -161,45 +162,52 @@ class Material(object):
self._density = density
self._density_units = units
@distrib_otf_file.setter
def distrib_otf_file(self, filename):
# TODO: remove this when distributed materials are merged
warnings.warn('This feature is not yet implemented in a release ' \
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
if not is_string(filename) and not filename is None:
if not isinstance(filename, basestring) and filename is not None:
msg = 'Unable to add OTF material file to Material ID={0} with a ' \
'non-string name {1}'.format(self._id, filename)
raise ValueError(msg)
self._distrib_otf_file = filename
@convert_to_distrib_comps.setter
def convert_to_distrib_comps(self):
# TODO: remove this when distributed materials are merged
warnings.warn('This feature is not yet implemented in a release ' \
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
self._convert_to_distrib_comps = True
def add_nuclide(self, nuclide, percent, percent_type='ao'):
"""Add a nuclide to the material
Parameters
----------
nuclide : str or openmc.nuclide.Nuclide
Nuclide to add
percent : float
Atom or weight percent
percent_type : str
'ao' for atom percent and 'wo' for weight percent
"""
if not isinstance(nuclide, (openmc.Nuclide, str)):
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'non-Nuclide value {1}'.format(self._id, nuclide)
raise ValueError(msg)
elif not is_float(percent):
elif not isinstance(percent, Real):
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
raise ValueError(msg)
elif not percent_type in ['ao', 'wo', 'at/g-cm']:
elif percent_type not in ['ao', 'wo', 'at/g-cm']:
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
raise ValueError(msg)
@ -213,8 +221,15 @@ class Material(object):
self._nuclides[nuclide._name] = (nuclide, percent, percent_type)
def remove_nuclide(self, nuclide):
"""Remove a nuclide from the material
Parameters
----------
nuclide : openmc.nuclide.Nuclide
Nuclide to remove
"""
if not isinstance(nuclide, openmc.Nuclide):
msg = 'Unable to remove a Nuclide {0} in Material ID={1} ' \
@ -225,20 +240,31 @@ class Material(object):
if nuclide._name in self._nuclides:
del self._nuclides[nuclide._name]
def add_element(self, element, percent, percent_type='ao'):
"""Add a natural element to the material
Parameters
----------
element : openmc.element.Element
Element to add
percent : float
Atom or weight percent
percent_type : str
'ao' for atom percent and 'wo' for weight percent
"""
if not isinstance(element, openmc.Element):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-Element value {1}'.format(self._id, element)
raise ValueError(msg)
if not is_float(percent):
if not isinstance(percent, Real):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
raise ValueError(msg)
if not percent_type in ['ao', 'wo']:
if percent_type not in ['ao', 'wo']:
msg = 'Unable to add an Element to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
raise ValueError(msg)
@ -248,36 +274,58 @@ class Material(object):
self._elements[element._name] = (element, percent, percent_type)
def remove_element(self, element):
"""Remove a natural element from the material
Parameters
----------
element : openmc.element.Element
Element to remove
"""
# If the Material contains the Element, delete it
if element._name in self._elements:
del self._elements[element._name]
def add_s_alpha_beta(self, name, xs):
r"""Add an :math:`S(\alpha,\beta)` table to the material
if not is_string(name):
Parameters
----------
name : str
Name of the :math:`S(\alpha,\beta)` table
xs : str
Cross section identifier, e.g. '71t'
"""
if not isinstance(name, basestring):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string table name {1}'.format(self._id, name)
raise ValueError(msg)
if not is_string(xs):
if not isinstance(xs, basestring):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string cross-section identifier {1}'.format(self._id, xs)
raise ValueError(msg)
self._sab.append((name, xs))
def make_isotropic_in_lab(self):
for nuclide_name in self._nuclides:
self._nuclides[nuclide_name][0].make_isotropic_in_lab()
def get_all_nuclides(self):
"""Returns all nuclides in the material
Returns
-------
nuclides : dict
Dictionary whose keys are nuclide names and values are 2-tuples of
(nuclide, density)
"""
nuclides = {}
@ -288,9 +336,7 @@ class Material(object):
return nuclides
def _repr__(self):
string = 'Material\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
@ -322,9 +368,7 @@ class Material(object):
return string
def get_nuclide_xml(self, nuclide, distrib=False):
def _get_nuclide_xml(self, nuclide, distrib=False):
xml_element = ET.Element("nuclide")
xml_element.set("name", nuclide[0]._name)
@ -334,7 +378,7 @@ class Material(object):
else:
xml_element.set("wo", str(nuclide[1]))
if not nuclide[0]._xs is None:
if nuclide[0]._xs is not None:
xml_element.set("xs", nuclide[0]._xs)
if not nuclide[0]._scattering is None:
@ -342,9 +386,7 @@ class Material(object):
return xml_element
def get_element_xml(self, element, distrib=False):
def _get_element_xml(self, element, distrib=False):
xml_element = ET.Element("element")
xml_element.set("name", str(element[0]._name))
@ -356,28 +398,31 @@ class Material(object):
return xml_element
def get_nuclides_xml(self, nuclides, distrib=False):
def _get_nuclides_xml(self, nuclides, distrib=False):
xml_elements = []
for nuclide in nuclides.values():
xml_elements.append(self.get_nuclide_xml(nuclide, distrib))
xml_elements.append(self._get_nuclide_xml(nuclide, distrib))
return xml_elements
def get_elements_xml(self, elements, distrib=False):
def _get_elements_xml(self, elements, distrib=False):
xml_elements = []
for element in elements.values():
xml_elements.append(self.get_element_xml(element, distrib))
xml_elements.append(self._get_element_xml(element, distrib))
return xml_elements
def get_material_xml(self):
"""Return XML representation of the material
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing material data
"""
# Create Material XML element
element = ET.Element("material")
@ -393,19 +438,17 @@ class Material(object):
subelement.set("units", self._density_units)
if not self._convert_to_distrib_comps:
# Create nuclide XML subelements
subelements = self.get_nuclides_xml(self._nuclides)
subelements = self._get_nuclides_xml(self._nuclides)
for subelement in subelements:
element.append(subelement)
# Create element XML subelements
subelements = self.get_elements_xml(self._elements)
subelements = self._get_elements_xml(self._elements)
for subelement in subelements:
element.append(subelement)
else:
subelement = ET.SubElement(element, "compositions")
comps = []
@ -418,29 +461,24 @@ class Material(object):
raise ValueError(msg)
comps.append(per)
if self._distrib_otf_file is None:
# Create values and units subelements
subsubelement = ET.SubElement(subelement, "values")
subsubelement.text = ' '.join([str(c) for c in comps])
subsubelement = ET.SubElement(subelement, "units")
subsubelement.text = dist_per_type
else:
# Specify the materials file
subsubelement = ET.SubElement(subelement, "otf_file_path")
subsubelement.text = self._distrib_otf_file
# Create nuclide XML subelements
subelements = self.get_nuclides_xml(self._nuclides, distrib=True)
for subelement_nuc in subelements:
subelement.append(subelement_nuc)
# Create element XML subelements
subelements = self.get_elements_xml(self._elements, distrib=True)
subelements = self._get_elements_xml(self._elements, distrib=True)
for subelement_ele in subelements:
subelement.append(subelement_ele)
@ -454,31 +492,41 @@ class Material(object):
class MaterialsFile(object):
"""Materials file used for an OpenMC simulation. Corresponds directly to the
materials.xml input file.
Attributes
----------
default_xs : str
The default cross section identifier applied to a nuclide when none is
specified
"""
def __init__(self):
# Initialize MaterialsFile class attributes
self._materials = []
self._default_xs = None
self._materials_file = ET.Element("materials")
@property
def default_xs(self):
return self._default_xs
@default_xs.setter
def default_xs(self, xs):
if not is_string(xs):
msg = 'Unable to set default xs to a non-string value'.format(xs)
raise ValueError(msg)
check_type('default xs', xs, basestring)
self._default_xs = xs
def add_material(self, material):
"""Add a material to the file.
Parameters
----------
material : Material
Material to add
"""
if not isinstance(material, Material):
msg = 'Unable to add a non-Material {0} to the ' \
@ -487,19 +535,33 @@ class MaterialsFile(object):
self._materials.append(material)
def add_materials(self, materials):
"""Add multiple materials to the file.
if not isinstance(materials, (tuple, list, MappingView)):
Parameters
----------
materials : tuple or list of Material
Materials to add
"""
if not isinstance(materials, Iterable):
msg = 'Unable to create OpenMC materials.xml file from {0} which ' \
'is not a Python tuple/list'.format(materials)
'is not iterable'.format(materials)
raise ValueError(msg)
for material in materials:
self.add_material(material)
def remove_material(self, material):
"""Remove a material from the file
def remove_materials(self, material):
Parameters
----------
material : Material
Material to remove
"""
if not isinstance(material, Material):
msg = 'Unable to remove a non-Material {0} from the ' \
@ -509,26 +571,25 @@ class MaterialsFile(object):
self._materials.remove(material)
def make_isotropic_in_lab(self):
for material in self._materials:
materials.make_isotropic_in_lab()
def create_material_subelements(self):
def _create_material_subelements(self):
subelement = ET.SubElement(self._materials_file, "default_xs")
if not self._default_xs is None:
if self._default_xs is not None:
subelement.text = self._default_xs
for material in self._materials:
xml_element = material.get_material_xml()
self._materials_file.append(xml_element)
def export_to_xml(self):
"""Create a materials.xml file that can be used for a simulation.
self.create_material_subelements()
"""
self._create_material_subelements()
# Clean the indentation in the file to be user-readable
sort_xml_elements(self._materials_file)

252
openmc/mesh.py
View file

@ -1,8 +1,14 @@
from collections import Iterable
import copy
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
from openmc.checkvalue import *
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than)
if sys.version_info[0] >= 3:
basestring = str
# "Static" variable for auto-generated and Mesh IDs
AUTO_MESH_ID = 10000
@ -14,9 +20,37 @@ def reset_auto_mesh_id():
class Mesh(object):
"""A structured Cartesian mesh in two or three dimensions
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
type : str
Type of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
lower_left : Iterable of float
The lower-left corner of the structured mesh. If only two coordinate are
given, it is assumed that the mesh is an x-y mesh.
upper_right : Iterable of float
The upper-right corner of the structrued mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
width : Iterable of float
The width of mesh cells in each direction.
"""
def __init__(self, mesh_id=None, name=''):
# Initialize Mesh class attributes
self.id = mesh_id
self.name = name
@ -26,9 +60,7 @@ class Mesh(object):
self._upper_right = None
self._width = None
def __eq__(self, mesh2):
# Check type
if self._type != mesh2._type:
return False
@ -49,14 +81,11 @@ 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
@ -74,201 +103,87 @@ class Mesh(object):
else:
return existing
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def type(self):
return self._type
@property
def dimension(self):
return self._dimension
@property
def lower_left(self):
return self._lower_left
@property
def upper_right(self):
return self._upper_right
@property
def width(self):
return self._width
@property
def num_mesh_cells(self):
return np.prod(self._dimension)
@id.setter
def id(self, mesh_id):
if mesh_id is None:
global AUTO_MESH_ID
self._id = AUTO_MESH_ID
AUTO_MESH_ID += 1
# Check that the ID is an integer and wasn't already used
elif not is_integer(mesh_id):
msg = 'Unable to set a non-integer Mesh ID {0}'.format(mesh_id)
raise ValueError(msg)
elif mesh_id < 0:
msg = 'Unable to set Mesh ID to {0} since it must be a ' \
'non-negative integer'.format(mesh_id)
raise ValueError(msg)
else:
check_type('mesh ID', mesh_id, Integral)
check_greater_than('mesh ID', mesh_id, 0)
self._id = mesh_id
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Mesh ID={0} with a non-string ' \
'value {1}'.format(self._id, name)
raise ValueError(msg)
else:
self._name = name
check_type('name for mesh ID={0}'.format(self._id), name, basestring)
self._name = name
@type.setter
def type(self, type):
if not is_string(type):
msg = 'Unable to set Mesh ID={0} for type {1} which is not ' \
'a string'.format(self._id, type)
raise ValueError(msg)
elif not type in ['rectangular', 'hexagonal']:
msg = 'Unable to set Mesh ID={0} for type {1} which since ' \
'only rectangular and hexagonal meshes are ' \
'supported '.format(self._id, type)
raise ValueError(msg)
self._type = type
def type(self, meshtype):
check_type('type for mesh ID={0}'.format(self._id),
meshtype, basestring)
check_value('type for mesh ID={0}'.format(self._id),
meshtype, ['rectangular', 'hexagonal'])
self._type = meshtype
@dimension.setter
def dimension(self, dimension):
if not isinstance(dimension, (tuple, list, np.ndarray)):
msg = 'Unable to set Mesh ID={0} with dimension {1} which is ' \
'not a Python list, tuple or NumPy ' \
'array'.format(self._id, dimension)
raise ValueError(msg)
elif len(dimension) != 2 and len(dimension) != 3:
msg = 'Unable to set Mesh ID={0} with dimension {1} since it ' \
'must include 2 or 3 dimensions'.format(self._id, dimension)
raise ValueError(msg)
for dim in dimension:
if not is_integer(dim):
msg = 'Unable to set Mesh ID={0} with dimension {1} which ' \
'is a non-integer'.format(self._id, dim)
raise ValueError(msg)
check_type('mesh dimension', dimension, Iterable, Integral)
check_length('mesh dimension', dimension, 2, 3)
self._dimension = dimension
@lower_left.setter
def lower_left(self, lower_left):
if not isinstance(lower_left, (tuple, list, np.ndarray)):
msg = 'Unable to set Mesh ID={0} with lower_left {1} which is ' \
'not a Python list, tuple or NumPy ' \
'array'.format(self._id, lower_left)
raise ValueError(msg)
elif len(lower_left) != 2 and len(lower_left) != 3:
msg = 'Unable to set Mesh ID={0} with lower_left {1} since it ' \
'must include 2 or 3 dimensions'.format(self._id, lower_left)
raise ValueError(msg)
for coord in lower_left:
if not is_integer(coord) and not is_float(coord):
msg = 'Unable to set Mesh ID={0} with lower_left {1} which ' \
'is neither neither an integer nor a floating point ' \
'value'.format(self._id, coord)
raise ValueError(msg)
check_type('mesh lower_left', lower_left, Iterable, Real)
check_length('mesh lower_left', lower_left, 2, 3)
self._lower_left = lower_left
@upper_right.setter
def upper_right(self, upper_right):
if not isinstance(upper_right, (tuple, list, np.ndarray)):
msg = 'Unable to set Mesh ID={0} with upper_right {1} which ' \
'is not a Python list, tuple or NumPy ' \
'array'.format(self._id, upper_right)
raise ValueError(msg)
if len(upper_right) != 2 and len(upper_right) != 3:
msg = 'Unable to set Mesh ID={0} with upper_right {1} since it ' \
'must include 2 or 3 dimensions'.format(self._id, upper_right)
raise ValueError(msg)
for coord in upper_right:
if not is_integer(coord) and not is_float(coord):
msg = 'Unable to set Mesh ID={0} with upper_right {1} which ' \
'is neither an integer nor a floating point ' \
'value'.format(self._id, coord)
raise ValueError(msg)
check_type('mesh upper_right', upper_right, Iterable, Real)
check_length('mesh upper_right', upper_right, 2, 3)
self._upper_right = upper_right
@width.setter
def width(self, width):
if not width is None:
if not isinstance(width, (tuple, list, np.ndarray)):
msg = 'Unable to set Mesh ID={0} with width {1} which ' \
'is not a Python list, tuple or NumPy ' \
'array'.format(self._id, width)
raise ValueError(msg)
if len(width) != 2 and len(width) != 3:
msg = 'Unable to set Mesh ID={0} with width {1} since it must ' \
'include 2 or 3 dimensions'.format(self._id, width)
raise ValueError(msg)
for dim in width:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to set Mesh ID={0} with width {1} which is ' \
'neither an integer nor a floating point ' \
'value'.format(self._id, width)
raise ValueError(msg)
check_type('mesh width', width, Iterable, Real)
check_length('mesh width', width, 2, 3)
self._width = width
def __repr__(self):
string = 'Mesh\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
@ -279,53 +194,32 @@ class Mesh(object):
string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._width)
return string
def get_mesh_xml(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", self._type)
if len(self._dimension) == 2:
subelement = ET.SubElement(element, "dimension")
subelement.text = '{0} {1}'.format(self._dimension[0],
self._dimension[1])
else:
subelement = ET.SubElement(element, "dimension")
subelement.text = '{0} {1} {2}'.format(self._dimension[0],
self._dimension[1],
self._dimension[2])
subelement = ET.SubElement(element, "dimension")
subelement.text = ' '.join(map(str, self._dimension))
if len(self._lower_left) == 2:
subelement = ET.SubElement(element, "lower_left")
subelement.text = '{0} {1}'.format(self._lower_left[0],
self._lower_left[1])
else:
subelement = ET.SubElement(element, "lower_left")
subelement.text = '{0} {1} {2}'.format(self._lower_left[0],
self._lower_left[1],
self._lower_left[2])
subelement = ET.SubElement(element, "lower_left")
subelement.text = ' '.join(map(str, self._lower_left))
if not self._upper_right is None:
if len(self._upper_right) == 2:
subelement = ET.SubElement(element, "upper_right")
subelement.text = '{0} {1}'.format(self._upper_right[0],
self._upper_right[1])
else:
subelement = ET.SubElement(element, "upper_right")
subelement.text = '{0} {1} {2}'.format(self._upper_right[0],
self._upper_right[1],
self._upper_right[2])
if self._upper_right is not None:
subelement = ET.SubElement(element, "upper_right")
subelement.text = ' '.join(map(str, self._upper_right))
if not self._width is None:
if len(self._width) == 2:
subelement = ET.SubElement(element, "width")
subelement.text = '{0} {1}'.format(self._width[0],
self._width[1])
else:
subelement = ET.SubElement(element, "width")
subelement.text = '{0} {1} {2}'.format(self._width[0],
self._width[1],
self._width[2])
if self._width is not None:
subelement = ET.SubElement(element, "width")
subelement.text = ' '.join(map(str, self._width))
return element
return element

69
openmc/nuclide.py
View file

@ -1,10 +1,35 @@
from openmc.checkvalue import *
from numbers import Integral
import sys
from openmc.checkvalue import check_type
if sys.version_info[0] >= 3:
basestring = str
class Nuclide(object):
"""A nuclide that can be used in a material.
Parameters
----------
name : str
Name of the nuclide, e.g. U-235
xs : str
Cross section identifier, e.g. 71c
Attributes
----------
name : str
Name of the nuclide, e.g. U-235
xs : str
Cross section identifier, e.g. 71c
zaid : int
1000*(atomic number) + mass number. As an example, the zaid of U-235
would be 92235.
"""
def __init__(self, name='', xs=None):
# Initialize class attributes
self._name = ''
self._xs = None
@ -14,12 +39,10 @@ class Nuclide(object):
# Set the Material class attributes
self.name = name
if not xs is None:
if xs is not None:
self.xs = xs
def __eq__(self, nuclide2):
# Check type
if not isinstance(nuclide2, Nuclide):
return False
@ -35,63 +58,41 @@ class Nuclide(object):
else:
return True
def __hash__(self):
return hash((self._name, self._xs))
@property
def name(self):
return self._name
@property
def xs(self):
return self._xs
@property
def zaid(self):
return self._zaid
@property
def scattering(self):
return self._scattering
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Nuclide with a non-string ' \
'value {0}'.format(name)
raise ValueError(msg)
check_type('name', name, basestring)
self._name = name
@xs.setter
def xs(self, xs):
if not is_string(xs):
msg = 'Unable to set cross-section identifier xs for Nuclide ' \
'with a non-string value {0}'.format(xs)
raise ValueError(msg)
check_type('cross-section identifier', xs, basestring)
self._xs = xs
@zaid.setter
def zaid(self, zaid):
if not is_integer(zaid):
msg = 'Unable to set zaid for Nuclide ' \
'with a non-integer {0}'.format(zaid)
raise ValueError(msg)
check_type('zaid', zaid, Integral)
self._zaid = zaid
<<<<<<< HEAD
@scattering.setter
def scattering(self, scattering):
@ -104,13 +105,17 @@ class Nuclide(object):
self._scattering = scattering
=======
>>>>>>> develop
def __repr__(self):
string = 'Nuclide - {0}\n'.format(self._name)
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
if self._zaid is not None:
string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid)
<<<<<<< HEAD
if self._scattering is not None:
string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t',
self._scattering)
=======
>>>>>>> develop
return string

249
openmc/opencg_compatible.py
View file

@ -62,8 +62,20 @@ OPENMC_LATTICES = {}
OPENCG_LATTICES = {}
def get_opencg_material(openmc_material):
"""Return an OpenCG material corresponding to an OpenMC material.
Parameters
----------
openmc_material : openmc.material.Material
OpenMC material
Returns
-------
opencg_material : opencg.Material
Equivalent OpenCG material
"""
if not isinstance(openmc_material, openmc.Material):
msg = 'Unable to create an OpenCG Material from {0} ' \
@ -91,6 +103,19 @@ def get_opencg_material(openmc_material):
def get_openmc_material(opencg_material):
"""Return an OpenMC material corresponding to an OpenCG material.
Parameters
----------
opencg_material : opencg.Material
OpenCG material
Returns
-------
openmc_material : openmc.material.Material
Equivalent OpenMC material
"""
if not isinstance(opencg_material, opencg.Material):
msg = 'Unable to create an OpenMC Material from {0} ' \
@ -118,6 +143,25 @@ def get_openmc_material(opencg_material):
def is_opencg_surface_compatible(opencg_surface):
"""Determine whether OpenCG surface is compatible with OpenMC geometry.
A surface is considered compatible if there is a one-to-one correspondence
between OpenMC and OpenCG surface types. Note that some OpenCG surfaces,
e.g. SquarePrism, do not have a one-to-one correspondence with OpenMC
surfaces but can still be converted into an equivalent collection of OpenMC
surfaces.
Parameters
----------
opencg_surface : opencg.Surface
OpenCG surface
Returns
-------
bool
Whether OpenCG surface is compatible with OpenMC
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to check if OpenCG Surface is compatible' \
@ -132,6 +176,19 @@ def is_opencg_surface_compatible(opencg_surface):
def get_opencg_surface(openmc_surface):
"""Return an OpenCG surface corresponding to an OpenMC surface.
Parameters
----------
openmc_surface : openmc.surface.Surface
OpenMC surface
Returns
-------
opencg_surface : opencg.Surface
Equivalent OpenCG surface
"""
if not isinstance(openmc_surface, openmc.Surface):
msg = 'Unable to create an OpenCG Surface from {0} ' \
@ -205,6 +262,19 @@ def get_opencg_surface(openmc_surface):
def get_openmc_surface(opencg_surface):
"""Return an OpenMC surface corresponding to an OpenCG surface.
Parameters
----------
opencg_surface : opencg.Surface
OpenCG surface
Returns
-------
openmc_surface : openmc.surface.Surface
Equivalent OpenMC surface
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
@ -269,7 +339,6 @@ def get_openmc_surface(opencg_surface):
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
# Add the OpenMC Surface to the global collection of all OpenMC Surfaces
OPENMC_SURFACES[surface_id] = openmc_surface
@ -280,6 +349,23 @@ def get_openmc_surface(opencg_surface):
def get_compatible_opencg_surfaces(opencg_surface):
"""Generate OpenCG surfaces that are compatible with OpenMC equivalent to an
OpenCG surface that is not compatible. For example, this method may be used
to convert a ZSquarePrism OpenCG surface into a collection of equivalent
XPlane and YPlane OpenCG surfaces.
Parameters
----------
opencg_surface : opencg.Surface
OpenCG surface that is incompatible with OpenMC
Returns
-------
surfaces : list of opencg.Surface
Collection of surfaces equivalent to the original one but compatible
with OpenMC
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
@ -349,6 +435,19 @@ def get_compatible_opencg_surfaces(opencg_surface):
def get_opencg_cell(openmc_cell):
"""Return an OpenCG cell corresponding to an OpenMC cell.
Parameters
----------
openmc_cell : openmc.universe.Cell
OpenMC cell
Returns
-------
opencg_cell : opencg.Cell
Equivalent OpenCG cell
"""
if not isinstance(openmc_cell, openmc.Cell):
msg = 'Unable to create an OpenCG Cell from {0} which ' \
@ -398,7 +497,26 @@ def get_opencg_cell(openmc_cell):
def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
"""Generate OpenCG cells that are compatible with OpenMC equivalent to an OpenCG
cell that is not compatible.
Parameters
----------
opencg_cell : opencg.Cell
OpenCG cell
opencg_surface : opencg.Surface
OpenCG surface that causes the incompatibility, e.g. an instance of
XSquarePrism
halfspace : {-1, 1}
Which halfspace defined by the surface is contained in the cell
Returns
-------
compatible_cells : list of opencg.Cell
Collection of cells equivalent to the original one but compatible with
OpenMC
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create compatible OpenMC Cell from {0} which ' \
'is not an OpenCG Cell'.format(opencg_cell)
@ -409,7 +527,7 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
'not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
elif not halfspace in [-1, +1]:
elif halfspace not in [-1, +1]:
msg = 'Unable to create compatible Cell since {0}' \
'is not a +/-1 halfspace'.format(halfspace)
raise ValueError(msg)
@ -418,8 +536,8 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
compatible_cells = []
# SquarePrism Surfaces
if opencg_surface._type in ['x-squareprism',
'y-squareprism', 'z-squareprism']:
if opencg_surface._type in ['x-squareprism', 'y-squareprism',
'z-squareprism']:
# Get the compatible Surfaces (XPlanes and YPlanes)
compatible_surfaces = get_compatible_opencg_surfaces(opencg_surface)
@ -436,13 +554,11 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
# If Cell is outside SquarePrism, add "outside" of Surface halfspaces
else:
# Create 8 Cell clones to represent each of the disjoint planar
# Surface halfspace intersections
num_clones = 8
for clone_id in range(num_clones):
# Create a cloned OpenCG Cell with Surfaces compatible with OpenMC
clone = opencg_cell.clone()
compatible_cells.append(clone)
@ -500,6 +616,14 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
def make_opencg_cells_compatible(opencg_universe):
"""Make all cells in an OpenCG universe compatible with OpenMC.
Parameters
----------
opencg_universe : opencg.Universe
Universe to check
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to make compatible OpenCG Cells for {0} which ' \
@ -545,8 +669,20 @@ def make_opencg_cells_compatible(opencg_universe):
return
def get_openmc_cell(opencg_cell):
"""Return an OpenMC cell corresponding to an OpenCG cell.
Parameters
----------
opencg_cell : opencg.Cell
OpenCG cell
Returns
-------
openmc_cell : openmc.universe.Cell
Equivalent OpenMC cell
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create an OpenMC Cell from {0} which ' \
@ -597,8 +733,20 @@ def get_openmc_cell(opencg_cell):
return openmc_cell
def get_opencg_universe(openmc_universe):
"""Return an OpenCG universe corresponding to an OpenMC universe.
Parameters
----------
openmc_universe : openmc.universe.Universe
OpenMC universe
Returns
-------
opencg_universe : opencg.Universe
Equivalent OpenCG universe
"""
if not isinstance(openmc_universe, openmc.Universe):
msg = 'Unable to create an OpenCG Universe from {0} which ' \
@ -633,6 +781,19 @@ def get_opencg_universe(openmc_universe):
def get_openmc_universe(opencg_universe):
"""Return an OpenMC universe corresponding to an OpenCG universe.
Parameters
----------
opencg_universe : opencg.Universe
OpenCG universe
Returns
-------
openmc_universe : openmc.universe.Universe
Equivalent OpenMC universe
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to create an OpenMC Universe from {0} which ' \
@ -670,6 +831,19 @@ def get_openmc_universe(opencg_universe):
def get_opencg_lattice(openmc_lattice):
"""Return an OpenCG lattice corresponding to an OpenMC lattice.
Parameters
----------
openmc_lattice : openmc.universe.Lattice
OpenMC lattice
Returns
-------
opencg_lattice : opencg.Lattice
Equivalent OpenCG lattice
"""
if not isinstance(openmc_lattice, openmc.Lattice):
msg = 'Unable to create an OpenCG Lattice from {0} which ' \
@ -701,7 +875,7 @@ def get_opencg_lattice(openmc_lattice):
lower_left = new_lower_left
# Initialize an empty array for the OpenCG nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), \
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
dtype=opencg.Universe)
# Create OpenCG Universes for each unique nested Universe in this Lattice
@ -723,8 +897,8 @@ def get_opencg_lattice(openmc_lattice):
opencg_lattice.setUniverses(universe_array)
offset = np.array(lower_left, dtype=np.float64) - \
((np.array(pitch, dtype=np.float64) * \
np.array(dimension, dtype=np.float64))) / -2.0
((np.array(pitch, dtype=np.float64) *
np.array(dimension, dtype=np.float64))) / -2.0
opencg_lattice.setOffset(offset)
# Add the OpenMC Lattice to the global collection of all OpenMC Lattices
@ -737,6 +911,19 @@ def get_opencg_lattice(openmc_lattice):
def get_openmc_lattice(opencg_lattice):
"""Return an OpenMC lattice corresponding to an OpenCG lattice.
Parameters
----------
opencg_lattice : opencg.Lattice
OpenCG lattice
Returns
-------
openmc_lattice : openmc.universe.Lattice
Equivalent OpenMC lattice
"""
if not isinstance(opencg_lattice, opencg.Lattice):
msg = 'Unable to create an OpenMC Lattice from {0} which ' \
@ -756,7 +943,7 @@ def get_openmc_lattice(opencg_lattice):
universes = opencg_lattice._universes
# Initialize an empty array for the OpenMC nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)), \
universe_array = np.ndarray(tuple(np.array(dimension)),
dtype=openmc.Universe)
# Create OpenMC Universes for each unique nested Universe in this Lattice
@ -773,11 +960,11 @@ def get_openmc_lattice(opencg_lattice):
universe_array[x][y][z] = unique_universes[universe_id]
# Reverse y-dimension in array to match ordering in OpenCG
universe_array = universe_array[:,::-1,:]
universe_array = universe_array[:, ::-1, :]
lower_left = np.array(offset, dtype=np.float64) + \
((np.array(width, dtype=np.float64) * \
np.array(dimension, dtype=np.float64))) / -2.0
((np.array(width, dtype=np.float64) *
np.array(dimension, dtype=np.float64))) / -2.0
openmc_lattice = openmc.RectLattice(lattice_id=lattice_id)
openmc_lattice.dimension = dimension
@ -795,6 +982,19 @@ def get_openmc_lattice(opencg_lattice):
def get_opencg_geometry(openmc_geometry):
"""Return an OpenCG geometry corresponding to an OpenMC geometry.
Parameters
----------
openmc_geometry : openmc.universe.Geometry
OpenMC geometry
Returns
-------
opencg_geometry : opencg.Geometry
Equivalent OpenCG geometry
"""
if not isinstance(openmc_geometry, openmc.Geometry):
msg = 'Unable to get OpenCG geometry from {0} which is ' \
@ -822,6 +1022,19 @@ def get_opencg_geometry(openmc_geometry):
def get_openmc_geometry(opencg_geometry):
"""Return an OpenMC geometry corresponding to an OpenCG geometry.
Parameters
----------
opencg_geometry : opencg.Geometry
OpenCG geometry
Returns
-------
openmc_geometry : openmc.universe.Geometry
Equivalent OpenMC geometry
"""
if not isinstance(opencg_geometry, opencg.Geometry):
msg = 'Unable to get OpenMC geometry from {0} which is ' \
@ -849,8 +1062,8 @@ def get_openmc_geometry(opencg_geometry):
# Make the entire geometry "compatible" before assigning auto IDs
universes = opencg_geometry.getAllUniverses()
for universe_id, universe in universes.items():
if not isinstance(universe, opencg.Lattice):
make_opencg_cells_compatible(universe)
if not isinstance(universe, opencg.Lattice):
make_opencg_cells_compatible(universe)
opencg_geometry.assignAutoIds()

39
openmc/particle_restart.py
View file

@ -1,9 +1,44 @@
#!/usr/bin/env python
import struct
class Particle(object):
"""Information used to restart a specific particle that caused a simulation to
fail.
Parameters
----------
filename : str
Path to the particle restart file
Attributes
----------
filetype : int
Integer indicating the file type
revision : int
Revision of the particle restart format
current_batch : int
The batch containing the particle
gen_per_batch : int
Number of generations per batch
current_gen : int
The generation containing the particle
n_particles : int
Number of particles per generation
run_mode : int
Type of simulation (criticality or fixed source)
id : long
Identifier of the particle
weight : float
Weight of the particle
energy : float
Energy of the particle in MeV
xyz : list of float
Position of the particle
uvw : list of float
Directional cosines of the particle
"""
def __init__(self, filename):
if filename.endswith('.h5'):
import h5py

409
openmc/plots.py
View file

@ -1,14 +1,21 @@
from collections import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
from openmc.checkvalue import (check_type, check_value, check_length,
check_greater_than, check_less_than)
if sys.version_info[0] >= 3:
basestring = str
# A static variable for auto-generated Plot IDs
AUTO_PLOT_ID = 10000
def reset_auto_plot_id():
global AUTO_PLOT_ID
AUTO_PLOT_ID = 10000
@ -18,9 +25,50 @@ BASES = ['xy', 'xz', 'yz']
class Plot(object):
"""Definition of a finite region of space to be plotted, either as a slice plot
in two dimensions or as a voxel plot in three dimensions.
Parameters
----------
plot_id : int
Unique identifier for the plot
name : str
Name of the plot
Attributes
----------
id : int
Unique identifier
name : str
Name of the plot
width : Iterable of float
Width of the plot in each basis direction
pixels : Iterable of int
Number of pixels to use in each basis direction
origin : tuple or list of ndarray
Origin (center) of the plot
filename :
Path to write the plot to
color : {'cell', 'mat'}
Indicate whether the plot should be colored by cell or by material
type : {'slice', 'voxel'}
The type of the plot
basis : {'xy', 'xz', 'yz'}
The basis directions for the plot
background : tuple or list of ndarray
Color of the background defined by RGB
mask_components : Iterable of int
Unique id numbers of the cells or materials to plot
mask_background : Iterable of int
Color to apply to all cells/materials not listed in mask_components
defined by RGB
col_spec : dict
Dictionary indicating that certain cells/materials (keys) should be
colored with a specific RGB (values)
"""
def __init__(self, plot_id=None, name=''):
# Initialize Plot class attributes
self.id = plot_id
self.name = name
@ -36,295 +84,139 @@ class Plot(object):
self._mask_background = None
self._col_spec = None
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def width(self):
return self._width
@property
def pixels(self):
return self._pixels
@property
def origin(self):
return self._origin
@property
def filename(self):
return self._filename
@property
def color(self):
return self._color
@property
def type(self):
return self._type
@property
def basis(self):
return self._basis
@property
def background(self):
return self._background
@property
def mask_componenets(self):
return self._mask_components
@property
def mask_background(self):
return self._mask_background
@property
def col_spec(self):
return self._col_spec
@id.setter
def id(self, plot_id):
if plot_id is None:
global AUTO_PLOT_ID
self._id = AUTO_PLOT_ID
AUTO_PLOT_ID += 1
# Check that the ID is an integer and wasn't already used
elif not is_integer(plot_id):
msg = 'Unable to set a non-integer Plot ID {0}'.format(plot_id)
raise ValueError(msg)
elif plot_id < 0:
msg = 'Unable to set Plot ID to {0} since it must be a ' \
'non-negative integer'.format(plot_id)
raise ValueError(msg)
else:
check_type('plot ID', plot_id, Integral)
check_greater_than('plot ID', plot_id, 0)
self._id = plot_id
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Plot ID={0} with a non-string ' \
'value {1}'.format(self._id, name)
raise ValueError(msg)
else:
self._name = name
check_type('plot name', name, basestring)
self._name = name
@width.setter
def width(self, width):
if not isinstance(width, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with width {1} which is not ' \
'a Python tuple/list or NumPy array'.format(self._id, width)
raise ValueError(msg)
elif len(width) != 2 and len(width) != 3:
msg = 'Unable to create Plot ID={0} with width {1} since only 2D ' \
'and 3D plots are supported'.format(self._id, width)
raise ValueError(msg)
for dim in width:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to create Plot ID={0} with width {1} since ' \
'each element must be a floating point value or ' \
'integer'.format(self._id, width)
raise ValueError(msg)
check_type('plot width', width, Iterable, Real)
check_length('plot width', width, 2, 3)
self._width = width
@origin.setter
def origin(self, origin):
if not isinstance(origin, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with origin {1} which is not ' \
'a Python tuple/list or NumPy array'.format(self._id, origin)
raise ValueError(msg)
elif len(origin) != 3:
msg = 'Unable to create Plot ID={0} with origin {1} since only ' \
'a 3D coordinate must be input'.format(self._id, origin)
raise ValueError(msg)
for dim in origin:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to create Plot ID={0} with origin {1} since ' \
'each element must be a floating point value or ' \
'integer'.format(self._id, origin)
raise ValueError(msg)
check_type('plot origin', origin, Iterable, Real)
check_length('plot origin', origin, 3)
self._origin = origin
@pixels.setter
def pixels(self, pixels):
if not isinstance(pixels, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with pixels {1} which is not ' \
'a Python tuple/list or NumPy array'.format(self._id, pixels)
raise ValueError(msg)
elif len(pixels) != 2 and len(pixels) != 3:
msg = 'Unable to create Plot ID={0} with pixels {1} since ' \
'only 2D and 3D plots are supported'.format(self._id, pixels)
raise ValueError(msg)
check_type('plot pixels', pixels, Iterable, Integral)
check_length('plot pixels', pixels, 2, 3)
for dim in pixels:
if not is_integer(dim):
msg = 'Unable to create Plot ID={0} with pixel value {1} ' \
'which is not an integer'.format(self._id, dim)
raise ValueError(msg)
elif dim < 0:
msg = 'Unable to create Plot ID={0} with pixel value {1} ' \
'which is less than 0'.format(self._id, dim)
raise ValueError(msg)
check_greater_than('plot pixels', dim, 0)
self._pixels = pixels
@filename.setter
def filename(self, filename):
if not is_string(filename):
msg = 'Unable to create Plot ID={0} with filename {1} which is ' \
'not a string'.format(self._id, filename)
raise ValueError(msg)
check_type('filename', filename, basestring)
self._filename = filename
@color.setter
def color(self, color):
if not is_string(color):
msg = 'Unable to create Plot ID={0} with color {1} which is not ' \
'a string'.format(self._id, color)
raise ValueError(msg)
elif not color in ['cell', 'mat']:
msg = 'Unable to create Plot ID={0} with color {1} which is not ' \
'a cell or mat'.format(self._id, color)
raise ValueError(msg)
check_type('plot color', color, basestring)
check_value('plot color', color, ['cell', 'mat'])
self._color = color
@type.setter
def type(self, type):
if not is_string(type):
msg = 'Unable to create Plot ID={0} with type {1} which is not ' \
'a string'.format(self._id, type)
raise ValueError(msg)
elif not type in ['slice', 'voxel']:
msg = 'Unable to create Plot ID={0} with type {1} which is not ' \
'slice or voxel'.format(self._id, type)
raise ValueError(msg)
self._type = type
def type(self, plottype):
check_type('plot type', plottype, basestring)
check_value('plot type', plottype, ['slice', 'voxel'])
self._type = plottype
@basis.setter
def basis(self, basis):
if not is_string(basis):
msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \
'a string'.format(self._id, basis)
raise ValueError(msg)
elif not basis in ['xy', 'xz', 'yz']:
msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \
'xy, xz, or yz'.format(self._id, basis)
raise ValueError(msg)
check_type('plot basis', basis, basestring)
check_value('plot basis', basis, ['xy', 'xz', 'yz'])
self._basis = basis
@background.setter
def background(self, background):
if not isinstance(background, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with background {1} ' \
'which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, background)
raise ValueError(msg)
elif len(background) != 3:
msg = 'Unable to create Plot ID={0} with background {1} ' \
'which is not 3 integer RGB ' \
'values'.format(self._id, background)
raise ValueError(msg)
check_type('plot background', background, Iterable, Integral)
check_length('plot background', background, 3)
for rgb in background:
if not is_integer(rgb):
msg = 'Unable to create Plot ID={0} with background RGB ' \
'value {1} which is not an integer'.format(self._id, rgb)
raise ValueError(msg)
elif rgb < 0 or rgb > 255:
msg = 'Unable to create Plot ID={0} with background RGB value ' \
'{1} which is not between 0 and 255'.format(self._id, rgb)
raise ValueError(msg)
check_greater_than('plot background',rgb, 0, True)
check_less_than('plot background', rgb, 256)
self._background = background
@col_spec.setter
def col_spec(self, col_spec):
if not isinstance(col_spec, dict):
msg= 'Unable to create Plot ID={0} with col_spec parameter {1} ' \
'which is not a Python dictionary of IDs to ' \
'pixels'.format(self._id, col_spec)
raise ValueError(msg)
check_type('plot col_spec parameter', col_spec, dict, Integral)
for key in col_spec:
if not is_integer(key):
msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
'which is not an integer'.format(self._id, key)
raise ValueError(msg)
elif key < 0:
if key < 0:
msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
'which is less than 0'.format(self._id, key)
raise ValueError(msg)
elif not isinstance(col_spec[key], (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with col_spec RGB ' \
'values {1} which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, col_spec[key])
elif not isinstance(col_spec[key], Iterable):
msg = 'Unable to create Plot ID={0} with col_spec RGB values' \
' {1} which is not iterable'.format(self._id, col_spec[key])
raise ValueError(msg)
elif len(col_spec[key]) != 3:
@ -335,63 +227,23 @@ class Plot(object):
self._col_spec = col_spec
@mask_componenets.setter
def mask_components(self, mask_components):
if not isinstance(mask_components, (list, tuple, np.ndarray)):
msg = 'Unable to create Plot ID={0} with mask components {1} ' \
'which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, mask_components)
raise ValueError(msg)
check_type('plot mask_components', mask_components, Iterable, Integral)
for component in mask_components:
if not is_integer(component):
msg = 'Unable to create Plot ID={0} with mask component {1} ' \
'which is not an integer'.format(self._id, component)
raise ValueError(msg)
elif component < 0:
msg = 'Unable to create Plot ID={0} with mask component {1} ' \
'which is less than 0'.format(self._id, component)
raise ValueError(msg)
check_greater_than('plot mask_components', component, 0, True)
self._mask_components = mask_components
@mask_background.setter
def mask_background(self, mask_background):
if not isinstance(mask_background, (list, tuple, np.ndarray)):
msg = 'Unable to create Plot ID={0} with mask background {1} ' \
'which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, mask_background)
raise ValueError(msg)
elif len(mask_background) != 3 and len(mask_background) != 0:
msg = 'Unable to create Plot ID={0} with mask background ' \
'{1} since 3 RGB values must be ' \
'input'.format(self._id, mask_background)
raise ValueError(msg)
check_type('plot mask background', mask_background, Iterable, Integral)
check_length('plot mask background', mask_background, 3)
for rgb in mask_background:
if not is_integer(rgb):
msg = 'Unable to create Plot ID={0} with mask background RGB ' \
'value {1} which is not an integer'.format(self._id, rgb)
raise ValueError(msg)
elif rgb < 0 or rgb > 255:
msg = 'Unable to create Plot ID={0} with mask bacground ' \
'RGB value {1} which is not between 0 and ' \
'255'.format(self._id, rgb)
raise ValueError(msg)
check_greater_than('plot mask background', rgb, 0, True)
check_less_than('plot mask background', rgb, 256)
self._mask_background = mask_background
def __repr__(self):
string = 'Plot\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
@ -409,8 +261,15 @@ class Plot(object):
string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec)
return string
def get_plot_xml(self):
"""Return XML representation of the plot
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing plot data
"""
element = ET.Element("plot")
element.set("id", str(self._id))
@ -422,64 +281,55 @@ class Plot(object):
element.set("basis", self._basis)
subelement = ET.SubElement(element, "origin")
text = ''
for coord in self._origin:
text += str(coord) + ' '
subelement.text = text.rstrip(' ')
subelement.text = ' '.join(map(str, self._origin))
subelement = ET.SubElement(element, "width")
text = ''
for dim in self._width:
text += str(dim) + ' '
subelement.text = text.rstrip(' ')
subelement.text = ' '.join(map(str, self._width))
subelement = ET.SubElement(element, "pixels")
text = ''
for dim in self._pixels:
text += str(dim) + ' '
subelement.text = text.rstrip(' ')
subelement.text = ' '.join(map(str, self._pixels))
if not self._mask_background is None:
if self._mask_background is not None:
subelement = ET.SubElement(element, "background")
text = ''
for rgb in self._background:
text += str(rgb) + ' '
subelement.text = text.rstrip(' ')
if not self._col_spec is None:
subelement.text = ' '.join(map(str, self._background))
if self._col_spec is not None:
for key in self._col_spec:
subelement = ET.SubElement(element, "col_spec")
subelement.set("id", '{0}'.format(key))
value = self._col_spec[key]
subelement.set("rgb",'{0} {1} {2}'.format(value[0],
value[1], value[2]))
subelement.set("id", str(key))
subelement.set("rgb", ' '.join(map(
str, self._col_spec[key])))
if not self._mask_components is None:
if self._mask_components is not None:
subelement = ET.SubElement(element, "mask")
text = ''
for id in self._mask_components:
text += str(id) + ' '
subelement.set("components", text.rstrip(' '))
rgb = self._mask_background
subelement.set("background", '{0} {1} {2}'.format(rgb[0],
rgb[1], rgb[2]))
subelement.set("components", ' '.join(map(
str, self._mask_components)))
subelement.set("background", ' '.join(map(
str, self._mask_background)))
return element
class PlotsFile(object):
"""Plots file used for an OpenMC simulation. Corresponds directly to the
plots.xml input file.
"""
def __init__(self):
# Initialize PlotsFile class attributes
self._plots = []
self._plots_file = ET.Element("plots")
def add_plot(self, plot):
"""Add a plot to the file.
Parameters
----------
plot : Plot
Plot to add
"""
if not isinstance(plot, Plot):
msg = 'Unable to add a non-Plot {0} to the PlotsFile'.format(plot)
@ -487,15 +337,20 @@ class PlotsFile(object):
self._plots.append(plot)
def remove_plot(self, plot):
"""Remove a plot from the file.
Parameters
----------
plot : Plot
Plot to remove
"""
self._plots.remove(plot)
def create_plot_subelements(self):
def _create_plot_subelements(self):
for plot in self._plots:
xml_element = plot.get_plot_xml()
if len(plot._name) > 0:
@ -503,10 +358,12 @@ class PlotsFile(object):
self._plots_file.append(xml_element)
def export_to_xml(self):
"""Create a plots.xml file that can be used by OpenMC.
self.create_plot_subelements()
"""
self._create_plot_subelements()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._plots_file)

1360
openmc/settings.py
View file

File diff suppressed because it is too large Load diff

180
openmc/statepoint.py
View file

@ -13,17 +13,28 @@ if sys.version > '3':
class SourceSite(object):
"""A single source site produced from fission.
Attributes
----------
weight : float
Weight of the particle arising from the site
xyz : list of float
Cartesian coordinates of the site
uvw : list of float
Directional cosines for particles emerging from the site
E : float
Energy of the emerging particle in MeV
"""
def __init__(self):
self._weight = None
self._xyz = None
self._uvw = None
self._E = None
def __repr__(self):
string = 'SourceSite\n'
string += '{0: <16}{1}{2}\n'.format('\tweight', '=\t', self._weight)
string += '{0: <16}{1}{2}\n'.format('\tE', '=\t', self._E)
@ -31,32 +42,56 @@ class SourceSite(object):
string += '{0: <16}{1}{2}\n'.format('\t(u,v,w)', '=\t', self._uvw)
return string
@property
def weight(self):
return self._weight
@property
def xyz(self):
return self._xyz
@property
def uvw(self):
return self._uvw
@property
def E(self):
return self._E
class StatePoint(object):
"""State information on a simulation at a certain point in time (at the end of a
given batch). Statepoints can be used to analyze tally results as well as
restart a simulation.
Attributes
----------
k_combined : list
Combined estimator for k-effective and its uncertainty
n_particles : int
Number of particles per generation
n_batches : int
Number of batches
current_batch :
Number of batches simulated
results : bool
Indicate whether tally results have been read
source : ndarray of SourceSite
Array of source sites
with_summary : bool
Indicate whether statepoint data has been linked against a summary file
tallies : dict
Dictionary whose keys are tally IDs and whose values are Tally objects
tallies_present : bool
Indicate whether user-defined tallies are present
global_tallies : ndarray
Global tallies and their uncertainties
n_realizations : int
Number of tally realizations
"""
def __init__(self, filename):
if filename.endswith('.h5'):
import h5py
self._f = h5py.File(filename, 'r')
@ -79,7 +114,6 @@ class StatePoint(object):
# Read tally metadata
self._read_tallies()
def close(self):
self._f.close()
@ -128,7 +162,6 @@ class StatePoint(object):
return self._n_realizations
def _read_metadata(self):
# Read filetype
self._filetype = self._get_int(path='filetype')[0]
@ -169,9 +202,7 @@ class StatePoint(object):
if self._run_mode == 2:
self._read_criticality()
def _read_criticality(self):
# Read criticality information
if self._run_mode == 2:
@ -191,9 +222,7 @@ class StatePoint(object):
# Read CMFD information (if used)
self._read_cmfd()
def _read_cmfd(self):
base = 'cmfd'
# Read CMFD information
@ -217,9 +246,7 @@ class StatePoint(object):
self._cmfd_srccmp = self._get_double(self._current_batch,
path='{0}/cmfd_srccmp'.format(base))
def _read_meshes(self):
# Initialize dictionaries for the Meshes
# Keys - Mesh IDs
# Values - Mesh objects
@ -282,9 +309,7 @@ class StatePoint(object):
# Add mesh to the global dictionary of all Meshes
self._meshes[mesh_id] = mesh
def _read_tallies(self):
# Initialize dictionaries for the Tallies
# Keys - Tally IDs
# Values - Tally objects
@ -423,7 +448,6 @@ class StatePoint(object):
# Add the scores to the Tally
for j, score in enumerate(scores):
# If this is a scattering moment, insert the scattering order
if '-n' in score:
score = score.replace('-n', '-' + str(moments[j]))
@ -437,10 +461,13 @@ class StatePoint(object):
# Add Tally to the global dictionary of all Tallies
self.tallies[tally_key] = tally
def read_results(self):
"""Read tally results and store them in the ``tallies`` attribute. No results
are read when the statepoint is instantiated.
# Number of realizations for global Tallies
"""
# Number of realizations for global Tallies
self._n_realizations = self._get_int(path='n_realizations')[0]
# Read global Tallies
@ -483,7 +510,8 @@ class StatePoint(object):
sum_sq = results[1::2]
# Define a routine to convert 0 to 1
nonzero = lambda val: 1 if not val else val
def nonzero(val):
return 1 if not val else val
# Reshape the results arrays
new_shape = (nonzero(tally.num_filter_bins),
@ -494,13 +522,17 @@ class StatePoint(object):
sum_sq = np.reshape(sum_sq, new_shape)
# Set the data for this Tally
tally.set_results(sum=sum, sum_sq=sum_sq)
tally.sum = sum
tally.sum_sq = sum_sq
# Indicate that Tally results have been read
self._results = True
def read_source(self):
"""Read and store source sites from the statepoint file. By default, source
sites are not loaded upon initialization.
"""
# Check whether Tally results have been read
if not self._results:
@ -520,7 +552,6 @@ class StatePoint(object):
# Initialize SourceSite object for each particle
for i in range(self._n_particles):
# Initialize new source site
site = SourceSite()
@ -536,9 +567,18 @@ class StatePoint(object):
# Store the source site in the NumPy array
self._source[i] = site
def compute_ci(self, confidence=0.95):
"""Computes confidence intervals for each Tally bin."""
"""Computes confidence intervals for each Tally bin.
This method is equivalent to calling compute_stdev(...) when the
confidence is known as opposed to its corresponding t value.
Parameters
----------
confidence : float, optional
Confidence level. Defaults to 0.95.
"""
# Determine significance level and percentile for two-sided CI
alpha = 1 - confidence
@ -548,11 +588,16 @@ class StatePoint(object):
t_value = scipy.stats.t.ppf(percentile, self._n_realizations - 1)
self.compute_stdev(t_value)
def compute_stdev(self, t_value=1.0):
"""
Computes the sample mean and the standard deviation of the mean
"""Computes the sample mean and the standard deviation of the mean
for each Tally bin.
Parameters
----------
t_value : float, optional
Student's t-value applied to the uncertainty. Defaults to 1.0,
meaning the reported value is the sample standard deviation.
"""
# Determine number of realizations
@ -572,49 +617,45 @@ class StatePoint(object):
if s != 0.0:
self._global_tallies[i, 1] = t_value * np.sqrt((s2 / n - s**2) / (n-1))
# Calculate sample mean and standard deviation for user-defined Tallies
for tally_id, tally in self.tallies.items():
tally.compute_std_dev(t_value)
def get_tally(self, scores=[], filters=[], nuclides=[],
name=None, id=None, estimator=None):
"""Finds and returns a Tally object with certain properties.
This routine searches the list of Tallies and returns the first Tally
found it finds which satisfieds all of the input parameters.
found it finds which satisfies all of the input parameters.
NOTE: The input parameters do not need to match the complete Tally
specification and may only represent a subset of the Tallies properties.
specification and may only represent a subset of the Tally's properties.
Parameters
----------
scores : list
A list of one or more score strings (default is []).
filters : list
A list of Filter objects (default is []).
nuclides : list
A list of Nuclide objects (default is []).
name : str
The name specified for the Tally (default is None).
id : int
The id specified for the Tally (default is None).
estimator: str
The type of estimator ('tracklength', 'analog'; default is None).
scores : list, optional
A list of one or more score strings (default is []).
filters : list, optional
A list of Filter objects (default is []).
nuclides : list, optional
A list of Nuclide objects (default is []).
name : str, optional
The name specified for the Tally (default is None).
id : int, optional
The id specified for the Tally (default is None).
estimator: str, optional
The type of estimator ('tracklength', 'analog'; default is None).
Returns
-------
A Tally object.
tally : Tally
A tally matching the specified criteria
Raises
------
LookupError : An error when a Tally meeting all of the input
parameters cannot be found in the statepoint.
LookupError
If a Tally meeting all of the input parameters cannot be found in
the statepoint.
"""
tally = None
@ -640,7 +681,7 @@ class StatePoint(object):
# Iterate over the scores requested by the user
for score in scores:
if not score in test_tally.scores:
if score not in test_tally.scores:
contains_scores = False
break
@ -653,7 +694,7 @@ class StatePoint(object):
# Iterate over the Filters requested by the user
for filter in filters:
if not filter in test_tally.filters:
if filter not in test_tally.filters:
contains_filters = False
break
@ -666,7 +707,7 @@ class StatePoint(object):
# Iterate over the Nuclides requested by the user
for nuclide in nuclides:
if not nuclide in test_tally.nuclides:
if nuclide not in test_tally.nuclides:
contains_nuclides = False
break
@ -683,14 +724,15 @@ class StatePoint(object):
return tally
def link_with_summary(self, summary):
"""Links Tallies and Filters with Summary model information.
This routine retrieves model information (materials, geometry) from a
Summary object populated with an HDF5 'summary.h5' file and inserts
it into the Tally objects. This can be helpful when viewing and
manipulating large scale Tally data.
Summary object populated with an HDF5 'summary.h5' file and inserts it
into the Tally objects. This can be helpful when viewing and
manipulating large scale Tally data. Note that it is necessary to link
against a summary to populate the Tallies with any user-specified "name"
XML tags.
Parameters
----------
@ -699,8 +741,10 @@ class StatePoint(object):
Raises
------
ValueError : An error when the argument passed to the 'summary'
parameter is not an openmc.Summary object.
ValueError
An error when the argument passed to the 'summary' parameter is not
an openmc.Summary object.
"""
if not isinstance(summary, openmc.summary.Summary):
@ -709,7 +753,6 @@ class StatePoint(object):
raise ValueError(msg)
for tally_id, tally in self.tallies.items():
# Get the Tally name from the summary file
tally.name = summary.tallies[tally_id].name
tally.with_summary = True
@ -717,7 +760,6 @@ class StatePoint(object):
nuclide_zaids = copy.deepcopy(tally.nuclides)
for nuclide_zaid in nuclide_zaids:
tally.remove_nuclide(nuclide_zaid)
if nuclide_zaid == -1:
tally.add_nuclide(openmc.Nuclide('total'))
@ -725,7 +767,6 @@ class StatePoint(object):
tally.add_nuclide(summary.nuclides[nuclide_zaid])
for filter in tally.filters:
if filter.type == 'surface':
surface_ids = []
for bin in filter.bins:
@ -749,12 +790,11 @@ class StatePoint(object):
for bin in filter.bins:
material_ids.append(summary.materials[bin].id)
filter.bins = material_ids
self._with_summary = True
def _get_data(self, n, typeCode, size):
return list(struct.unpack('={0}{1}'.format(n,typeCode),
return list(struct.unpack('={0}{1}'.format(n, typeCode),
self._f.read(n*size)))
def _get_int(self, n=1, path=None):

155
openmc/summary.py
View file

@ -2,16 +2,18 @@ import numpy as np
import openmc
try:
import h5py
except ImportError:
msg = 'Unable to import h5py which is needed by openmc.summary'
raise ImportError(msg)
class Summary(object):
"""Information summarizing the geometry, materials, and tallies used in a
simulation.
"""
def __init__(self, filename):
# A user may not have h5py, but they can still use the rest of the
# Python API so we'll only try to import h5py if the user actually inits
# a Summary object.
import h5py
openmc.reset_auto_ids()
@ -27,9 +29,7 @@ class Summary(object):
self._read_geometry()
self._read_tallies()
def _read_metadata(self):
# Read OpenMC version
self.version = [self._f['version_major'][0],
self._f['version_minor'][0],
@ -44,9 +44,7 @@ class Summary(object):
self.gen_per_batch = self._f['gen_per_batch'][0]
self.n_procs = self._f['n_procs'][0]
def _read_geometry(self):
# Read in and initialize the Materials and Geometry
self._read_nuclides()
self._read_materials()
@ -56,9 +54,7 @@ class Summary(object):
self._read_lattices()
self._finalize_geometry()
def _read_nuclides(self):
self.n_nuclides = self._f['nuclides/n_nuclides'][0]
# Initialize dictionary for each Nuclide
@ -67,7 +63,6 @@ class Summary(object):
self.nuclides = {}
for key in self._f['nuclides'].keys():
if key == 'n_nuclides':
continue
@ -88,9 +83,7 @@ class Summary(object):
self.nuclides[zaid] = openmc.Nuclide(name=name, xs=xs)
self.nuclides[zaid].zaid = zaid
def _read_materials(self):
self.n_materials = self._f['materials/n_materials'][0]
# Initialize dictionary for each Material
@ -99,7 +92,6 @@ class Summary(object):
self.materials = {}
for key in self._f['materials'].keys():
if key == 'n_materials':
continue
@ -116,7 +108,6 @@ class Summary(object):
# Read the names of the S(a,b) tables for this Material
for i in range(1, n_sab+1):
sab_table = self._f['materials'][key]['sab_tables'][str(i)][0]
# Read the cross-section identifiers for each S(a,b) table
@ -148,9 +139,7 @@ class Summary(object):
# Add the Material to the global dictionary of all Materials
self.materials[index] = material
def _read_surfaces(self):
self.n_surfaces = self._f['geometry/n_surfaces'][0]
# Initialize dictionary for each Surface
@ -159,7 +148,6 @@ class Summary(object):
self.surfaces = {}
for key in self._f['geometry/surfaces'].keys():
if key == 'n_surfaces':
continue
@ -171,7 +159,6 @@ class Summary(object):
coeffs = self._f['geometry/surfaces'][key]['coefficients'][...]
# Create the Surface based on its type
if surf_type == 'X Plane':
x0 = coeffs[0]
surface = openmc.XPlane(surface_id, bc, x0, name)
@ -232,9 +219,7 @@ class Summary(object):
# Add Surface to global dictionary of all Surfaces
self.surfaces[index] = surface
def _read_cells(self):
self.n_cells = self._f['geometry/n_cells'][0]
# Initialize dictionary for each Cell
@ -251,7 +236,6 @@ class Summary(object):
self._cell_fills = {}
for key in self._f['geometry/cells'].keys():
if key == 'n_cells':
continue
@ -310,9 +294,7 @@ class Summary(object):
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell
def _read_universes(self):
self.n_universes = self._f['geometry/n_universes'][0]
# Initialize dictionary for each Universe
@ -321,7 +303,6 @@ class Summary(object):
self.universes = {}
for key in self._f['geometry/universes'].keys():
if key == 'n_universes':
continue
@ -340,9 +321,7 @@ class Summary(object):
# Add the Universe to the global list of Universes
self.universes[index] = universe
def _read_lattices(self):
self.n_lattices = self._f['geometry/n_lattices'][0]
# Initialize lattices for each Lattice
@ -351,7 +330,6 @@ class Summary(object):
self.lattices = {}
for key in self._f['geometry/lattices'].keys():
if key == 'n_lattices':
continue
@ -383,7 +361,6 @@ class Summary(object):
lattice.lower_left = lower_left
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void (-22)
if outer != -22:
lattice.outer = self.universes[outer]
@ -395,14 +372,14 @@ class Summary(object):
for x in range(universe_ids.shape[0]):
for y in range(universe_ids.shape[1]):
for z in range(universe_ids.shape[2]):
universes[x,y,z] = \
self.get_universe_by_id(universe_ids[x,y,z])
universes[x, y, z] = \
self.get_universe_by_id(universe_ids[x, y, z])
# Transpose, reverse y-dimension for appropriate ordering
shape = universes.shape
universes = np.transpose(universes, (1,0,2))
universes = np.transpose(universes, (1, 0, 2))
universes.shape = shape
universes = universes[:,::-1,:]
universes = universes[:, ::-1, :]
lattice.universes = universes
if offset_size > 0:
@ -420,8 +397,8 @@ class Summary(object):
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'][0]
universe_ids = \
self._f['geometry/lattices'][key]['universes'][...]
universe_ids = self._f[
'geometry/lattices'][key]['universes'][...]
# Create the Lattice
lattice = openmc.HexLattice(lattice_id=lattice_id, name=name)
@ -441,9 +418,9 @@ class Summary(object):
for i in range(universe_ids.shape[0]):
for j in range(universe_ids.shape[1]):
for k in range(universe_ids.shape[2]):
if universe_ids[i,j,k] != -1:
universes[i,j,k] = \
self.get_universe_by_id(universe_ids[i,j,k])
if universe_ids[i, j, k] != -1:
universes[i, j, k] = self.get_universe_by_id(
universe_ids[i, j, k])
if offset_size > 0:
lattice.offsets = offsets
@ -451,15 +428,12 @@ class Summary(object):
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
def _finalize_geometry(self):
# Initialize Geometry object
self.openmc_geometry = openmc.Geometry()
# Iterate over all Cells and add fill Materials, Universes and Lattices
for cell_key in self._cell_fills.keys():
# Determine fill type ('normal', 'universe', or 'lattice') and ID
fill_type = self._cell_fills[cell_key][0]
fill_id = self._cell_fills[cell_key][1]
@ -482,16 +456,14 @@ class Summary(object):
root_universe = self.get_universe_by_id(0)
self.openmc_geometry.root_universe = root_universe
def _read_tallies(self):
# Initialize dictionaries for the Tallies
# Keys - Tally IDs
# Values - Tally objects
self.tallies = {}
# Read the number of tallies
if not 'tallies' in self._f.keys():
if 'tallies' not in self._f:
self.n_tallies = 0
return
@ -554,8 +526,11 @@ class Summary(object):
# Add Tally to the global dictionary of all Tallies
self.tallies[tally_id] = tally
def make_opencg_geometry(self):
"""Create OpenCG geometry based on the information contained in the summary
file. The geometry is stored as the 'opencg_geometry' attribute.
"""
try:
from openmc.opencg_compatible import get_opencg_geometry
@ -567,8 +542,22 @@ class Summary(object):
if self.opencg_geometry is None:
self.opencg_geometry = get_opencg_geometry(self.openmc_geometry)
def get_nuclide_by_zaid(self, zaid):
"""Return a Nuclide object given the 'zaid' identifier for the nuclide.
Parameters
----------
zaid : int
1000*Z + A, where Z is the atomic number of the nuclide and A is the
mass number. For example, the zaid for U-235 is 92235.
Returns
-------
nuclide : openmc.nuclide.Nuclide or None
Nuclide matching the specified zaid, or None if no matching object
is found.
"""
for index, nuclide in self.nuclides.items():
if nuclide._zaid == zaid:
@ -576,8 +565,20 @@ class Summary(object):
return None
def get_material_by_id(self, material_id):
"""Return a Material object given the material id
Parameters
----------
id : int
Unique identifier for the material
Returns
-------
material : openmc.material.Material
Material with given id
"""
for index, material in self.materials.items():
if material._id == material_id:
@ -585,8 +586,20 @@ class Summary(object):
return None
def get_surface_by_id(self, surface_id):
"""Return a Surface object given the surface id
Parameters
----------
id : int
Unique identifier for the surface
Returns
-------
surface : openmc.surface.Surface
Surface with given id
"""
for index, surface in self.surfaces.items():
if surface._id == surface_id:
@ -594,8 +607,20 @@ class Summary(object):
return None
def get_cell_by_id(self, cell_id):
"""Return a Cell object given the cell id
Parameters
----------
id : int
Unique identifier for the cell
Returns
-------
cell : openmc.universe.Cell
Cell with given id
"""
for index, cell in self.cells.items():
if cell._id == cell_id:
@ -603,8 +628,20 @@ class Summary(object):
return None
def get_universe_by_id(self, universe_id):
"""Return a Universe object given the universe id
Parameters
----------
id : int
Unique identifier for the universe
Returns
-------
universe : openmc.universe.Universe
Universe with given id
"""
for index, universe in self.universes.items():
if universe._id == universe_id:
@ -612,8 +649,20 @@ class Summary(object):
return None
def get_lattice_by_id(self, lattice_id):
"""Return a Lattice object given the lattice id
Parameters
----------
id : int
Unique identifier for the lattice
Returns
-------
lattice : openmc.universe.Lattice
Lattice with given id
"""
for index, lattice in self.lattices.items():
if lattice._id == lattice_id:

767
openmc/surface.py
View file

File diff suppressed because it is too large Load diff

596
openmc/tallies.py
View file

File diff suppressed because it is too large Load diff

84
openmc/trigger.py
View file

@ -1,25 +1,47 @@
from numbers import Real
from xml.etree import ElementTree as ET
import sys
from openmc.checkvalue import *
from openmc.checkvalue import check_type, check_value
if sys.version_info[0] >= 3:
basestring = str
class Trigger(object):
"""A criterion for when to finish a simulation based on tally uncertainties.
Parameters
----------
trigger_type : {'variance', 'std_dev', 'rel_err'}
Determine whether to trigger on the variance, standard deviation, or
relative error of scores.
threshold : float
The threshold for the trigger type.
Attributes
----------
trigger_type : {'variance', 'std_dev', 'rel_err'}
Determine whether to trigger on the variance, standard deviation, or
relative error of scores.
threshold : float
The threshold for the trigger type.
scores : list of str
Scores which should be checked against the trigger
"""
def __init__(self, trigger_type, threshold):
# Initialize Mesh class attributes
self.trigger_type = trigger_type
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
@ -36,47 +58,40 @@ class Trigger(object):
else:
return existing
@property
def trigger_type(self):
return self._trigger_type
@property
def threshold(self):
return self._threshold
@property
def scores(self):
return self._scores
@trigger_type.setter
def trigger_type(self, trigger_type):
if not trigger_type in ['variance', 'std_dev', 'rel_err']:
msg = 'Unable to create a tally trigger with ' \
'type "{0}"'.format(trigger_type)
raise ValueError(msg)
check_value('tally trigger type', trigger_type,
['variance', 'std_dev', 'rel_err'])
self._trigger_type = trigger_type
@threshold.setter
def threshold(self, threshold):
if not is_float(threshold):
msg = 'Unable to set a tally trigger threshold with ' \
'threshold "{0}"'.format(threshold)
raise ValueError(msg)
check_type('tally trigger threshold', threshold, Real)
self._threshold = threshold
def add_score(self, score):
"""Add a score to the list of scores to be checked against the trigger.
if not is_string(score):
Parameters
----------
score : str
Score to append
"""
if not isinstance(score, basestring):
msg = 'Unable to add score "{0}" to tally trigger since ' \
'it is not a string'.format(score)
raise ValueError(msg)
@ -87,28 +102,25 @@ class Trigger(object):
else:
self._scores.append(score)
def __repr__(self):
string = 'Trigger\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._trigger_type)
string += '{0: <16}{1}{2}\n'.format('\tThreshold', '=\t', self._threshold)
string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self._scores)
return string
def get_trigger_xml(self, element):
"""Return XML representation of the trigger
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing trigger data
"""
subelement = ET.SubElement(element, "trigger")
subelement.set("type", self._trigger_type)
subelement.set("threshold", str(self._threshold))
# Scores
if len(self._scores) != 0:
scores = ''
for score in self._scores:
scores += '{0} '.format(score)
scores.rstrip(' ')
subelement.set("scores", scores)
subelement.set("scores", ' '.join(map(str, self._scores)))

889
openmc/universe.py
View file

File diff suppressed because it is too large Load diff