Added python api to build mgxs library, functionality should be there, but i havent done enough testing to be confident yet.

This commit is contained in:
Adam Nelson 2015-11-15 15:03:23 -05:00
parent d125bf5bd5
commit 683f782744
3 changed files with 681 additions and 4 deletions

View file

@ -177,7 +177,7 @@ attributes/sub-elements required to describe the meta-data:
*Default*: ``33``
The following attributes/sub-elements are the actual cross section values to
The following attributes/sub-elements are the cross section values to
be used during the transport process.
:total:

View file

@ -503,10 +503,10 @@ class Material(object):
def _get_macroscopic_xml(self, macroscopic, distrib=False):
xml_element = ET.Element("macroscopic")
xml_element.set("name", macroscopic[0]._name)
xml_element.set("name", macroscopic._name)
if macroscopic[0].xs is not None:
xml_element.set("xs", macroscopic[0].xs)
if macroscopic.xs is not None:
xml_element.set("xs", macroscopic.xs)
return xml_element

677
openmc/mgxs_library.py Normal file
View file

@ -0,0 +1,677 @@
from collections import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import warnings
import sys
if sys.version_info[0] >= 3:
basestring = str
import numpy as np
import openmc
from openmc.mgxs import EnergyGroups
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.clean_xml import *
# MGXS Representations supported by OpenMC
REPRESENTATIONS = ['isotropic', 'angle']
def ndarray_to_string(arr):
"""Converts a numpy ndarray in to a join with spaces between entries
similar to ' '.join(map(str,arr)) but applied to all sub-dimensions.
"""
shape = arr.shape
ndim = arr.ndim
text = ''
if ndim == 1:
text += ' '.join(map(str, arr[:]))
elif ndim == 2:
for i in xrange(shape[0]):
text += ' '.join(map(str, arr[i,:]))
text += '\n'
elif ndim == 3:
for i in xrange(shape[0]):
for j in xrange(shape[1]):
text += ' '.join(map(str, arr[i,j,:]))
text += '\n'
elif ndim == 4:
for i in xrange(shape[0]):
for j in xrange(shape[1]):
for k in xrange(shape[2]):
text += ' '.join(map(str, arr[i,j,k,:]))
text += '\n'
elif ndim == 5:
for i in xrange(shape[0]):
for j in xrange(shape[1]):
for k in xrange(shape[2]):
for l in xrange(shape[3]):
text += ' '.join(map(str, arr[i,j,k,l,:]))
text += '\n'
return text
class Xsdata(object):
"""A multi-group cross section data set (xsdata) providing all the
multi-group data necessary for a multi-group OpenMC calculation.
Parameters
----------
name : str, optional
Name of the mgxs data set.
representation : str
Method used in generating the MGXS (isotropic or angle-dependent flux
weighting). Defaults to 'isotropic'
Attributes
----------
name : str
Unique identifier for the xsdata object
alias : str
Separate unique identifier for the xsdata object
kT : float
Temperature (in units of MeV) of this data set.
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure
fissionable : boolean
Whether or not this is a fissionable data set.
scatt_type : str
Angular distribution representation (legendre, histogram, or tabular)
order : int
Either the Legendre order, number of bins, or number of points used to
describe the angular distribution associated with each group-to-group
transfer probability.
tabular_legendre : dict
Set how to treat the Legendre scattering kernel (tabular or leave in
Legendre polynomial form). Dict contains two keys: ``enable`` and
``num_points``. ``enable`` is a boolean and ``num_points`` is the
number of points to use, if ``enable`` is True.
"""
def __init__(self, name, energy_groups, representation="isotropic"):
# Initialize class attributes
self._name = name
self._energy_groups = energy_groups
self._representation = representation
self._alias = None
self._kT = None
self._fissionable = False
self._scatt_type = 'legendre'
self._order = None
self._tabular_legendre = None
self._num_polar = None
self._num_azimuthal = None
self._total = None
self._absorption = None
self._scatter = None
self._multiplicity = None
self._fission = None
self._nu_fission = None
self._k_fission = None
self._chi = None
@property
def name(self):
return self._name
@property
def energy_groups(self):
return self._energy_groups
@property
def alias(self):
return self._alias
@property
def kT(self):
return self._kT
@property
def scatt_type(self):
return self._scatt_type
@property
def order(self):
return self._order
@property
def tabular_legendre(self):
return self._tabular_legendre
@property
def num_polar(self):
return self._num_polar
@property
def num_azimuthal(self):
return self._num_azimuthal
@property
def total(self):
return self._total
@property
def absorption(self):
return self._absorption
@property
def scatter(self):
return self._scatter
@property
def multiplicity(self):
return self._multiplicity
@property
def fission(self):
return self._fission
@property
def nu_fission(self):
return self._nu_fission
@property
def k_fission(self):
return self._k_fission
@property
def chi(self):
return self._chi
@property
def num_orders(self):
if (self._order is not None) and (self._scatt_type is not None):
if self._scatt_type is 'legendre':
return self._order + 1
else:
return self._order
@name.setter
def name(self, name):
check_type('name for Xsdata', name, basestring)
self._name = name
@energy_groups.setter
def energy_groups(self, energy_groups):
# Check validity of energy_groups
check_type("energy_groups", energy_groups, EnergyGroups)
# Check that there is one or more groups
if (energy_groups.num_energy_groups.num_group is None) or (energy_groups.num_energy_groups.num_group < 1):
msg = 'energy_groups object incorrectly initialized.'
raise ValueError(msg)
self._energy_groups = energy_groups
@representation.setter
def representation(self, representation):
# Check it is of valid type.
check_value('representation', representation, REPRESENTATIONS)
self._representation = representation
@alias.setter
def alias(self, alias):
if alias is not None:
check_type('alias for Xsdata', alias, basestring)
self._alias = alias
else:
self._alias = self._name
@kT.setter
def kT(self, kT):
# Check validity of type and that the kT value is >= 0
check_type("kT", kT, Real)
check_greater_than("kT", kT, 0.0, equality=True)
self._kT = kT
@scatt_type.setter
def scatt_type(self, scatt_type):
# check to see it is of a valid type and value
check_value("scatt_type", scatt_type, ['legendre', 'histogram',
'tabular'])
self._scatt_type = scatt_type
@order.setter
def order(self, order):
# Check type and value
check_type("order", order, Integral)
check_greater_than("order", order, 0, equality=True)
self._order = order
@tabular_legendre.setter
def tabular_legendre(self, tabular_legendre):
# Check to make sure this is a dict and it has our keys with the
# right values.
check_type("tabular_legendre", tabular_legendre, dict)
if 'enable' in tabular_legendre:
enable = tabular_legendre['enable']
check_type('enable', enable, bool)
else:
msg = "enable must be provided in tabular_legendre"
raise ValueError(msg)
if 'num_points' in tabular_legendre:
num_points = tabular_legendre['num_points']
check_value('num_points', num_points, Integral)
check_greater_than('num_points', num_points, 0)
else:
num_points = 33
self._tabular_legendre = {'enable': enable, 'num_points': num_points}
@num_polar.setter(self, num_polar):
# Make sure we have positive ints
check_value("num_polar", num_polar, Integral)
check_greater_than("num_polar", num_polar, 0)
self._num_polar = num_polar
@num_azimuthal.setter(self, num_azimuthal):
check_value("num_azimuthal", num_azimuthal, Integral)
check_greater_than("num_azimuthal", num_azimuthal, 0)
self._num_azimuthal = num_azimuthal
@total.setter
def total(self, total):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
# check we have a numpy list
check_type("total", total, np.ndarray, expected_iter_type=Real)
if total.shape == shape:
self._total = np.copy(total)
else:
msg = 'Shape of provided total "{0}" does not match shape ' \
'required, "{1}"'.format(total.shape, shape)
raise ValueError(msg)
@absorption.setter
def absorption(self, absorption):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
# check we have a numpy list
check_type("absorption", absorption, np.ndarray, expected_iter_type=Real)
if absorption.shape == shape:
self._absorption = np.copy(absorption)
else:
msg = 'Shape of provided absorption "{0}" does not match shape ' \
'required, "{1}"'.format(absorption.shape, shape)
raise ValueError(msg)
@fission.setter
def fission(self, fission):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
# check we have a numpy list
check_type("fission", fission, np.ndarray, expected_iter_type=Real)
if fission.shape == shape:
self._fission = np.copy(fission)
if np.sum(self._fission) > 0.0:
self._fissionable = True
else:
msg = 'Shape of provided fission "{0}" does not match shape ' \
'required, "{1}"'.format(fission.shape, shape)
raise ValueError(msg)
@k_fission.setter
def k_fission(self, k_fission):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
# check we have a numpy list
check_type("k_fission", k_fission, np.ndarray, expected_iter_type=Real)
if k_fission.shape == shape:
self._k_fission = np.copy(k_fission)
if np.sum(self._k_fission) > 0.0:
self._fissionable = True
else:
msg = 'Shape of provided k_fission "{0}" does not match shape ' \
'required, "{1}"'.format(k_fission.shape, shape)
raise ValueError(msg)
@chi.setter
def chi(self, chi):
if self._use_chi is not None:
msg = 'Providing chi when nu_fission already provided as matrix!'
raise ValueError(msg)
if self._representation is 'isotropic':
shape = (self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
# check we have a numpy list
check_type("chi", chi, np.ndarray, expected_iter_type=Real)
if chi.shape == shape:
self._chi = np.copy(chi)
else:
msg = 'Shape of provided chi "{0}" does not match shape ' \
'required, "{1}"'.format(chi.shape, shape)
raise ValueError(msg)
if self._use_chi is not None:
self._use_chi = True
@scatter.setter
def scatter(self, scatter):
if self._representation is 'isotropic':
shape = (self.num_orders, self._energy_groups.num_group,
self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal, self.num_orders,
self._energy_groups.num_group,
self._energy_groups.num_group)
# check we have a numpy list
check_type("scatter", scatter, np.ndarray, expected_iter_type=Real)
if scatter.shape == shape:
self._scatter = np.copy(scatter)
else:
msg = 'Shape of provided scatter "{0}" does not match shape ' \
'required, "{1}"'.format(scatter.shape, shape)
raise ValueError(msg)
@multiplicity.setter
def multiplicity(self, multiplicity):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_group,
self._energy_groups.num_group)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group,
self._energy_groups.num_group)
# check we have a numpy list
check_type("multiplicity", multiplicity, np.ndarray,
expected_iter_type=Real)
if multiplicity.shape == shape:
self._multiplicity = np.copy(multiplicity)
else:
msg = 'Shape of provided multiplicity "{0}" does not match shape ' \
'required, "{1}"'.format(multiplicity.shape, shape)
raise ValueError(msg)
@nu_fission.setter
def nu_fission(self, nu_fission):
# nu_fission ca nbe given as a vector or a matrix
# Vector is used when chi also exists.
# Matrix is used when chi does not exist.
# We have to check that the correct form is given, but only if
# chi already has been set. If not, we just check that this is OK
# and set the use_chi flag.
# First lets set our dimensions here since they get used repeatedly
# throughout this code.
if self._representation is 'isotropic':
shape_vec = (self._energy_groups.num_group)
shape_mat = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
elif self._representation is 'angle':
shape_vec = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group)
shape_mat = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_group,
self._energy_groups.num_group)
# Begin by checking the case when chi has already been given and thus
# the rules for filling in nu_fission are set.
if self._use_chi is not None:
if self._use_chi:
shape = shape_vec
else:
shape = shape_mat
if nu_fission.shape /= shape:
msg = "Invalid Shape of Nu_fission!"
raise ValueError(msg)
else:
# Get shape of nu_fission so we can figure if we need chi or not
if nu_fission.shape == shape_vec:
self._use_chi = True
shape = shape_vec
elif nu_fission.shape = shape_mat:
self._use_chi = False
shape = shape_mat
else:
msg = "Invalid Shape of Nu_fission!"
raise ValueError(msg)
# check we have a numpy list
check_type("nu_fission", nu_fission, np.ndarray, expected_iter_type=Real)
self._nu_fission = np.copy(nu_fission)
def _get_xsdata_xml(self):
element = ET.Element("xsdata")
element.set("name", xsdata._name)
if xsdata._alias is not None:
subelement = ET.SubElement(element, 'alias')
subelement.text(xsdata.alias)
if xsdata._kT is not None:
subelement = ET.SubElement(element, 'kT')
subelement.text(str(self._kT))
if xsdata._fissionable is not None:
subelement = ET.SubElement(element, 'fissionable')
subelement.text(str(self._fissionable))
if xsdata._representation is not None:
subelement = ET.SubElement(element, 'representation')
subelement.text(self._representation)
if xsdata._representation == 'angle':
if xsdata._num_azimuthal is not None:
subelement = ET.SubElement(element, 'num_azimuthal')
subelement.text(str(self._num_azimuthal))
if xsdata._num_polar is not None:
subelement = ET.SubElement(element, 'num_polar')
subelement.text(str(self._num_polar))
if xsdata._scatt_type is not None:
subelement = ET.SubElement(element, 'scatt_type')
subelement.text(self._scatt_type)
if xsdata._order is not None:
subelement = ET.SubElement(element, 'order')
subelement.text(str(self._order))
if xsdata._tabular_legendre is not None:
subelement = ET.SubElement(element, 'tabular_legendre')
subelement.set('enable', str(xsdata._tabular_legendre['enable']))
subelement.set('num_points', str(xsdata._tabular_legendre['num_points']))
if self._total is not None:
subelement = ET.SubElement(element, 'total')
subelement.text(ndarray_to_string(self._total))
if self._absorption is not None:
subelement = ET.SubElement(element, 'absorption')
subelement.text(ndarray_to_string(self._absorption))
if self._scatter is not None:
subelement = ET.SubElement(element, 'scatter')
subelement.text(ndarray_to_string(self._scatter))
if self._multiplicity is not None:
subelement = ET.SubElement(element, 'multiplicity')
subelement.text(ndarray_to_string(self._multiplicity))
if self._fissionable:
if self._fission is not None:
subelement = ET.SubElement(element, 'fission')
subelement.text(ndarray_to_string(self._fission))
if self._k_fission is not None:
subelement = ET.SubElement(element, 'k_fission')
subelement.text(ndarray_to_string(self._k_fission))
if self._nu_fission is not None:
subelement = ET.SubElement(element, 'nu_fission')
subelement.text(ndarray_to_string(self._nu_fission))
if self._chi is not None:
subelement = ET.SubElement(element, 'chi')
subelement.text(ndarray_to_string(self._chi))
return element
class MGXSLibraryFile(object):
"""Multi-Group Cross Sections file used for an OpenMC simulation.
Corresponds directly to the MG version of the cross_sections.xml input file.
Attributes
----------
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure.
inverse_velocities : Iterable of Real
Inverse of velocities, units of sec/cm
filename : str
XML file to write to.
"""
def __init__(self, energy_groups):
# Initialize MGXSLibraryFile class attributes
self._xsdatas = []
self._energy_groups = energy_groups
self._inverse_velocities = None
self._cross_sections_file = ET.Element("cross_sections")
@property
def inverse_velocities(self):
return self._inverse_velocities
@property
def energy_groups(self):
return self._energy_groups
@inverse_velocities.setter
def inverse_velocities(self, inverse_velocities):
cv.check_type('inverse_velocities', inverse_velocities, Iterable, Real)
cv.check_greater_than('number of inverse_velocities',
len(inverse_velocities), 0.0)
self._inverse_velocities = np.array(inverse_velocities)
@energy_groups.setter
def energy_groups(self, energy_groups):
check_type("energy groups", energy_groups, EnergyGroups)
self._energy_groups = energy_groups
def add_xsdata(self, xsdata):
"""Add an xsdata entry to the file.
Parameters
----------
xsdata : Xsdata
MGXS information to add
"""
# Check the type
if not isinstance(xsdata, Xsdata):
msg = 'Unable to add a non-Xsdata "{0}" to the ' \
'MGXSLibraryFile'.format(xsdata)
raise ValueError(msg)
# Make sure energy groups match.
if xsdata.energy_groups /= self._energy_groups:
msg = 'Energy groups of Xsdata do not match that of MGXSLibraryFile!'
raise ValueError(msg)
self._xsdatas.append(xsdata)
def add_xsdatas(self, xsdatas):
"""Add multiple xsdatas to the file.
Parameters
----------
xsdatas : tuple or list of Xsdata
Xsdatas to add
"""
if not isinstance(xsdatas, Iterable):
msg = 'Unable to create OpenMC xsdatas.xml file from "{0}" which ' \
'is not iterable'.format(xsdatas)
raise ValueError(msg)
for xsdata in xsdatas:
self.add_xsdata(xsdata)
def remove_xsdata(self, xsdata):
"""Remove a xsdata from the file
Parameters
----------
xsdata : Xsdata
Xsdata to remove
"""
if not isinstance(xsdata, Xsdata):
msg = 'Unable to remove a non-Xsdata "{0}" from the ' \
'XsdatasFile'.format(xsdata)
raise ValueError(msg)
self._xsdatas.remove(xsdata)
def _create_groups_subelement(self):
if self._energy_groups is not None:
element = ET.SubElement(self._cross_sections_file, "groups")
element.text = str(self._energy_groups.num_group)
def _create_group_structure_subelement(self):
if self._energy_groups is not None:
element = ET.SubElement(self._cross_sections_file,
"group_structure")
element.text = ' '.join(map(str, self._energy_groups.group_edges))
def _create_inverse_velocities_subelement(self):
if self._inverse_velocities is not None:
element = ET.SubElement(self._cross_sections_file,
"inverse_velocities")
element.text = ' '.join(map(str, self._inverse_velocities))
def _create_xsdata_subelements(self):
for xsdata in self._xsdatas:
xml_element = xsdata.get_xsdata_xml()
self._cross_sections_file.append(xml_element)
def export_to_xml(self, filename='mg_cross_sections.xml'):
"""Create an mg_cross_sections.xml file that can be used for a
simulation.
Parameters
----------
filename : str, optional
filename of file, default is mg_cross_sections.xml
"""
# Reset xml element tree
self._cross_sections_file.clear()
self._create_groups_subelement()
self._create_group_structure_subelement()
self._create_inverse_velocities_subelement()
self._create_xsdata_subelements()
# Clean the indentation in the file to be user-readable
sort_xml_elements(self._cross_sections_file)
clean_xml_indentation(self._cross_sections_file)
# Write the XML Tree to the xsdatas.xml file
tree = ET.ElementTree(self._cross_sections_file)
tree.write(filename, xml_declaration=True,
encoding='utf-8', method="xml")