diff --git a/src/utils/openmc/__init__.py b/src/utils/openmc/__init__.py new file mode 100644 index 0000000000..7902f18e43 --- /dev/null +++ b/src/utils/openmc/__init__.py @@ -0,0 +1,14 @@ +#!/usr/bin/env python + +from openmc.element import * +from openmc.geometry import * +from openmc.nuclide import * +from openmc.material import * +from openmc.opencsg_compatible import * +from openmc.plots import * +from openmc.settings import * +from openmc.surface import * +from openmc.universe import * +from openmc.tallies import * +from openmc.executor import * +#from statepoint import * diff --git a/src/utils/openmc/checkvalue.py b/src/utils/openmc/checkvalue.py new file mode 100644 index 0000000000..b99ed2cb48 --- /dev/null +++ b/src/utils/openmc/checkvalue.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python + +import numpy as np + + +def is_integer(val): + return isinstance(val, (int, np.int32, np.int64)) + + +def is_float(val): + return isinstance(val, (float, np.float32, np.float64)) + + +def is_string(val): + return isinstance(val, (str, np.str)) + diff --git a/src/utils/openmc/clean_xml.py b/src/utils/openmc/clean_xml.py new file mode 100644 index 0000000000..b63b62ea2f --- /dev/null +++ b/src/utils/openmc/clean_xml.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python + +def sort_xml_elements(tree): + + # Retrieve all children of the root XML node in the tree + elements = tree.getchildren() + + # Initialize empty lists for the sorted and comment elements + sorted_elements = list() + + # Initialize an empty set of tags (e.g., Surface, Cell, and Lattice) + tags = set() + + # Find the unique tags in the tree + for element in elements: + tags.add(element.tag) + + # Initialize an empty list for the comment elements + comment_elements = list() + + # Find the comment elements and record their ordering within the + # tree using a precedence with respect to the subsequent nodes + for index, element in enumerate(elements): + next_element = None + + if 'Comment' in str(element.tag): + + if index < len(elements)-1: + next_element = elements[index+1] + + comment_elements.append((element, next_element)) + + # Now iterate over all tags and order the elements within each tag + for tag in tags: + + # Retrieve all of the elements for this tag + try: + tagged_elements = tree.findall(tag) + except: + continue + + # Initialize an empty list of tuples to sort (id, element) + tagged_data = list() + + # Retrieve the IDs for each of the elements + for element in tagged_elements: + key = element.get('id') + + # If this element has an "ID" tag, append it to the list to sort + if not key is None: + tagged_data.append((key, element)) + + # Sort the elements according to the IDs for this tag + tagged_data.sort() + sorted_elements.extend(list(item[-1] for item in tagged_data)) + + # Add the comment elements while preserving the original precedence + for element, next_element in comment_elements: + index = sorted_elements.index(next_element) + sorted_elements.insert(index, element) + + # Remove all of the sorted elements from the tree + for element in sorted_elements: + tree.remove(element) + + # Add the sorted elements back to the tree in the proper order + tree.extend(sorted_elements) + + +def clean_xml_indentation(element, level=0): + ''' + copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint + it basically walks your tree and adds spaces and newlines so the tree is + printed in a nice way + ''' + + i = "\n" + level*" " + + if len(element): + + if not element.text or not element.text.strip(): + element.text = i + " " + + if not element.tail or not element.tail.strip(): + element.tail = i + + for element in element: + clean_xml_indentation(element, level+1) + + if not element.tail or not element.tail.strip(): + element.tail = i + + else: + if level and (not element.tail or not element.tail.strip()): + element.tail = i \ No newline at end of file diff --git a/src/utils/openmc/constants.py b/src/utils/openmc/constants.py new file mode 100644 index 0000000000..4ed83d235d --- /dev/null +++ b/src/utils/openmc/constants.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python + +"""Dictionaries of integer-to-string mapppings from openmc/src/constants.F90""" + +SURFACE_TYPES = {1: 'x-plane', + 2: 'y-plane', + 3: 'z-plane', + 4: 'plane', + 5: 'x-cylinder', + 6: 'y-cylinder', + 7: 'z-cylinder', + 8: 'sphere', + 9: 'x-cone', + 10: 'y-cone', + 11: 'z-cone'} + +BC_TYPES = {0: 'transmission', + 1: 'vacuum', + 2: 'reflective', + 3: 'periodic'} + +FILL_TYPES = {1: 'normal', + 2: 'fill', + 3: 'lattice'} + +LATTICE_TYPES = {1: 'rectangular', + 2: 'hexagonal'} + +ESTIMATOR_TYPES = {1: 'analog', + 2: 'tracklength'} + +FILTER_TYPES = {1: 'universe', + 2: 'material', + 3: 'cell', + 4: 'cellborn', + 5: 'surface', + 6: 'mesh', + 7: 'energy', + 8: 'energyout', + 9: 'distribcell'} + +SCORE_TYPES = {-1: 'flux', + -2: 'total', + -3: 'scatter', + -4: 'nu-scatter', + -5: 'scatter-n', + -6: 'scatter-pn', + -7: 'nu-scatter-n', + -8: 'nu-scatter-pn', + -9: 'transport', + -10: 'n1n', + -11: 'absorption', + -12: 'fission', + -13: 'nu-fission', + -14: 'kappa-fission', + -15: 'current', + -16: 'flux-yn', + -17: 'total-yn', + -18: 'scatter-yn', + -19: 'nu-scatter-yn', + -20: 'events', + 1: '(n,total)', + 2: '(n,elastic)', + 4: '(n,level)', + 11: '(n,2nd)', + 16: '(n,2n)', + 17: '(n,3n)', + 18: '(n,fission)', + 19: '(n,f)', + 20: '(n,nf)', + 21: '(n,2nf)', + 22: '(n,na)', + 23: '(n,n3a)', + 24: '(n,2na)', + 25: '(n,3na)', + 28: '(n,np)', + 29: '(n,n2a)', + 30: '(n,2n2a)', + 32: '(n,nd)', + 33: '(n,nt)', + 34: '(n,nHe-3)', + 35: '(n,nd2a)', + 36: '(n,nt2a)', + 37: '(n,4n)', + 38: '(n,3nf)', + 41: '(n,2np)', + 42: '(n,3np)', + 44: '(n,n2p)', + 45: '(n,npa)', + 91: '(n,nc)', + 101: '(n,disappear)', + 102: '(n,gamma)', + 103: '(n,p)', + 104: '(n,d)', + 105: '(n,t)', + 106: '(n,3He)', + 107: '(n,a)', + 108: '(n,2a)', + 109: '(n,3a)', + 111: '(n,2p)', + 112: '(n,pa)', + 113: '(n,t2a)', + 114: '(n,d2a)', + 115: '(n,pd)', + 116: '(n,pt)', + 117: '(n,da)', + 201: '(n,Xn)', + 202: '(n,Xgamma)', + 203: '(n,Xp)', + 204: '(n,Xd)', + 205: '(n,Xt)', + 206: '(n,X3He)', + 207: '(n,Xa)', + 444: '(damage)', + 649: '(n,pc)', + 699: '(n,dc)', + 749: '(n,tc)', + 799: '(n,3Hec)', + 849: '(n,tc)'} +SCORE_TYPES.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)}) +SCORE_TYPES.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)}) +SCORE_TYPES.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)}) +SCORE_TYPES.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)}) +SCORE_TYPES.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)}) +SCORE_TYPES.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)}) diff --git a/src/utils/openmc/element.py b/src/utils/openmc/element.py new file mode 100644 index 0000000000..202d22fbd3 --- /dev/null +++ b/src/utils/openmc/element.py @@ -0,0 +1,69 @@ +#!/usr/bin/env python + +from openmc.checkvalue import * + +class Element(object): + + def __init__(self, name='', xs=None): + + # Initialize class attributes + self._name = '' + self._xs = None + + # Set the Material class attributes + self.setName(name) + + if not xs is None: + self.setXS(xs) + + + def __eq__(self, element2): + + # check type + if not isinstance(element2, Element): + return False + + # Check name + if self._name != element2._name: + return False + + # Check xs + elif self._xs != element2._xs: + return False + + else: + return True + + + def __hash__(self): + hashable = list() + hashable.append(self._name) + hashable.append(self._xs) + return hash(tuple(hashable)) + + + def setName(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) + + self._name = name + + + def setXS(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) + + self._xs = xs + + + def __repr__(self): + + string = 'Element - {0}\n'.format(self._name) + string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) + return string diff --git a/src/utils/openmc/executor.py b/src/utils/openmc/executor.py new file mode 100644 index 0000000000..c6697ac92f --- /dev/null +++ b/src/utils/openmc/executor.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python + +from openmc.checkvalue import * +import subprocess +import os + + +FNULL = open(os.devnull, 'w') + + +class Executor(object): + + + def __init__(self): + self._working_directory = '.' + + + def setWorkingDirectory(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): + 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 plotGeometry(self, output=True): + + if output: + subprocess.check_call('openmc -p', shell=True, + cwd=self._working_directory) + else: + subprocess.check_call('openmc -p', shell=True, stdout=FNULL, + cwd=self._working_directory) + + + def runSimulation(self, particles=None, threads=None, geometry_debug=False, + restart_file=None, tracks=False, mpi_procs=1, output=True): + + post_args = ' ' + pre_args = '' + + if is_integer(particles) and particles > 0: + post_args += '-n {0} '.format(particles) + + if is_integer(threads) and threads > 0: + post_args += '-s {0} '.format(threads) + + if geometry_debug: + post_args += '-g ' + + if is_string(restart_file): + 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) + + command = pre_args + 'openmc ' + post_args + + if output: + subprocess.check_call(command, shell=True, + cwd=self._working_directory) + else: + subprocess.check_call(command, shell=True, stdout=FNULL, + cwd=self._working_directory) \ No newline at end of file diff --git a/src/utils/openmc/geometry.py b/src/utils/openmc/geometry.py new file mode 100644 index 0000000000..223981fa3f --- /dev/null +++ b/src/utils/openmc/geometry.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python + +import openmc +from openmc.clean_xml import * +from xml.etree import ElementTree as ET + +def reset_auto_ids(): + openmc.reset_auto_material_id() + openmc.reset_auto_surface_id() + openmc.reset_auto_cell_id() + openmc.reset_auto_universe_id() + + +class Geometry(object): + + def __init__(self): + + # Initialize Geometry class attributes + self._root_universe = None + self._offsets = dict() + + + def getOffset(self, path, filter_offset): + """ + Returns the corresponding location in the results array for a given + path and filter number. + + 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) + + filter_offset : int + An integer that specifies which offset map the filter is using + + """ + + # Return memoize'd offset if possible + if self._offsets.has_key((path, filter_offset)): + offset = self._offsets[(path, filter_offset)] + + # Begin recursive call to compute the offset starting with the base Universe + else: + offset = self._root_universe.getOffset(path, filter_offset) + self._offsets[(path, filter_offset)] = offset + + # Return the final offset + return offset + + + def getAllCells(self): + return self._root_universe.getAllCells() + + + def getAllUniverses(self): + return self._root_universe.getAllUniverses() + + + def getAllNuclides(self): + + nuclides = dict() + materials = self.getAllMaterials() + + for material in materials: + nuclides.update(material.getAllNuclides()) + + return nuclides + + + def getAllMaterials(self): + + material_cells = self.getAllMaterialCells() + materials = set() + + for cell in material_cells: + materials.add(cell._fill) + + return list(materials) + + + def getAllMaterialCells(self): + + all_cells = self.getAllCells() + material_cells = set() + + for cell_id, cell in all_cells.items(): + if cell._type == 'normal': + material_cells.add(cell) + + return list(material_cells) + + + def getAllMaterialUniverses(self): + + all_universes = self.getAllUniverses() + 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) + + return list(material_universes) + + + def setRootUniverse(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: + msg = 'Unable to add root Universe {0} to Geometry since it has ' \ + 'ID={1} instead of ID=0'.format(root_universe, root_universe._id) + raise ValueError(msg) + + self._root_universe = root_universe + + + +class GeometryFile(object): + + def __init__(self): + + # Initialize GeometryFile class attributes + self._geometry = None + self._geometry_file = ET.Element("geometry") + + + def setGeometry(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) + + self._geometry = geometry + + + def exportToXML(self): + + root_universe = self._geometry._root_universe + root_universe.createXMLSubElement(self._geometry_file) + + # Clean the indentation in the file to be user-readable + sort_xml_elements(self._geometry_file) + clean_xml_indentation(self._geometry_file) + + # Write the XML Tree to the materials.xml file + tree = ET.ElementTree(self._geometry_file) + tree.write("geometry.xml", xml_declaration=True, + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py new file mode 100644 index 0000000000..c6e9276afc --- /dev/null +++ b/src/utils/openmc/material.py @@ -0,0 +1,413 @@ +#!/usr/bin/env python + +import openmc +from openmc.checkvalue import * +from openmc.clean_xml import * +from xml.etree import ElementTree as ET +from collections import MappingView +from copy import deepcopy +import numpy as np + + +# A list of all IDs for all Materials created +MATERIAL_IDS = list() + +# 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 + MATERIAL_IDS = list() + + +# Units for density supported by OpenMC +DENSITY_UNITS = ['g/cm3', '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): + + def __init__(self, material_id=None, name=''): + + # Initialize class attributes + self._id = None + self._name = '' + self._density = None + self._density_units = '' + + # A dictionary of Nuclides + # Keys - Nuclide names + # Values - tuple (nuclide, percent, percent type) + self._nuclides = dict() + + # A dictionary of Elements + # Keys - Element names + # Values - tuple (element, percent, percent type) + self._elements = dict() + + # If specified, a list of tuples of (table name, xs identifier) + self._sab = list() + + # Set the Material class attributes + self.setId(material_id) + self.setName(name) + + + def setId(self, material_id=None): + + global MATERIAL_IDS + + # If the Material already has an ID, remove it from global list + if not self._id is None: + MATERIAL_IDS.remove(self._id) + + if material_id is None: + global AUTO_MATERIAL_ID + 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: + self._id = material_id + MATERIAL_IDS.append(material_id) + + + def setName(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 + + + def setDensity(self, units, density=NO_DENSITY): + + 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) + + 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) + + if density == NO_DENSITY and units is not 'sum': + msg = 'Unable to set the density Material ID={0} because '\ + 'a density must be set when not using sum unit'.format(self._id) + raise ValueError(msg) + + self._density = density + self._density_units = units + + + def addNuclide(self, nuclide, percent, percent_type='ao'): + + if not isinstance(nuclide, openmc.Nuclide): + msg = 'Unable to add an Nuclide to Material ID={0} with a ' \ + 'non-Nuclide value {1}'.format(self._id, nuclide) + raise ValueError(msg) + + elif not is_float(percent): + msg = 'Unable to add an Nuclide to Material ID={0} with a non-floating ' \ + 'non-floating point percent value {1}'.format(self._id, percent) + raise ValueError(msg) + + elif not percent_type in ['ao', 'wo', 'at/g-cm']: + msg = 'Unable to add an Nuclide to Material ID={0} with a ' \ + 'percent type {1}'.format(self._id, percent_type) + raise ValueError(msg) + + # Copy this Nuclide to separate it from the same Nuclide in other Materials + nuclide = deepcopy(nuclide) + + self._nuclides[nuclide._name] = (nuclide, percent, percent_type) + + + def removeNuclide(self, nuclide): + + if not isinstance(nuclide, openmc.Nuclide): + msg = 'Unable to remove a Nuclide {0} in Material ID={1} ' \ + 'since it is not a Nuclide'.format(self._id, nuclide) + raise ValueError(msg) + + # If the Material contains the Nuclide, delete it + if self._nuclides.has_key(nuclide._name): + del self._nuclides[nuclide._name] + + + def addElement(self, element, percent, percent_type='ao'): + + 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): + msg = 'Unable to add an Element to Material ID={0} with a non-floating ' \ + 'point percent value {1}'.format(self._id, percent) + raise ValueError(msg) + + if not percent_type in ['ao', 'wo']: + msg = 'Unable to add an Element to Material ID={0} with a percent ' \ + 'percent type {1}'.format(self._id, percent_type) + raise ValueError(msg) + + # Copy this Element to separate it from the same Element in other Materials + element = deepcopy(element) + + self._elements[element._name] = (element, percent, percent_type) + + + def removeElement(self, element): + + # If the Material contains the Element, delete it + if self._elements.has_key(element._name): + del self._elements[element._name] + + + def addSAlphaBeta(self, name, xs): + + if not is_string(name): + 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): + 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 getAllNuclides(self): + + nuclides = dict() + + for nuclide_name, nuclide_tuple in self._nuclides.items(): + nuclide = nuclide_tuple[0] + density = nuclide_tuple[1] + nuclides[nuclide._name] = (nuclide, density) + + 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) + + string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density) + string += ' [{0}]\n'.format(self._density_units) + + string += '{0: <16}'.format('\tS(a,b) Tables') + '\n' + + for sab in self._sab: + string += '{0: <16}{1}[{2}{3}]\n'.format('\tS(a,b)', '=\t', sab[0], sab[1]) + + string += '{0: <16}'.format('\tNuclides') + '\n' + + for nuclide in self._nuclides: + percent = self._nuclides[nuclide][1] + percent_type = self._nuclides[nuclide][2] + string += '{0: <16}'.format('\t{0}'.format(nuclide)) + string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type) + + string += '{0: <16}\n'.format('\tElements') + + for element in self._elements: + percent = self._nuclides[element][1] + percent_type = self._nuclides[element][2] + string += '{0: >16}'.format('\t{0}'.format(element)) + string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type) + + return string + + + def getNuclideXML(self, nuclide): + + xml_element = ET.Element("nuclide") + xml_element.set("name", nuclide[0]._name) + + if nuclide[2] is 'ao': + xml_element.set("ao", str(nuclide[1])) + else: + xml_element.set("wo", str(nuclide[1])) + + if not nuclide[0]._xs is None: + xml_element.set("xs", nuclide[0]._xs) + + return xml_element + + + def getElementXML(self, element): + + xml_element = ET.Element("element") + xml_element.set("name", str(element[0]._name)) + + if element[2] is 'ao': + xml_element.set("ao", str(element[1])) + else: + xml_element.set("wo", str(element[1])) + + return xml_element + + + def getNuclidesXML(self, nuclides): + + xml_elements = list() + + for nuclide in nuclides.values(): + xml_elements.append(self.getNuclideXML(nuclide)) + + return xml_elements + + + def getElementsXML(self, elements): + + xml_elements = list() + + for element in elements.values(): + xml_elements.append(self.getElementXML(element)) + + return xml_elements + + + def getMaterialXML(self): + + # Create Material XML element + element = ET.Element("material") + element.set("id", str(self._id)) + + # Create density XML subelement + subelement = ET.SubElement(element, "density") + if self._density_units is not 'sum': + subelement.set("value", str(self._density)) + subelement.set("units", self._density_units) + + # Create nuclide XML subelements + subelements = self.getNuclidesXML(self._nuclides) + for subelement in subelements: + element.append(subelement) + + # Create element XML subelements + subelements = self.getElementsXML(self._elements) + for subelement in subelements: + element.append(subelement) + + if len(self._sab) > 0: + for sab in self._sab: + subelement = ET.SubElement(element, "sab") + subelement.set("name", sab[0]) + subelement.set("xs", sab[1]) + + return element + + +class MaterialsFile(object): + + def __init__(self): + + # Initialize MaterialsFile class attributes + self._materials = list() + self._default_xs = None + self._materials_file = ET.Element("materials") + + + def addMaterial(self, material): + + if not isinstance(material, Material): + msg = 'Unable to add a non-Material {0} to the ' \ + 'MaterialsFile'.format(material) + raise ValueError(msg) + + self._materials.append(material) + + + def addMaterials(self, materials): + + if not isinstance(materials, (tuple, list, MappingView)): + msg = 'Unable to create OpenMC materials.xml file from {0} which ' \ + 'is not a Python tuple/list'.format(materials) + raise ValueError(msg) + + for material in materials: + self.addMaterial(material) + + + def removeMaterial(self, material): + + if not isinstance(material, Material): + msg = 'Unable to remove a non-Material {0} from the ' \ + 'MaterialsFile'.format(material) + raise ValueError(msg) + + self._materials.remove(material) + + + def setDefaultXS(self, xs): + + if not is_string(xs): + msg = 'Unable to set default xs to a non-string value'.format(xs) + raise ValueError(msg) + + self._default_xs = xs + + + def createMaterialSubelements(self): + + subelement = ET.SubElement(self._materials_file, "default_xs") + + if not self._default_xs is None: + subelement.text = self._default_xs + + for material in self._materials: + xml_element = material.getMaterialXML() + + if len(material._name) > 0: + self._materials_file.append(ET.Comment(material._name)) + + self._materials_file.append(xml_element) + + + def exportToXML(self): + + self.createMaterialSubelements() + + # Clean the indentation in the file to be user-readable + sort_xml_elements(self._materials_file) + clean_xml_indentation(self._materials_file) + + # Write the XML Tree to the materials.xml file + tree = ET.ElementTree(self._materials_file) + tree.write("materials.xml", xml_declaration=True, + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/nuclide.py b/src/utils/openmc/nuclide.py new file mode 100644 index 0000000000..1748acb46b --- /dev/null +++ b/src/utils/openmc/nuclide.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +from openmc.checkvalue import * + + +class Nuclide(object): + + def __init__(self, name='', xs=None): + + # Initialize class attributes + self._name = '' + self._xs = None + + # Set the Material class attributes + self.setName(name) + + if not xs is None: + self.setXS(xs) + + + def __eq__(self, nuclide2): + + # Check type + if not isinstance(nuclide2, Nuclide): + return False + + # Check name + elif self._name != nuclide2._name: + return False + + # Check xs + elif self._xs != nuclide2._xs: + return False + + else: + return True + + + def __hash__(self): + hashable = list() + hashable.append(self._name) + hashable.append(self._xs) + return hash(tuple(hashable)) + + + def setName(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) + + self._name = name + + + def setXS(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) + + self._xs = xs + + + def __repr__(self): + + string = 'Nuclide - {0}\n'.format(self._name) + string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) + return string diff --git a/src/utils/openmc/opencsg_compatible.py b/src/utils/openmc/opencsg_compatible.py new file mode 100644 index 0000000000..85d900927d --- /dev/null +++ b/src/utils/openmc/opencsg_compatible.py @@ -0,0 +1,820 @@ +#!/usr/bin/env python + +import openmc +import opencsg +import copy +import numpy as np + + +# A dictionary of all OpenMC Materials created +# Keys - Material IDs +# Values - Materials +OPENMC_MATERIALS = dict() + +# A dictionary of all OpenCSG Materials created +# Keys - Material IDs +# Values - Materials +OPENCSG_MATERIALS = dict() + +# A dictionary of all OpenMC Surfaces created +# Keys - Surface IDs +# Values - Surfaces +OPENMC_SURFACES = dict() + +# A dictionary of all OpenCSG Surfaces created +# Keys - Surface IDs +# Values - Surfaces +OPENCSG_SURFACES = dict() + +# A dictionary of all OpenMC Cells created +# Keys - Cell IDs +# Values - Cells +OPENMC_CELLS = dict() + +# A dictionary of all OpenCSG Cells created +# Keys - Cell IDs +# Values - Cells +OPENCSG_CELLS = dict() + +# A dictionary of all OpenMC Universes created +# Keys - Universes IDs +# Values - Universes +OPENMC_UNIVERSES = dict() + +# A dictionary of all OpenCSG Universes created +# Keys - Universes IDs +# Values - Universes +OPENCSG_UNIVERSES = dict() + +# A dictionary of all OpenMC Lattices created +# Keys - Lattice IDs +# Values - Lattices +OPENMC_LATTICES = dict() + +# A dictionary of all OpenCSG Lattices created +# Keys - Lattice IDs +# Values - Lattices +OPENCSG_LATTICES = dict() + + + +def get_opencsg_material(openmc_material): + + if not isinstance(openmc_material, openmc.Material): + msg = 'Unable to create an OpenCSG Material from {0} ' \ + 'which is not an OpenMC Material'.format(openmc_material) + raise ValueError(msg) + + global OPENCSG_MATERIALS + material_id = openmc_material._id + + # If this Material was already created, use it + if material_id in OPENCSG_MATERIALS.keys(): + return OPENCSG_MATERIALS[material_id] + + # Create an OpenCSG Material to represent this OpenMC Material + name = openmc_material._name + opencsg_material = opencsg.Material(material_id=material_id, name=name) + + # Add the OpenMC Material to the global collection of all OpenMC Materials + OPENMC_MATERIALS[material_id] = openmc_material + + # Add the OpenCSG Material to the global collection of all OpenCSG Materials + OPENCSG_MATERIALS[material_id] = opencsg_material + + return opencsg_material + + +def get_openmc_material(opencsg_material): + + if not isinstance(opencsg_material, opencsg.Material): + msg = 'Unable to create an OpenMC Material from {0} ' \ + 'which is not an OpenCSG Material'.format(opencsg_material) + raise ValueError(msg) + + global OPENMC_MATERIALS + material_id = opencsg_material._id + + # If this Material was already created, use it + if material_id in OPENMC_MATERIALS.keys(): + return OPENMC_MATERIALS[material_id] + + # Create an OpenMC Material to represent this OpenCSG Material + name = opencsg_material._name + openmc_material = openmc.Material(material_id=material_id, name=name) + + # Add the OpenMC Material to the global collection of all OpenMC Materials + OPENMC_MATERIALS[material_id] = openmc_material + + # Add the OpenCSG Material to the global collection of all OpenCSG Materials + OPENCSG_MATERIALS[material_id] = opencsg_material + + return openmc_material + + +def is_opencsg_surface_compatible(opencsg_surface): + + if not isinstance(opencsg_surface, opencsg.Surface): + msg = 'Unable to check if OpenCSG Surface is compatible' \ + 'since {0} is not a Surface'.format(opencsg_surface) + raise ValueError(msg) + + if opencsg_surface._type in ['x-squareprism', 'y-squareprism', 'z-squareprism']: + return False + else: + return True + + +def get_opencsg_surface(openmc_surface): + + if not isinstance(openmc_surface, openmc.Surface): + msg = 'Unable to create an OpenCSG Surface from {0} ' \ + 'which is not an OpenMC Surface'.format(openmc_surface) + raise ValueError(msg) + + global OPENCSG_SURFACES + surface_id = openmc_surface._id + + # If this Material was already created, use it + if surface_id in OPENCSG_SURFACES.keys(): + return OPENCSG_SURFACES[surface_id] + + # Create an OpenCSG Surface to represent this OpenMC Surface + name = openmc_surface._name + + # Correct for OpenMC's syntax for Surfaces dividing Cells + boundary = openmc_surface._bc_type + if boundary == 'transmission': + boundary = 'interface' + + opencsg_surface = None + + if openmc_surface._type == 'plane': + A = openmc_surface._coeffs['A'] + B = openmc_surface._coeffs['B'] + C = openmc_surface._coeffs['C'] + D = openmc_surface._coeffs['D'] + opencsg_surface = opencsg.Plane(surface_id, name, boundary, A, B, C, D) + + elif openmc_surface._type == 'x-plane': + x0 = openmc_surface._coeffs['x0'] + opencsg_surface = opencsg.XPlane(surface_id, name, boundary, x0) + + elif openmc_surface._type == 'y-plane': + y0 = openmc_surface._coeffs['y0'] + opencsg_surface = opencsg.YPlane(surface_id, name, boundary, y0) + + elif openmc_surface._type == 'z-plane': + z0 = openmc_surface._coeffs['z0'] + opencsg_surface = opencsg.ZPlane(surface_id, name, boundary, z0) + + elif openmc_surface._type == 'x-cylinder': + y0 = openmc_surface._coeffs['y0'] + z0 = openmc_surface._coeffs['z0'] + R = openmc_surface._coeffs['R'] + opencsg_surface = opencsg.XCylinder(surface_id, name, boundary, y0, z0, R) + + elif openmc_surface._type == 'y-cylinder': + x0 = openmc_surface._coeffs['x0'] + z0 = openmc_surface._coeffs['z0'] + R = openmc_surface._coeffs['R'] + opencsg_surface = opencsg.YCylinder(surface_id, name, boundary, x0, z0, R) + + elif openmc_surface._type == 'z-cylinder': + x0 = openmc_surface._coeffs['x0'] + y0 = openmc_surface._coeffs['y0'] + R = openmc_surface._coeffs['R'] + opencsg_surface = opencsg.ZCylinder(surface_id, name, boundary, x0, y0, R) + + # Add the OpenMC Surface to the global collection of all OpenMC Surfaces + OPENMC_SURFACES[surface_id] = openmc_surface + + # Add the OpenCSG Surface to the global collection of all OpenCSG Surfaces + OPENCSG_SURFACES[surface_id] = opencsg_surface + + return opencsg_surface + + +def get_openmc_surface(opencsg_surface): + + if not isinstance(opencsg_surface, opencsg.Surface): + msg = 'Unable to create an OpenMC Surface from {0} which ' \ + 'is not an OpenCSG Surface'.format(opencsg_surface) + raise ValueError(msg) + + global openmc_surface + surface_id = opencsg_surface._id + + # If this Surface was already created, use it + if surface_id in OPENMC_SURFACES.keys(): + return OPENMC_SURFACES[surface_id] + + # Create an OpenMC Surface to represent this OpenCSG Surface + name = opencsg_surface._name + + # Correct for OpenMC's syntax for Surfaces dividing Cells + boundary = opencsg_surface._boundary_type + if boundary == 'interface': + boundary = 'transmission' + + if opencsg_surface._type == 'plane': + A = opencsg_surface._coeffs['A'] + B = opencsg_surface._coeffs['B'] + C = opencsg_surface._coeffs['C'] + D = opencsg_surface._coeffs['D'] + openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name) + + elif opencsg_surface._type == 'x-plane': + x0 = opencsg_surface._coeffs['x0'] + openmc_surface = openmc.XPlane(surface_id, boundary, x0, name) + + elif opencsg_surface._type == 'y-plane': + y0 = opencsg_surface._coeffs['y0'] + openmc_surface = openmc.YPlane(surface_id, boundary, y0, name) + + elif opencsg_surface._type == 'z-plane': + z0 = opencsg_surface._coeffs['z0'] + openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name) + + elif opencsg_surface._type == 'x-cylinder': + y0 = opencsg_surface._coeffs['y0'] + z0 = opencsg_surface._coeffs['z0'] + R = opencsg_surface._coeffs['R'] + openmc_surface = openmc.XCylinder(surface_id, boundary, y0, z0, R, name) + + elif opencsg_surface._type == 'y-cylinder': + x0 = opencsg_surface._coeffs['x0'] + z0 = opencsg_surface._coeffs['z0'] + R = opencsg_surface._coeffs['R'] + openmc_surface = openmc.YCylinder(surface_id, boundary, x0, z0, R, name) + + elif opencsg_surface._type == 'z-cylinder': + x0 = opencsg_surface._coeffs['x0'] + y0 = opencsg_surface._coeffs['y0'] + R = opencsg_surface._coeffs['R'] + openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name) + + else: + msg = 'Unable to create an OpenMC Surface from an OpenCSG ' \ + 'Surface of type {0} since it is not a compatible ' \ + 'Surface type in OpenMC'.format(opencsg_surface._type) + raise ValueError(msg) + + + # Add the OpenMC Surface to the global collection of all OpenMC Surfaces + OPENMC_SURFACES[surface_id] = openmc_surface + + # Add the OpenCSG Surface to the global collection of all OpenCSG Surfaces + OPENCSG_SURFACES[surface_id] = opencsg_surface + + return openmc_surface + + +def get_compatible_opencsg_surfaces(opencsg_surface): + + if not isinstance(opencsg_surface, opencsg.Surface): + msg = 'Unable to create an OpenMC Surface from {0} which ' \ + 'is not an OpenCSG Surface'.format(opencsg_surface) + raise ValueError(msg) + + global OPENMC_SURFACES + surface_id = opencsg_surface._id + + # If this Surface was already created, use it + if surface_id in OPENMC_SURFACES.keys(): + return OPENMC_SURFACES[surface_id] + + # Create an OpenMC Surface to represent this OpenCSG Surface + name = opencsg_surface._name + boundary = opencsg_surface._boundary_type + + if opencsg_surface._type == 'x-squareprism': + y0 = opencsg_surface._coeffs['y0'] + z0 = opencsg_surface._coeffs['z0'] + R = opencsg_surface._coeffs['R'] + + # Create a list of the four planes we need + left = opencsg.YPlane(name=name, boundary=boundary, y0=y0-R) + right = opencsg.YPlane(name=name, boundary=boundary, y0=y0+R) + bottom = opencsg.ZPlane(name=name, boundary=boundary, z0=z0-R) + top = opencsg.ZPlane(name=name, boundary=boundary, z0=z0+R) + surfaces = [left, right, bottom, top] + + elif opencsg_surface._type == 'y-squareprism': + x0 = opencsg_surface._coeffs['x0'] + z0 = opencsg_surface._coeffs['z0'] + R = opencsg_surface._coeffs['R'] + + # Create a list of the four planes we need + left = opencsg.XPlane(name=name, boundary=boundary, x0=x0-R) + right = opencsg.XPlane(name=name, boundary=boundary, x0=x0+R) + bottom = opencsg.ZPlane(name=name, boundary=boundary, z0=z0-R) + top = opencsg.ZPlane(name=name, boundary=boundary, z0=z0+R) + surfaces = [left, right, bottom, top] + + elif opencsg_surface._type == 'z-squareprism': + x0 = opencsg_surface._coeffs['x0'] + y0 = opencsg_surface._coeffs['y0'] + R = opencsg_surface._coeffs['R'] + + # Create a list of the four planes we need + left = opencsg.XPlane(name=name, boundary=boundary, x0=x0-R) + right = opencsg.XPlane(name=name, boundary=boundary, x0=x0+R) + bottom = opencsg.YPlane(name=name, boundary=boundary, y0=y0-R) + top = opencsg.YPlane(name=name, boundary=boundary, y0=y0+R) + surfaces = [left, right, bottom, top] + + else: + msg = 'Unable to create a compatible OpenMC Surface an OpenCSG ' \ + 'Surface of type {0} since it already a compatible ' \ + 'Surface type in OpenMC'.format(opencsg_surface._type) + raise ValueError(msg) + + # Add the OpenMC Surface(s) to the global collection of all OpenMC Surfaces + OPENMC_SURFACES[surface_id] = surfaces + + # Add the OpenCSG Surface to the global collection of all OpenCSG Surfaces + OPENCSG_SURFACES[surface_id] = opencsg_surface + + return surfaces + + +def get_opencsg_cell(openmc_cell): + + if not isinstance(openmc_cell, openmc.Cell): + msg = 'Unable to create an OpenCSG Cell from {0} which ' \ + 'is not an OpenMC Cell'.format(openmc_cell) + raise ValueError(msg) + + global OPENCSG_CELLS + cell_id = openmc_cell._id + + # If this Cell was already created, use it + if cell_id in OPENCSG_CELLS.keys(): + return OPENCSG_CELLS[cell_id] + + # Create an OpenCSG Cell to represent this OpenMC Cell + name = openmc_cell._name + opencsg_cell = opencsg.Cell(cell_id, name) + + fill = openmc_cell._fill + + if (openmc_cell._type == 'normal'): + opencsg_cell.setFill(get_opencsg_material(fill)) + elif (openmc_cell._type == 'fill'): + opencsg_cell.setFill(get_opencsg_universe(fill)) + else: + opencsg_cell.setFill(get_opencsg_lattice(fill)) + + surfaces = openmc_cell._surfaces + + for surface_id in surfaces.keys(): + surface = surfaces[surface_id][0] + halfspace = surfaces[surface_id][1] + opencsg_cell.addSurface(get_opencsg_surface(surface), halfspace) + + # Add the OpenMC Cell to the global collection of all OpenMC Cells + OPENMC_CELLS[cell_id] = openmc_cell + + # Add the OpenCSG Cell to the global collection of all OpenCSG Cells + OPENCSG_CELLS[cell_id] = opencsg_cell + + return opencsg_cell + + +def get_compatible_opencsg_cells(opencsg_cell, opencsg_surface, halfspace): + + if not isinstance(opencsg_cell, opencsg.Cell): + msg = 'Unable to create compatible OpenMC Cell from {0} which ' \ + 'is not an OpenCSG Cell'.format(opencsg_cell) + raise ValueError(msg) + + elif not isinstance(opencsg_surface, opencsg.Surface): + msg = 'Unable to create compatible OpenMC Cell since {0} is ' \ + 'not an OpenCSG Surface'.format(opencsg_surface) + raise ValueError(msg) + + elif not halfspace in [-1, +1]: + msg = 'Unable to create compatible Cell since {0}' \ + 'is not a +/-1 halfspace'.format(halfspace) + raise ValueError(msg) + + # Initialize an empty list for the new compatible cells + compatible_cells = list() + + # SquarePrism Surfaces + if opencsg_surface._type in ['x-squareprism', 'y-squareprism', 'z-squareprism']: + + # Get the compatible Surfaces (XPlanes and YPlanes) + compatible_surfaces = get_compatible_opencsg_surfaces(opencsg_surface) + + opencsg_cell.removeSurface(opencsg_surface) + + # If Cell is inside SquarePrism, add "inside" of Surface halfspaces + if halfspace == -1: + opencsg_cell.addSurface(compatible_surfaces[0], +1) + opencsg_cell.addSurface(compatible_surfaces[1], -1) + opencsg_cell.addSurface(compatible_surfaces[2], +1) + opencsg_cell.addSurface(compatible_surfaces[3], -1) + compatible_cells.append(opencsg_cell) + + # 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 OpenCSG Cell with Surfaces compatible with OpenMC + clone = opencsg_cell.clone() + compatible_cells.append(clone) + + # Top left subcell - add left XPlane, top YPlane + if clone_id == 0: + clone.addSurface(compatible_surfaces[0], -1) + clone.addSurface(compatible_surfaces[3], +1) + + # Top center subcell - add top YPlane, left/right XPlanes + elif clone_id == 1: + clone.addSurface(compatible_surfaces[0], +1) + clone.addSurface(compatible_surfaces[1], -1) + clone.addSurface(compatible_surfaces[3], +1) + + # Top right subcell - add top YPlane, right XPlane + elif clone_id == 2: + clone.addSurface(compatible_surfaces[1], +1) + clone.addSurface(compatible_surfaces[3], +1) + + # Right center subcell - add right XPlane, top/bottom YPlanes + elif clone_id == 3: + clone.addSurface(compatible_surfaces[1], +1) + clone.addSurface(compatible_surfaces[3], -1) + clone.addSurface(compatible_surfaces[2], +1) + + # Bottom right subcell - add right XPlane, bottom YPlane + elif clone_id == 4: + clone.addSurface(compatible_surfaces[1], +1) + clone.addSurface(compatible_surfaces[2], -1) + + # Bottom center subcell - add bottom YPlane, left/right XPlanes + elif clone_id == 5: + clone.addSurface(compatible_surfaces[0], +1) + clone.addSurface(compatible_surfaces[1], -1) + clone.addSurface(compatible_surfaces[2], -1) + + # Bottom left subcell - add bottom YPlane, left XPlane + elif clone_id == 6: + clone.addSurface(compatible_surfaces[0], -1) + clone.addSurface(compatible_surfaces[2], -1) + + # Left center subcell - add left XPlane, top/bottom YPlanes + elif clone_id == 7: + clone.addSurface(compatible_surfaces[0], -1) + clone.addSurface(compatible_surfaces[3], -1) + clone.addSurface(compatible_surfaces[2], +1) + + # Remove redundant Surfaces from the Cells + for cell in compatible_cells: + cell.removeRedundantSurfaces() + + # Return the list of OpenMC compatible OpenCSG Cells + return compatible_cells + + +def make_opencsg_cells_compatible(opencsg_universe): + + if not isinstance(opencsg_universe, opencsg.Universe): + msg = 'Unable to make compatible OpenCSG Cells for {0} which ' \ + 'is not an OpenCSG Universe'.format(opencsg_universe) + raise ValueError(msg) + + # Check all OpenCSG Cells in this Universe for compatibility with OpenMC + opencsg_cells = opencsg_universe._cells + + for cell_id, opencsg_cell in opencsg_cells.items(): + + # Check each of the OpenCSG Surfaces for OpenMC compatibility + surfaces = opencsg_cell._surfaces + + for surface_id in surfaces.keys(): + surface = surfaces[surface_id][0] + halfspace = surfaces[surface_id][1] + + # If this Surface is not compatible with OpenMC, create compatible + # OpenCSG cells with a compatible version of this OpenCSG Surface + if not is_opencsg_surface_compatible(surface): + + # Get the one or more OpenCSG Cells that are compatible with OpenMC + # NOTE: This does not necessarily make the OpenCSG fully compatible. + # It only removes the incompatible Surface and replaces it with + # compatible OpenCSG Surface(s). The recursive call at the end + # of this block is necessary in the event that there are more + # incompatible Surfaces in this Cell that are not accounted for. + cells = get_compatible_opencsg_cells(opencsg_cell, surface, halfspace) + + # Remove the non-compatible OpenCSG Cell from the Universe + opencsg_universe.removeCell(opencsg_cell) + + # Add the compatible OpenCSG Cells to the Universe + opencsg_universe.addCells(cells) + + # Make recursive call to look at the updated state of the + # OpenCSG Universe and return + return make_opencsg_cells_compatible(opencsg_universe) + + # If all OpenCSG Cells in the OpenCSG Universe are compatible, return + return + + + +def get_openmc_cell(opencsg_cell): + + if not isinstance(opencsg_cell, opencsg.Cell): + msg = 'Unable to create an OpenMC Cell from {0} which ' \ + 'is not an OpenCSG Cell'.format(opencsg_cell) + raise ValueError(msg) + + global OPENMC_CELLS + cell_id = opencsg_cell._id + + # If this Cell was already created, use it + if cell_id in OPENMC_CELLS.keys(): + return OPENMC_CELLS[cell_id] + + # Create an OpenCSG Cell to represent this OpenMC Cell + name = opencsg_cell._name + openmc_cell = openmc.Cell(cell_id, name) + + fill = opencsg_cell._fill + + if (opencsg_cell._type == 'universe'): + openmc_cell.setFill(get_openmc_universe(fill)) + elif (opencsg_cell._type == 'lattice'): + openmc_cell.setFill(get_openmc_lattice(fill)) + else: + openmc_cell.setFill(get_openmc_material(fill)) + + surfaces = opencsg_cell._surfaces + + for surface_id in surfaces.keys(): + surface = surfaces[surface_id][0] + halfspace = surfaces[surface_id][1] + openmc_cell.addSurface(get_openmc_surface(surface), halfspace) + + # Add the OpenMC Cell to the global collection of all OpenMC Cells + OPENMC_CELLS[cell_id] = openmc_cell + + # Add the OpenCSG Cell to the global collection of all OpenCSG Cells + OPENCSG_CELLS[cell_id] = opencsg_cell + + return openmc_cell + + + +def get_opencsg_universe(openmc_universe): + + if not isinstance(openmc_universe, openmc.Universe): + msg = 'Unable to create an OpenCSG Universe from {0} which ' \ + 'is not an OpenMC Universe'.format(openmc_universe) + raise ValueError(msg) + + global OPENCSG_UNIVERSES + universe_id = openmc_universe._id + + # If this Universe was already created, use it + if universe_id in OPENCSG_UNIVERSES.keys(): + return OPENCSG_UNIVERSES[universe_id] + + # Create an OpenCSG Universe to represent this OpenMC Universe + name = openmc_universe._name + opencsg_universe = opencsg.Universe(universe_id, name) + + # Convert all OpenMC Cells in this Universe to OpenCSG Cells + openmc_cells = openmc_universe._cells + + for cell_id, openmc_cell in openmc_cells.items(): + opencsg_cell = get_opencsg_cell(openmc_cell) + opencsg_universe.addCell(opencsg_cell) + + # Add the OpenMC Universe to the global collection of all OpenMC Universes + OPENMC_UNIVERSES[universe_id] = openmc_universe + + # Add the OpenCSG Universe to the global collection of all OpenCSG Universes + OPENCSG_UNIVERSES[universe_id] = opencsg_universe + + return opencsg_universe + + +def get_openmc_universe(opencsg_universe): + + if not isinstance(opencsg_universe, opencsg.Universe): + msg = 'Unable to create an OpenMC Universe from {0} which ' \ + 'is not an OpenCSG Universe'.format(opencsg_universe) + raise ValueError(msg) + + global OPENMC_UNIVERSES + universe_id = opencsg_universe._id + + # If this Universe was already created, use it + if universe_id in OPENMC_UNIVERSES.keys(): + return OPENMC_UNIVERSES[universe_id] + + # Make all OpenCSG Cells and Surfaces in this Universe compatible with OpenMC + make_opencsg_cells_compatible(opencsg_universe) + + # Create an OpenMC Universe to represent this OpenCSg Universe + name = opencsg_universe._name + openmc_universe = openmc.Universe(universe_id, name) + + # Convert all OpenCSG Cells in this Universe to OpenMC Cells + opencsg_cells = opencsg_universe._cells + + for cell_id, opencsg_cell in opencsg_cells.items(): + openmc_cell = get_openmc_cell(opencsg_cell) + openmc_universe.addCell(openmc_cell) + + # Add the OpenMC Universe to the global collection of all OpenMC Universes + OPENMC_UNIVERSES[universe_id] = openmc_universe + + # Add the OpenCSG Universe to the global collection of all OpenCSG Universes + OPENCSG_UNIVERSES[universe_id] = opencsg_universe + + return openmc_universe + + +def get_opencsg_lattice(openmc_lattice): + + if not isinstance(openmc_lattice, openmc.Lattice): + msg = 'Unable to create an OpenCSG Lattice from {0} which ' \ + 'is not an OpenMC Lattice'.format(openmc_lattice) + raise ValueError(msg) + + global OPENCSG_LATTICES + lattice_id = openmc_lattice._id + + # If this Lattice was already created, use it + if lattice_id in OPENCSG_LATTICES.keys(): + return OPENCSG_LATTICES[lattice_id] + + # Create an OpenCSG Lattice to represent this OpenMC Lattice + name = openmc_lattice._name + dimension = openmc_lattice._dimension + width = openmc_lattice._width + lower_left = openmc_lattice._lower_left + universes = openmc_lattice._universes + + # Initialize an empty array for the OpenCSG nested Universes in this Lattice + universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), \ + dtype=opencsg.Universe) + + # Create OpenCSG Universes for each unique nested Universe in this Lattice + unique_universes = openmc_lattice.getUniqueUniverses() + + for universe_id, universe in unique_universes.items(): + unique_universes[universe_id] = get_opencsg_universe(universe) + + # Build the nested Universe array + for z in range(dimension[2]): + for y in range(dimension[1]): + for x in range(dimension[0]): + universe_id = universes[x][y][z]._id + universe_array[z][y][x] = unique_universes[universe_id] + + opencsg_lattice = opencsg.Lattice(lattice_id, name) + opencsg_lattice.setDimension(dimension) + opencsg_lattice.setWidth(width) + opencsg_lattice.setUniverses(universe_array) + + offset = np.array(lower_left, dtype=np.float64) - \ + ((np.array(width, dtype=np.float64) * \ + np.array(dimension, dtype=np.float64))) / -2.0 + opencsg_lattice.setOffset(offset) + + # Add the OpenMC Lattice to the global collection of all OpenMC Lattices + OPENMC_LATTICES[lattice_id] = openmc_lattice + + # Add the OpenCSG Lattice to the global collection of all OpenCSG Lattices + OPENCSG_LATTICES[lattice_id] = opencsg_lattice + + return opencsg_lattice + + +def get_openmc_lattice(opencsg_lattice): + + if not isinstance(opencsg_lattice, opencsg.Lattice): + msg = 'Unable to create an OpenMC Lattice from {0} which ' \ + 'is not an OpenCSG Lattice'.format(opencsg_lattice) + raise ValueError(msg) + + global OPENMC_LATTICES + lattice_id = opencsg_lattice._id + + # If this Lattice was already created, use it + if lattice_id in OPENMC_LATTICES.keys(): + return OPENMC_LATTICES[lattice_id] + + dimension = opencsg_lattice._dimension + width = opencsg_lattice._width + offset = opencsg_lattice._offset + universes = opencsg_lattice._universes + + # Initialize an empty array for the OpenMC nested Universes in this Lattice + universe_array = np.ndarray(tuple(np.array(dimension)), \ + dtype=openmc.Universe) + + # Create OpenMC Universes for each unique nested Universe in this Lattice + unique_universes = opencsg_lattice.getUniqueUniverses() + + for universe_id, universe in unique_universes.items(): + unique_universes[universe_id] = get_openmc_universe(universe) + + # Build the nested Universe array + for z in range(dimension[2]): + for y in range(dimension[1]): + for x in range(dimension[0]): + universe_id = universes[z][y][x]._id + universe_array[x][y][z] = unique_universes[universe_id] + + openmc_lattice = openmc.Lattice(lattice_id=lattice_id) + openmc_lattice.setDimension(dimension) + openmc_lattice.setWidth(width) + openmc_lattice.setUniverses(universe_array) + + lower_left = np.array(offset, dtype=np.float64) + \ + ((np.array(width, dtype=np.float64) * \ + np.array(dimension, dtype=np.float64))) / -2.0 + openmc_lattice.setLowerLeft(lower_left) + + # Add the OpenMC Lattice to the global collection of all OpenMC Lattices + OPENMC_LATTICES[lattice_id] = openmc_lattice + + # Add the OpenCSG Lattice to the global collection of all OpenCSG Lattices + OPENCSG_LATTICES[lattice_id] = opencsg_lattice + + return openmc_lattice + + +def get_opencsg_geometry(openmc_geometry): + + if not isinstance(openmc_geometry, openmc.Geometry): + msg = 'Unable to get OpenCSG geometry from {0} which is ' \ + 'not an OpenMC Geometry object'.format(openmc_geometry) + raise ValueError(msg) + + # Clear dictionaries and auto-generated IDs + OPENMC_SURFACES.clear() + OPENCSG_SURFACES.clear() + OPENMC_CELLS.clear() + OPENCSG_CELLS.clear() + OPENMC_UNIVERSES.clear() + OPENCSG_UNIVERSES.clear() + OPENMC_LATTICES.clear() + OPENCSG_LATTICES.clear() + + opencsg.geometry.reset_auto_ids() + + openmc_root_universe = openmc_geometry._root_universe + opencsg_root_universe = get_opencsg_universe(openmc_root_universe) + + opencsg_geometry = opencsg.Geometry() + opencsg_geometry.setRootUniverse(opencsg_root_universe) + opencsg_geometry.initializeCellOffsets() + + return opencsg_geometry + + +def get_openmc_geometry(opencsg_geometry): + + if not isinstance(opencsg_geometry, opencsg.Geometry): + msg = 'Unable to get OpenMC geometry from {0} which is ' \ + 'not an OpenCSG Geometry object'.format(opencsg_geometry) + raise ValueError(msg) + + # Deep copy the goemetry since it may be modified to make all Surfaces + # compatible with OpenMC's specifications + opencsg_geometry = copy.deepcopy(opencsg_geometry) + + # Update Cell bounding boxes in Geometry + opencsg_geometry.updateBoundingBoxes() + + # Clear dictionaries and auto-generated ID + OPENMC_SURFACES.clear() + OPENCSG_SURFACES.clear() + OPENMC_CELLS.clear() + OPENCSG_CELLS.clear() + OPENMC_UNIVERSES.clear() + OPENCSG_UNIVERSES.clear() + OPENMC_LATTICES.clear() + OPENCSG_LATTICES.clear() + + openmc.reset_auto_ids() + + opencsg_root_universe = opencsg_geometry._root_universe + openmc_root_universe = get_openmc_universe(opencsg_root_universe) + + openmc_geometry = openmc.Geometry() + openmc_geometry.setRootUniverse(openmc_root_universe) + + return openmc_geometry diff --git a/src/utils/openmc/plots.py b/src/utils/openmc/plots.py new file mode 100644 index 0000000000..dab68095ee --- /dev/null +++ b/src/utils/openmc/plots.py @@ -0,0 +1,434 @@ +#!/usr/bin/env python + +from openmc.checkvalue import * +from openmc.clean_xml import * +from xml.etree import ElementTree as ET +import numpy as np + + +# 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 + + +BASES = ['xy', 'xz', 'yz'] + + +class Plot(object): + + def __init__(self, plot_id=None, name=''): + + # Initialize Plot class attributes + self._id = None + self._name = '' + self._width = [4.0, 4.0] + self._pixels = [1000, 1000] + self._origin = [0., 0., 0.] + self._filename = 'plot' + self._color = 'cell' + self._type = 'slice' + self._basis = 'xy' + self._background = [0, 0, 0] + self._mask_components = None + self._mask_background = None + self._col_spec = None + + self.setId(plot_id) + self.setName(name) + + + def setId(self, plot_id=None): + + 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: + self._id = plot_id + + + def setName(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 + + + def setWidth(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) + + self._width = width + + + def setOrigin(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) + + self._origin = origin + + + def setPixels(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) + + 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) + + self._pixels = pixels + + + def setFilename(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) + + self._filename = filename + + + def setColor(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) + + self._color = color + + + def setType(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 setBasis(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) + + self._basis = basis + + + def setBackground(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) + + 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) + + self._background = background + + + def setColSpec(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) + + for key in col_spec.keys(): + + 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: + 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]) + raise ValueError(msg) + + elif len(col_spec[key]) != 3: + msg = 'Unable to create Plot ID={0} with col_spec RGB values of ' \ + 'length {1} since 3 values must be ' \ + 'input'.format(self._id, len(col_spec[key])) + raise ValueError(msg) + + self._col_spec = col_spec + + + def setMaskComponents(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) + + 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) + + self._mask_components = mask_components + + + def setMaskBackground(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) + + 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) + + 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) + string += '{0: <16}{1}{2}\n'.format('\tFilename', '=\t', self._filename) + string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) + string += '{0: <16}{1}{2}\n'.format('\tBasis', '=\t', self._basis) + string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._width) + string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin) + string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin) + string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color) + string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t', self._mask_components) + string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t', self._mask_background) + string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec) + return string + + + def getPlotXML(self): + + element = ET.Element("plot") + element.set("id", str(self._id)) + element.set("filename", self._filename) + element.set("color", self._color) + element.set("type", self._type) + + if self._type is 'slice': + element.set("basis", self._basis) + + subelement = ET.SubElement(element, "origin") + text = '' + for coord in self._origin: + text += str(coord) + ' ' + subelement.text = text.rstrip(' ') + + subelement = ET.SubElement(element, "width") + text = '' + for dim in self._width: + text += str(dim) + ' ' + subelement.text = text.rstrip(' ') + + subelement = ET.SubElement(element, "pixels") + text = '' + for dim in self._pixels: + text += str(dim) + ' ' + subelement.text = text.rstrip(' ') + + if not self._mask_background is 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: + + for key in self._col_spec.keys(): + 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])) + + if not self._mask_components is 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])) + + return element + + +class PlotsFile(object): + + def __init__(self): + + # Initialize PlotsFile class attributes + self._plots = list() + self._plots_file = ET.Element("plots") + + + def addPlot(self, plot): + + if not isinstance(plot, Plot): + msg = 'Unable to add a non-Plot {0} to the PlotsFile'.format(plot) + raise ValueError(msg) + + self._plots.append(plot) + + + def removePlot(self, plot): + self._plots.remove(plot) + + + def createPlotSubelements(self): + + for plot in self._plots: + + xml_element = plot.getPlotXML() + + if len(plot._name) > 0: + self._plots_file.append(ET.Comment(plot._name)) + + self._plots_file.append(xml_element) + + + def exportToXML(self): + + self.createPlotSubelements() + + # Clean the indentation in the file to be user-readable + clean_xml_indentation(self._plots_file) + + # Write the XML Tree to the plots.xml file + tree = ET.ElementTree(self._plots_file) + tree.write("plots.xml", xml_declaration=True, + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/settings.py b/src/utils/openmc/settings.py new file mode 100644 index 0000000000..95124562e8 --- /dev/null +++ b/src/utils/openmc/settings.py @@ -0,0 +1,1053 @@ +#!/usr/bin/env python + +from openmc.checkvalue import * +from openmc.clean_xml import * +from xml.etree import ElementTree as ET +import numpy as np +import multiprocessing + + +class SettingsFile(object): + + def __init__(self): + + # Eigenvalue subelement + self._batches = None + self._generations_per_batch = None + self._inactive = None + self._particles = None + + # Source subelement + self._source_subelement = None + self._source_file = None + self._source_space_type = None + self._source_space_params = None + self._source_angle_type = None + self._source_angle_params = None + self._source_energy_type = None + self._source_energy_params = None + + self._confidence_intervals = None + self._cross_sections = None + self._energy_grid = None + self._ptables = None + self._run_cmfd = None + self._seed = None + self._survival_biasing = None + + # Entropy subelement + self._entropy_dimension = None + self._entropy_lower_left = None + self._entropy_upper_right = None + + self._output = None + self._output_path = None + + # Statepoint subelement + self._statepoint_batches = None + self._statepoint_interval = None + self._sourcepoint_batches = None + self._sourcepoint_interval = None + self._sourcepoint_separate = None + self._sourcepoint_write = None + self._sourcepoint_overwrite = None + + self._threads = None + self._no_reduce = None + + self._verbosity = None + + self._trace = None + self._track = None + + # Cutoff subelement + self._weight = None + self._weight_avg = None + + # Uniform fission source subelement + self._ufs_dimension = 1 + self._ufs_lower_left = None + self._ufs_upper_right = None + + self._settings_file = ET.Element("settings") + self._eigenvalue_element = None + self._source_element = None + + + def setBatches(self, batches): + + if not is_integer(batches): + msg = 'Unable to set batches to a non-integer value {0}'.format(batches) + raise ValueError(msg) + + if batches <= 0: + msg = 'Unable to set batches to a negative value {0}'.format(batches) + raise ValueError(msg) + + self._batches = batches + + + def setGenerationsPerBatch(self, generations_per_batch): + + if not is_integer(generations_per_batch): + msg = 'Unable to set generations per batch to a non-integer ' \ + 'value {0}'.format(generations_per_batch) + raise ValueError(msg) + + if generations_per_batch <= 0: + msg = 'Unable to set generations per batch to a negative ' \ + 'value {0}'.format(generations_per_batch) + raise ValueError(msg) + + self._generations_per_batch = generations_per_batch + + + def setInactive(self, inactive): + + if not is_integer(inactive): + msg = 'Unable to set inactive batches to a non-integer ' \ + 'value {0}'.format(inactive) + raise ValueError(msg) + + if inactive <= 0: + msg = 'Unable to set inactive batches to a negative ' \ + 'value {0}'.format(inactive) + raise ValueError(msg) + + self._inactive = inactive + + + def setParticles(self, particles): + + if not is_integer(particles): + msg = 'Unable to set particles to a non-integer ' \ + 'value {0}'.format(particles) + raise ValueError(msg) + + if particles <= 0: + msg = 'Unable to set particles to a negative value {0}'.format(particles) + raise ValueError(msg) + + self._particles = particles + + + def setSourceFile(self, source_file): + + if not is_string(source_file): + msg = 'Unable to set source file to a non-string ' \ + 'value {0}'.format(source_file) + raise ValueError(msg) + + self._source_file = source_file + + + def setSourceSpace(self, type, params): + + if not is_string(type): + msg = 'Unable to set source space type to a non-string ' \ + 'value {0}'.format(type) + raise ValueError(msg) + + elif not type in ['box', 'point']: + msg = 'Unable to set source space type to {0} since it is not ' \ + 'box or point'.format(type) + raise ValueError(msg) + + elif not isinstance(params, (tuple, list, np.ndarray)): + msg = 'Unable to set source space parameters to {0} since it is ' \ + 'not a Python tuple, list or NumPy array'.format(params) + raise ValueError(msg) + + elif len(params) != 6: + msg = 'Unable to set source space parameters to {0} since it does ' \ + 'not contain 6 values'.format(params) + raise ValueError(msg) + + for param in params: + + if not is_integer(param) and not is_float(param): + msg = 'Unable to set source space parameters to {0} since it ' \ + 'is not an integer or floating point value'.format(param) + raise ValueError(msg) + + self._source_space_type = type + self._source_space_params = params + + + def setSourceAngle(self, type, params=[]): + + if not is_string(type): + msg = 'Unable to set source angle type to a non-string ' \ + 'value {0}'.format(type) + raise ValueError(msg) + + elif not type in ['isotropic', 'monodirectional']: + msg = 'Unable to set source angle type to {0} since it is not ' \ + 'isotropic or monodirectional'.format(type) + raise ValueError(msg) + + elif not isinstance(params, (tuple, list, np.ndarray)): + msg = 'Unable to set source angle parameters to {0} since it is ' \ + 'not a Python list/tuple or NumPy array'.format(params) + raise ValueError(msg) + + elif type is 'isotropic' and not params is None: + msg = 'Unable to set source angle parameters since they are not ' \ + 'it is not supported for isotropic type sources' + raise ValueError(msg) + + elif type is 'monodirectional' and len(params) != 3: + msg = 'Unable to set source angle parameters to {0} since 3 ' \ + 'parameters are required for monodirectional sources'.format(params) + raise ValueError(msg) + + for param in params: + + if not is_integer(param) and not is_float(param): + msg = 'Unable to set source angle parameters to {0} since it ' \ + 'is not an integer or floating point value'.format(param) + raise ValueError(msg) + + self._source_angle_type = type + self._source_angle_params = params + + + def setSourceEnergy(self, type, params=[]): + + if not is_string(type): + msg = 'Unable to set source energy type to a non-string ' \ + 'value {0}'.format(type) + raise ValueError(msg) + + elif not type in ['monoenergetic', 'watt', 'maxwell']: + msg = 'Unable to set source energy type to {0} since it is not ' \ + 'monoenergetic, watt or maxwell'.format(type) + raise ValueError(msg) + + elif not isinstance(params, (tuple, list, np.ndarray)): + msg = 'Unable to set source energy parameters to {0} since it is not ' \ + 'a Python list/tuple or NumPy array'.format(params) + raise ValueError(msg) + + elif type is 'monoenergetic' and not len(params) != 1: + msg = 'Unable to set source energy parameters to {0} since 1 ' \ + 'paramater is required for monenergetic sources'.format(params) + raise ValueError(msg) + + elif type is 'watt' and len(params) != 2: + msg = 'Unable to set source energy parameters to {0} since 2 ' \ + 'parameters are required for monoenergetic sources'.format(params) + raise ValueError(msg) + + elif type is 'maxwell' and len(params) != 2: + msg = 'Unable to set source energy parameters to {0} since 1 ' \ + 'parameter is required for maxwell sources'.format(params) + raise ValueError(msg) + + for param in params: + + if not is_integer(param) and not is_float(param): + msg = 'Unable to set source energy parameters to {0} since it ' \ + 'is not an integer or floating point value'.format(param) + raise ValueError(msg) + + self._source_energy_type = type + self._source_energy_params = params + + + def setOutput(self, output): + + if not isinstance(output, dict): + msg = 'Unable to set output to {0} which is not a Python dictionary ' \ + 'of string keys and boolean values'.format(output) + raise ValueError(msg) + + + for element in output.keys(): + + if not element in ['summary', 'cross_sections', 'tallies']: + msg = 'Unable to set output to {0} which is unsupported by ' \ + 'OpenMC'.format(element) + raise ValueError(msg) + + if not isinstance(output[element], (bool, np.bool)): + msg = 'Unable to set output for {0} to a non-boolean ' \ + 'value {1}'.format(element, output[element]) + raise ValueError(msg) + + self._output = output + + + def setOutputPath(self, output_path): + + if not is_string(output_path): + msg = 'Unable to set output path to non-string ' \ + 'value {0}'.format(output_path) + raise ValueError(msg) + + self._output_path = output_path + + + def setVerbosity(self, verbosity): + + if not is_integer(verbosity): + msg = 'Unable to set verbosity to non-integer value {0}'.format(verbosity) + raise ValueError(msg) + + if verbosity < 1 or verbosity > 10: + msg = 'Unable to set verbosity to {0} which is not between ' \ + '1 and 10'.format(verbosity) + raise ValueError(msg) + + self._verbosity = verbosity + + + def setStatepointBatches(self, batches): + + if not isinstance(batches, (tuple, list, np.ndarray)): + msg = 'Unable to set statepoint batches to {0} which is not a Python ' \ + 'tuple/list or NumPy array'.format(batches) + raise ValueError(msg) + + for batch in batches: + + if not is_integer(batch): + msg = 'Unable to set statepoint batches with non-integer ' \ + 'value {0}'.format(batch) + raise ValueError(msg) + + if batch <= 0: + msg = 'Unable to set statepoint batches with {0} which is less than ' \ + 'or equal to zero'.format(batch) + raise ValueError(msg) + + self._statepoint_batches = batches + + + def setStatepointInterval(self, interval): + + if not is_integer(interval): + msg = 'Unable to set statepoint interval to non-integer ' \ + 'value {0}'.format(interval) + raise ValueError(msg) + + self._statepoint_interval = interval + + + def setSourcepointBatches(self, batches): + + if not isinstance(batches, (tuple, list, np.ndarray)): + msg = 'Unable to set sourcepoint batches to {0} which is not a Python ' \ + 'tuple/list or NumPy array'.format(batches) + raise ValueError(msg) + + for batch in batches: + + if not is_integer(batch): + msg = 'Unable to set sourcepoint batches with non-integer ' \ + 'value {0}'.format(batch) + raise ValueError(msg) + + if batch <= 0: + msg = 'Unable to set sourcepoint batches with {0} which is less than ' \ + 'or equal to zero'.format(batch) + raise ValueError(msg) + + self._sourcepoint_batches = batches + + + def setSourcepointInterval(self, interval): + + if not is_integer(interval): + msg = 'Unable to set sourcepoint interval to non-integer ' \ + 'value {0}'.format(interval) + raise ValueError(msg) + + self._sourcepoint_interval = interval + + + def setSourcepointSeparate(self, source_separate): + + if not isinstance(source_separate, (bool, np.bool)): + msg = 'Unable to set sourcepoint separate to non-boolean ' \ + 'value {0}'.format(source_separate) + raise ValueError(msg) + + self._sourcepoint_separate = source_separate + + + def setSourcepointWrite(self, source_write): + + if not isinstance(source_write, (bool, np.bool)): + msg = 'Unable to set sourcepoint write to non-boolean ' \ + 'value {0}'.format(source_write) + raise ValueError(msg) + + self._sourcepoint_write = source_write + + + def setSourcepointOverwrite(self, source_overwrite): + + if not isinstance(source_overwrite, (bool, np.bool)): + msg = 'Unable to set sourcepoint overwrite to non-boolean ' \ + 'value {0}'.format(source_overwrite) + raise ValueError(msg) + + self._sourcepoint_overwrite = source_overwrite + + + def setConfidenceIntervals(self, confidence_intervals): + + if not isinstance(confidence_intervals, (bool, np.bool)): + msg = 'Unable to set confidence interval to non-boolean ' \ + 'value {0}'.format(confidence_intervals) + raise ValueError(msg) + + self._confidence_intervals = confidence_intervals + + + def setCrossSections(self, cross_sections): + + if not is_string(cross_sections): + msg = 'Unable to set cross sections to non-string ' \ + 'value {0}'.format(cross_sections) + raise ValueError(msg) + + self._cross_sections = cross_sections + + + def setEnergyGrid(self, energy_grid): + + if not energy_grid in ['union', 'nuclide']: + msg = 'Unable to set energy grid to {0} which is neither ' \ + 'union nor nuclide'.format(energy_grid) + raise ValueError(msg) + + self._energy_grid = energy_grid + + + def setPTables(self, ptables): + + if not isinstance(ptables, (bool, np.bool)): + msg = 'Unable to set ptables to non-boolean value {0}'.format(ptables) + raise ValueError(msg) + + self._ptables = ptables + + + def setRunCMFD(self, run_cmfd): + + if not isinstance(run_cmfd, (bool, np.bool)): + msg = 'Unable to set run_cmfd to non-boolean value {0}'.format(run_cmfd) + raise ValueError(msg) + + self._run_cmfd = run_cmfd + + + def setSeed(self, seed): + + if not is_integer(seed): + msg = 'Unable to set seed to non-integer value {0}'.format(seed) + raise ValueError(msg) + + elif seed <= 0: + msg = 'Unable to set seed to non-positive integer {0}'.format(seed) + raise ValueError(msg) + + self._seed = seed + + + def setSurvivalBiasing(self, survival_biasing): + + if not isinstance(survival_biasing, (bool, np.bool)): + msg = 'Unable to set survival biasing to non-boolean ' \ + 'value {0}'.format(survival_biasing) + raise ValueError(msg) + + self._survival_biasing = survival_biasing + + + def setWeight(self, weight, weight_avg): + + if not is_float(weight): + msg = 'Unable to set weight cutoff to non-floating point ' \ + 'value {0}'.format(weight) + raise ValueError(msg) + + elif not is_float(weight_avg): + msg = 'Unable to set weight avg. to non-floating point ' \ + 'value {0}'.format(weight_avg) + raise ValueError(msg) + + elif weight < 0.0: + msg = 'Unable to set weight cutoff to negative value {0}'.format(weight) + raise ValueError(msg) + + elif weight_avg < 0.0: + msg = 'Unable to set weight avg. to negative value {0}'.format(weight_avg) + raise ValueError(msg) + + self._weight = weight + self._weight_avg = weight_avg + + + def setEntropyDimension(self, dimension): + + if not isinstance(dimension, (tuple, list)): + msg = 'Unable to set entropy mesh dimension to {0} which is ' \ + 'not a Python tuple or list'.format(dimension) + raise ValueError(msg) + + elif len(dimension) < 3 or len(dimension) > 3: + msg = 'Unable to set entropy mesh dimension to {0} which is ' \ + 'not a set of 3 integer dimensions'.format(dimension) + raise ValueError(msg) + + for dim in dimension: + + if not is_integer(dim) and not is_float(dim): + msg = 'Unable to set entropy mesh dimension to a ' \ + 'non-integer or floating point value {0}'.format(dim) + raise ValueError(msg) + + self._entropy_dimension = dimension + + + def setEntropyLowerLeft(self, lower_left): + + if not isinstance(lower_left, (tuple, list)): + msg = 'Unable to set entropy mesh lower left corner to {0} which is ' \ + 'not a Python tuple or list'.format(lower_left) + raise ValueError(msg) + + elif len(lower_left) < 3 or len(lower_left) > 3: + msg = 'Unable to set entropy mesh lower left corner to {0} which is ' \ + 'not a 3D point'.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 entropy mesh lower left corner to a ' \ + 'non-integer or floating point value {0}'.format(coord) + raise ValueError(msg) + + self._entropy_lower_left = lower_left + + + def setEntropyUpperRight(self, upper_right): + + if not isinstance(upper_right, (tuple, list)): + msg = 'Unable to set entropy mesh upper right corner to {0} which is ' \ + 'not a Python tuple or list'.format(upper_right) + raise ValueError(msg) + + elif len(upper_right) < 3 or len(upper_right) > 3: + msg = 'Unable to set entropy mesh upper right corner to {0} which is ' \ + 'not a 3D point'.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 entropy mesh upper right corner to a ' \ + 'non-integer or floating point value {0}'.format(coord) + raise ValueError(msg) + + self._entropy_upper_right = upper_right + + + def setNoReduce(self, no_reduce): + + if not isinstance(no_reduce, (bool, np.bool)): + msg = 'Unable to set the no_reduce to a non-boolean ' \ + 'value {0}'.format(no_reduce) + raise ValueError(msg) + + self._no_reduce = no_reduce + + + def setThreads(self, threads): + + if not is_integer(threads): + msg = 'Unable to set the threads to a non-integer ' \ + 'value {0}'.format(threads) + raise ValueError(msg) + + elif threads <=0: + msg = 'Unable to set the threads to a negative value {0}'.format(threads) + raise ValueError(msg) + + self._threads = threads + + + def setTrace(self, trace): + + if not isinstance(trace, (list, tuple)): + msg = 'Unable to set the trace to {0} which is not a Python tuple ' \ + 'or list'.format(trace) + raise ValueError(msg) + + elif len(trace) != 3: + msg = 'Unable to set the trace to {0} since it does not contain ' \ + '3 elements - batch, generation, and particle'.format(trace) + raise ValueError(msg) + + elif trace[0] < 1: + msg = 'Unable to set the trace batch to {0} since it must be ' \ + 'greater than or equal to 1'.format(trace[0]) + raise ValueError(msg) + + elif trace[1] < 1: + msg = 'Unable to set the trace generation to {0} since it must be ' \ + 'greater than or equal to 1'.format(trace[1]) + raise ValueError(msg) + + elif trace[2] < 1: + msg = 'Unable to set the trace particle to {0} since it must be ' \ + 'greater than or equal to 1'.format(trace[2]) + raise ValueError(msg) + + self._trace = trace + + + def setTrack(self, track): + + if not isinstance(track, (list, tuple)): + msg = 'Unable to set the track to {0} which is not a Python tuple ' \ + 'or list'.format(track) + raise ValueError(msg) + + elif len(track) != 3: + msg = 'Unable to set the track to {0} since it does not contain ' \ + '3 elements - batch, generation, and particle'.format(track) + raise ValueError(msg) + + + elif track[0] < 1: + msg = 'Unable to set the track batch to {0} since it must be ' \ + 'greater than or equal to 1'.format(track[0]) + raise ValueError(msg) + + elif track[1] < 1: + msg = 'Unable to set the track generation to {0} since it must be ' \ + 'greater than or equal to 1'.format(track[1]) + raise ValueError(msg) + + elif track[2] < 1: + msg = 'Unable to set the track particle to {0} since it must be ' \ + 'greater than or equal to 1'.format(track[2]) + raise ValueError(msg) + + self._track = track + + + def setUFSDimension(self, dimension): + + if not is_integer(dimension) and not is_float(dimension): + msg = 'Unable to set UFS dimension to non-integer or non-floating ' \ + 'point value {0}'.format(dimension) + raise ValueError(msg) + + elif dimension < 1: + msg = 'Unable to set UFS dimension to value {0} which is less than ' \ + 'one'.format(dimension) + raise ValueError(msg) + + self._ufs_dimension = dimension + + + def setUFSLowerLeft(self, lower_left): + + if not isinstance(lower_left, (tuple, list, np.ndarray)): + msg = 'Unable to set UFS mesh lower left corner to {0} which is ' \ + 'not a Python tuple or list'.format(lower_left) + raise ValueError(msg) + + elif len(lower_left) < 3 or len(lower_left) > 3: + msg = 'Unable to set UFS mesh lower left corner to {0} which is not ' \ + 'a 3D point'.format(lower_left) + raise ValueError(msg) + + self._ufs_lower_left = lower_left + + + def setUFSUpperRight(self, upper_right): + + if not isinstance(upper_right, tuple) and not isinstance(upper_right, list): + msg = 'Unable to set UFs mesh upper right corner to {0} which is ' \ + 'not a Python tuple or list'.format(upper_right) + raise ValueError(msg) + + if len(upper_right) < 3 or len(upper_right) > 3: + msg = 'Unable to set UFS mesh upper right corner to {0} which is not ' \ + 'a 3D point'.format(upper_right) + raise ValueError(msg) + + self._ufs_upper_right = upper_right + + + def createEigenvalueSubelement(self): + + self.createParticlesSubelement() + self.createBatchesSubelement() + self.createInactiveSubelement() + self.createGenerationsPerBatchSubelement() + + + def createBatchesSubelement(self): + + if not self._batches is None: + + if self._eigenvalue_element is None: + self._eigenvalue_element = ET.SubElement(self._settings_file, "eigenvalue") + + element = ET.SubElement(self._eigenvalue_element, "batches") + element.text = '{0}'.format(self._batches) + + + def createGenerationsPerBatchSubelement(self): + + if not self._generations_per_batch is None: + + if self._eigenvalue_element is None: + self._eigenvalue_element = ET.SubElement(self._settings_file, "eigenvalue") + + element = ET.SubElement(self._eigenvalue_element, "generations_per_batch") + element.text = '{0}'.format(self._generations_per_batch) + + + def createInactiveSubelement(self): + + if not self._inactive is None: + + if self._eigenvalue_element is None: + self._eigenvalue_element = ET.SubElement(self._settings_file, "eigenvalue") + + element = ET.SubElement(self._eigenvalue_element, "inactive") + element.text = '{0}'.format(self._inactive) + + + def createParticlesSubelement(self): + + if not self._particles is None: + + if self._eigenvalue_element is None: + self._eigenvalue_element = ET.SubElement(self._settings_file, "eigenvalue") + + element = ET.SubElement(self._eigenvalue_element, "particles") + element.text = '{0}'.format(self._particles) + + + def createSourceSubelement(self): + + self.createSourceSpaceSubelement() + self.createSourceEnergySubelement() + self.createSourceAngleSubelement() + + + def createSourceSpaceSubelement(self): + + + if not self._source_space_params is None: + + if self._source_subelement is None: + self._source_subelement = ET.SubElement(self._settings_file, "source") + + element = ET.SubElement(self._source_subelement, "space") + element.set("type", self._source_space_type) + + subelement = ET.SubElement(element, "parameters") + + text = '' + for param in self._source_space_params: + text += '{0} '.format(param) + subelement.text = text.rstrip(' ') + + + def createSourceAngleSubelement(self): + + if not self._source_angle_params is None: + + if self._source_subelement is None: + self._source_subelement = ET.SubElement(self._settings_file, "source") + + element = ET.SubElement(self._source_subelement, "angle") + element.set("type", self._source_angle_type) + + subelement = ET.SubElement(element, "parameters") + + text = '' + for param in self._source_angle_params: + text += '{0} '.format(param) + subelement.text = text.rstrip(' ') + + + def createSourceEnergySubelement(self): + + if not self._source_energy_params is None: + + if self._source_subelement is None: + self._source_subelement = ET.SubElement(self._settings_file, "source") + + element = ET.SubElement(self._source_subelement, "energy") + element.set("type", self._source_energy_type) + + subelement = ET.SubElement(element, "parameters") + + text = '' + for param in self._source_energy_params: + text += '{0} '.format(param) + subelement.text = text.rstrip(' ') + + + def createOutputSubelement(self): + + if not self._output is None: + element = ET.SubElement(self._settings_file, "output") + + for key in self._output.keys(): + subelement = ET.SubElement(element, key) + subelement.text = str(self._output[key]).lower() + + self.createOuputPathSubelement() + + + def createOuputPathSubelement(self): + + if not self._output_path is None: + element = ET.SubElement(self._settings_file, "output_path") + element.text = self._output_path + + + def createVerbositySubelement(self): + + if not self._verbosity is None: + element = ET.SubElement(self._settings_file, "verbosity") + element.text = '{0}'.format(self._verbosity) + + + def createStatepointSubelement(self): + + # Batches subelement + if not self._statepoint_batches is None: + element = ET.SubElement(self._settings_file, "state_point") + subelement = ET.SubElement(element, "batches") + text = '' + for batch in self._statepoint_batches: + text += '{0} '.format(batch) + subelement.text = text.rstrip(' ') + + # Interval subelement + elif not self._statepoint_interval is None: + element = ET.SubElement(self._settings_file, "state_point") + subelement = ET.SubElement(element, "interval") + subelement.text = '{0}'.format(self._statepoint_interval) + + + def createSourcepointSubelement(self): + + # Batches subelement + if not self._sourcepoint_batches is None: + element = ET.SubElement(self._settings_file, "source_point") + subelement = ET.SubElement(element, "batches") + text = '' + for batch in self._sourcepoint_batches: + text += '{0} '.format(batch) + subelement.text = text.rstrip(' ') + + # Interval subelement + elif not self._sourcepoint_interval is None: + element = ET.SubElement(self._settings_file, "source_point") + subelement = ET.SubElement(element, "interval") + subelement.text = '{0}'.format(self._sourcepoint_interval) + + # Separate subelement + if not self._sourcepoint_separate is None: + subelement = ET.SubElement(element, "separate") + subelement.text = '{0}'.format(str(self._sourcepoint_separate).lower()) + + # Write subelement + if not self._sourcepoint_write is None: + subelement = ET.SubElement(element, "write") + subelement.text = '{0}'.format(str(self._sourcepoint_write).lower()) + + # Overwrite latest subelement + if not self._sourcepoint_overwrite is None: + subelement = ET.SubElement(element, "overwrite_latest") + subelement.text = '{0}'.format(str(self._sourcepoint_overwrite).lower()) + + + def createConfidenceIntervalsSubelement(self): + + if not self._confidence_intervals is None: + element = ET.SubElement(self._settings_file, "confidence_intervals") + element.text = '{0}'.format(str(self._confidence_intervals).lower()) + + + def createCrossSectionsSubelement(self): + + if not self._cross_sections is None: + element = ET.SubElement(self._settings_file, "cross_sections") + element.text = '{0}'.format(self._cross_sections) + + + def createEnergyGridSubelement(self): + + if not self._energy_grid is None: + element = ET.SubElement(self._settings_file, "energy_grid") + element.text = '{0}'.format(self._energy_grid) + + + def createPTablesSubelement(self): + + if not self._ptables is None: + element = ET.SubElement(self._settings_file, "ptables") + element.text = '{0}'.format(str(self._ptables).lower()) + + + def createRunCMFDSubelement(self): + + if not self._run_cmfd is None: + element = ET.SubElement(self._settings_file, "run_cmfd") + element.text = '{0}'.format(str(self._run_cmfd).lower()) + + + def createSeedSubelement(self): + + if not self._seed is None: + element = ET.SubElement(self._settings_file, "seed") + element.text = '{0}'.format(self._seed) + + + def createSurvivalBiasingSubelement(self): + + if not self._survival_biasing is None: + element = ET.SubElement(self._settings_file, "survival_biasing") + element.text = '{0}'.format(str(self._survival_biasing).lower()) + + + def createCutoffSubelement(self): + + if not self._weight is None: + element = ET.SubElement(self._settings_file, "cutoff") + + subelement = ET.SubElement(element, "weight") + subelement.text = '{0}'.format(self._weight) + + subelement = ET.SubElement(element, "weight_avg") + subelement.text = '{0}'.format(self._weight_avg) + + + def createEntropySubelement(self): + + if not self._entropy_lower_left is None and \ + not self._entropy_upper_right is None: + + element = ET.SubElement(self._settings_file, "entropy") + + subelement = ET.SubElement(element, "dimension") + subelement.text = '{0} {1} {2}'.format(self._entropy_dimension[0], + self._entropy_dimension[1], + self._entropy_dimension[2]) + + subelement = ET.SubElement(element, "lower_left") + subelement.text = '{0} {1} {2}'.format(self._entropy_lower_left[0], + self._entropy_lower_left[1], + self._entropy_lower_left[2]) + + subelement = ET.SubElement(element, "upper_right") + subelement.text = '{0} {1} {2}'.format(self._entropy_upper_right[0], + self._entropy_upper_right[1], + self._entropy_upper_right[2]) + + + def createNoReduceSubelement(self): + + if not self._no_reduce is None: + element = ET.SubElement(self._settings_file, "no_reduce") + element.text = '{0}'.format(str(self._no_reduce).lower()) + + + def createThreadsSubelement(self): + + if not self._threads is None: + element = ET.SubElement(self._settings_file, "threads") + element.text = '{0}'.format(self._threads) + + + def createTraceSubelement(self): + + if not self._trace is None: + element = ET.SubElement(self._settings_file, "trace") + + text = '' + for item in self._trace: + text += '{0} '.format(item) + element.text = text.rstrip(' ') + + + def createTrackSubelement(self): + + if not self._track is None: + element = ET.SubElement(self._settings_file, "track") + + text = '' + for item in self._track: + text += '{0} '.format(item) + element.text = text.rstrip(' ') + + + def createUFSSubelement(self): + + if not self._ufs_lower_left is None and \ + not self._ufs_upper_right is None: + + element = ET.SubElement(self._settings_file, "uniform_fs") + + subelement = ET.SubElement(element, "dimension") + subelement.text = '{0}'.format(self._ufs_dimension) + + subelement = ET.SubElement(element, "lower_left") + subelement.text = '{0} {1} {2}'.format(self._ufs_lower_left[0], + self._ufs_lower_left[1], + self._ufs_lower_left[2]) + + subelement = ET.SubElement(element, "upper_right") + subelement.text = '{0} {1} {2}'.format(self._ufs_upper_right[0], + self._ufs_upper_right[1], + self._ufs_upper_right[2]) + + + def exportToXML(self): + + self.createEigenvalueSubelement() + self.createSourceSubelement() + self.createOutputSubelement() + self.createStatepointSubelement() + self.createSourcepointSubelement() + self.createConfidenceIntervalsSubelement() + self.createCrossSectionsSubelement() + self.createEnergyGridSubelement() + self.createPTablesSubelement() + self.createRunCMFDSubelement() + self.createSeedSubelement() + self.createSurvivalBiasingSubelement() + self.createCutoffSubelement() + self.createEntropySubelement() + self.createNoReduceSubelement() + self.createThreadsSubelement() + self.createVerbositySubelement() + self.createTraceSubelement() + self.createTrackSubelement() + self.createUFSSubelement() + + # Clean the indentation in the file to be user-readable + clean_xml_indentation(self._settings_file) + + # Write the XML Tree to the settings.xml file + tree = ET.ElementTree(self._settings_file) + tree.write("settings.xml", xml_declaration=True, + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py new file mode 100644 index 0000000000..d2e6467827 --- /dev/null +++ b/src/utils/openmc/statepoint.py @@ -0,0 +1,1196 @@ +#!/usr/bin/env python + +import struct +import numpy as np +import scipy.stats +import openmc +from openmc.constants import * + + + +class SourceSite(object): + + 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) + string += '{0: <16}{1}{2}\n'.format('\t(x,y,z)', '=\t', self._xyz) + string += '{0: <16}{1}{2}\n'.format('\t(u,v,w)', '=\t', self._uvw) + return string + + + +class StatePoint(object): + + def __init__(self, filename): + + openmc.reset_auto_ids() + + if filename.endswith('.h5'): + import h5py + self._f = h5py.File(filename, 'r') + self._hdf5 = True + else: + self._f = open(filename, 'rb') + self._hdf5 = False + + # Set flags for what data has been read + self._results = False + self._source = False + + # Read all metadata + self._read_metadata() + + # Read information needed to reconstruct the Materials and Geometry + self._read_geometry() + + # Read information about tally meshes + self._read_meshes() + + # Read tally metadata + self._read_tallies() + + + def _read_metadata(self): + + # Read filetype + self._filetype = self._get_int(path='metadata/filetype')[0] + + # Read statepoint revision + self._revision = self._get_int(path='metadata/revision')[0] + if self._revision != 12: + raise Exception('Statepoint Revision is not consistent.') + + # Read OpenMC version + if self._hdf5: + self._version = [self._get_int(path='metadata/version_major')[0], + self._get_int(path='metadata/version_minor')[0], + self._get_int(path='metadata/version_release')[0]] + else: + self._version = self._get_int(3) + + # Read date and time + self._date_and_time = self._get_string(19, path='metadata/date_and_time') + + # Read path + self._path = self._get_string(255, path='metadata/path').strip() + + # Read random number seed + self._seed = self._get_long(path='metadata/seed')[0] + + # Read run information + self._run_mode = self._get_int(path='metadata/run_mode')[0] + self._n_particles = self._get_long(path='metadata/n_particles')[0] + self._n_batches = self._get_int(path='metadata/n_batches')[0] + + # Read current batch + self._current_batch = self._get_int(path='metadata/current_batch')[0] + + # Read whether or not the source site distribution is present + self._source_present = self._get_int(path='metadata/source_present')[0] + + # Read criticality information + if self._run_mode == 2: + self._read_criticality() + + + def _read_criticality(self): + + # Read criticality information + if self._run_mode == 2: + + self._n_inactive = self._get_int(path='criticality/n_inactive')[0] + self._gen_per_batch = self._get_int(path='criticality/gen_per_batch')[0] + self._k_batch = self._get_double( + self._current_batch*self._gen_per_batch, + path='criticality/k_generation') + self._entropy = self._get_double( + self._current_batch*self._gen_per_batch, path='criticality/entropy') + + self._k_col_abs = self._get_double(path='criticality/k_col_abs')[0] + self._k_col_tra = self._get_double(path='criticality/k_col_tra')[0] + self._k_abs_tra = self._get_double(path='criticality/k_abs_tra')[0] + self._k_combined = self._get_double(2, path='criticality/k_combined') + + # Read CMFD information (if used) + self._read_cmfd() + + + def _read_cmfd(self): + + base = 'criticality/cmfd' + + # Read CMFD information + self._cmfd_on = self._get_int(path='{0}/cmfd_on'.format(base))[0] + + if self._cmfd_on == 1: + + self._cmfd_indices = self._get_int(4, path='{0}/indices'.format(base)) + self._k_cmfd = self._get_double(self._current_batch, + path='{0}/k_cmfd'.format(base)) + self._cmfd_src = self._get_double_array(np.product(self._cmfd_indices), + path='{0}/cmfd_src'.format(base)) + self._cmfd_src = np.reshape(self._cmfd_src, tuple(self._cmfd_indices), + order='F') + self._cmfd_entropy = self._get_double(self._current_batch, + path='{0}/cmfd_entropy'.format(base)) + self._cmfd_balance = self._get_double(self._current_batch, + path='{0}/cmfd_balance'.format(base)) + self._cmfd_dominance = self._get_double(self._current_batch, + path='{0}/cmfd_dominance'.format(base)) + self._cmfd_srccmp = self._get_double(self._current_batch, + path='{0}/cmfd_srccmp'.format(base)) + + + def _read_geometry(self): + + # Read in and initialize the Materials and Geometry + self._read_nuclides() + self._read_materials() + self._read_surfaces() + self._read_cells() + self._read_universes() + self._read_lattices() + self._finalize_geometry() + + + + def _read_nuclides(self): + + # Initialize dictionary for each Nuclide + # Keys - Nuclide IDs + # Values - Nuclide objects + self._nuclides = dict() + + self._n_nuclides = self._get_int(path='model/nuclides/n_nuclides')[0] + + # OpenMC Nuclide IDs (defined internally) + self._nuclide_ids = self._get_int( + self._n_nuclides, path='model/nuclides/ids') + + # User-defined Nuclide string names + self._nuclide_keys = dict() + for i in range(self._n_nuclides): + nuclide_name = self._get_string(11, + path='model/nuclides/keys/{0}'.format(i+1)) + + # Remove leading and trailing characters from string name + nuclide_name = nuclide_name.lstrip('[\'') + nuclide_name = nuclide_name.rstrip('\']') + + self._nuclide_keys[i+1] = nuclide_name + + # Iterate over and initalize all Nuclides + for i, nuclide_id in enumerate(self._nuclide_ids): + + # Read the Nuclide's name (e.g., 'H-1' or 'U-235') + name = self._nuclide_keys[i+1].split('.')[0] + + # Read the Nuclide's cross-section identifier (e.g., '70c') + xs = self._nuclide_keys[i+1].split('.')[1] + + # Initialize this Nuclide and add it to the global dictionary of Nuclides + self._nuclides[nuclide_id] = openmc.Nuclide(name=name, xs=xs) + + + + def _read_materials(self): + + # Initialize dictionary for each Material + # Keys - Material keys + # Values - Material objects + self._materials = dict() + + self._n_materials = self._get_int(path='model/materials/n_materials')[0] + + # OpenMC Material IDs (redefined internally from user definitions) + self._material_ids = self._get_int( + self._n_materials, path='model/materials/ids') + + # User-defined Material IDs + self._material_keys = self._get_int( + self._n_materials, path='model/materials/keys') + + # Build dictionary of Materials + base = 'model/materials/material ' + + # Iterate over all Materials + for material_key in self._material_keys: + + # Read the number of Nuclides + n_nuclides = self._get_int( + path='{0}{1}/n_nuclides'.format(base, material_key))[0] + + # Read the Nuclide IDs + nuclide_ids = self._get_int( + n_nuclides, path='{0}{1}/nuclide_ids'.format(base, material_key)) + + # Read the Material's density + density = self._get_double( + path='{0}{1}/density'.format(base, material_key))[0] + + # Read the atom/b-cm densities for each Nuclide + nuclide_densities = self._get_double( + n_nuclides, path='{0}{1}/nuclide_densities'.format(base, material_key)) + + # Read the number of nuclides with S(a,b) tables + n_sab = self._get_int(path='{0}{1}/n_sab'.format(base, material_key))[0] + + sab_names = list() + sab_xs = list() + + # Read the names of the S(a,b) tables for this Material + for i in range(1, n_sab+1): + + sab_table = self._get_string(12, + path='{0}{1}/sab_tables/{2}'.format(base, material_key, i)) + + # Read the cross-section identifiers for each S(a,b) table + sab_names.append(sab_table.split('.')[0]) + sab_xs.append(sab_table.split('.')[1]) + + # Create the Material + material = openmc.Material(material_id=material_key) + + # Set the Material's density to g/cm3 - this is what is used in OpenMC + material.setDensity(density=density, units='g/cm3') + + # Add all Nuclides to the Material + for i, nuclide_id in enumerate(nuclide_ids): + nuclide = self._nuclides[nuclide_id] + density = nuclide_densities[i] + material.addNuclide(nuclide, percent=density, percent_type='ao') + + # Add S(a,b) table(s?) to the Material + for i in range(n_sab): + name = sab_names[i] + xs = sab_xs[i] + material.setSAlphaBeta(name, xs) + + # Add the Material to the global dictionary of all Materials + self._materials[material_key] = material + + + def _read_surfaces(self): + + # Initialize dictionary for each Surface + # Keys - Surface keys + # Values - Surfacee objects + self._surfaces = dict() + + self._n_surfaces = self._get_int(path='model/surfaces/n_surfaces')[0] + + # OpenMC Surface IDs (redefined internally from user definitions) + self._surface_ids = self._get_int( + self._n_surfaces, path='model/surfaces/ids') + + # User-defined Surface IDs + self._surface_keys = self._get_int( + self._n_surfaces, path='model/surfaces/keys') + + #Build dictionary of Surfaces + base = 'model/surfaces/surface ' + + # Iterate over all Surfaces + for surface_key in self._surface_keys: + + # Read the Surface type and boundary condition type codes + type = self._get_int(path='{0}{1}/type'.format(base, surface_key))[0] + bc = self._get_int(path='{0}{1}/bc'.format(base, surface_key))[0] + + # Read the Surface coefficients + n_coeffs = self._get_int(path='{0}{1}/n_coeffs'.format(base, surface_key))[0] + coeffs = self._get_double( + n_coeffs, path='{0}{1}/coeffs'.format(base, surface_key)) + + # Create the Surface based on it's type + + if SURFACE_TYPES[type] == 'x-plane': + x0 = coeffs[0] + surface = openmc.XPlane(surface_key, BC_TYPES[bc], x0) + + elif SURFACE_TYPES[type] == 'y-plane': + y0 = coeffs[0] + surface = openmc.YPlane(surface_key, BC_TYPES[bc], y0) + + elif SURFACE_TYPES[type] == 'z-plane': + z0 = coeffs[0] + surface = openmc.ZPlane(surface_key, BC_TYPES[bc], z0) + + elif SURFACE_TYPES[type] == 'plane': + A = coeffs[0] + B = coeffs[1] + C = coeffs[2] + D = coeffs[3] + surface = openmc.Plane(surface_key, BC_TYPES[bc], A, B, C, D) + + elif SURFACE_TYPES[type] == 'x-cylinder': + y0 = coeffs[0] + z0 = coeffs[1] + R = coeffs[2] + surface = openmc.XCylinder(surface_key, BC_TYPES[bc], y0, z0, R) + + elif SURFACE_TYPES[type] == 'y-cylinder': + x0 = coeffs[0] + z0 = coeffs[1] + R = coeffs[2] + surface = openmc.YCylinder(surface_key, BC_TYPES[bc], x0, z0, R) + + elif SURFACE_TYPES[type] == 'z-cylinder': + x0 = coeffs[0] + y0 = coeffs[1] + R = coeffs[2] + surface = openmc.ZCylinder(surface_key, BC_TYPES[bc], x0, y0, R) + + elif SURFACE_TYPES[type] == 'sphere': + x0 = coeffs[0] + y0 = coeffs[1] + z0 = coeffs[2] + R = coeffs[3] + surface = openmc.Sphere(surface_key, BC_TYPES[bc], x0, y0, z0, R) + + elif SURFACE_TYPES[type] in ['x-cone', 'y-cone', 'z-cone']: + x0 = coeffs[0] + y0 = coeffs[1] + z0 = coeffs[2] + R2 = coeffs[3] + + if SURFACE_TYPES[type] == 'x-cone': + surface = openmc.XCone(surface_key, BC_TYPES[bc], x0, y0, z0, R2) + if SURFACE_TYPES[type] == 'y-cone': + surface = openmc.YCone(surface_key, BC_TYPES[bc], x0, y0, z0, R2) + if SURFACE_TYPES[type] == 'z-cone': + surface = openmc.ZCone(surface_key, BC_TYPES[bc], x0, y0, z0, R2) + + # Add Surface to global dictionary of all Surfaces + self._surfaces[surface_key] = surface + + + def _read_cells(self): + + # Initialize dictionary for each Cell + # Keys - Cell keys + # Values - Cell objects + self._cells = dict() + + self._n_cells = self._get_int(path='model/cells/n_cells')[0] + + # OpenMC Cell IDs (redefined internally from user definitions) + self._cell_ids = self._get_int(self._n_cells, path='model/cells/ids') + + # User-defined Cell IDs + self._cell_keys = self._get_int(self._n_cells, path='model/cells/keys') + + # Build dictionary of Cells + base = 'model/cells/cell ' + + # This is used later when the Universes and Lattices have been created + # Keys - Cell ID (user-defined) + # Values - (fill type string, fill ID) tuple + self._cell_fills = dict() + + # Iterate over all Cells + for cell_key in self._cell_keys: + + # Create this Cell + cell = openmc.Cell(cell_id=cell_key) + + # Read the integer fill type code + fill_type = self._get_int( + path='{0}{1}/fill_type'.format(base, cell_key))[0] + + if FILL_TYPES[fill_type] == 'normal': + fill_id = self._get_int( + path='{0}{1}/material'.format(base, cell_key))[0] + + elif FILL_TYPES[fill_type] == 'fill': + fill_id = self._get_int(path='{0}{1}/fill'.format(base, cell_key))[0] + maps = self._get_int(path='{0}{1}/maps'.format(base, cell_key))[0] + + if maps > 0: + cell.setOffset(self._get_int( + maps, path='{0}{1}/offset'.format(base, cell_key))) + + else: + fill_id = self._get_int(path='{0}{1}/lattice'.format(base, cell_key))[0] + + # Store the Cell fill information for use after Universe/Lattice creation + self._cell_fills[cell_key] = (FILL_TYPES[fill_type], fill_id) + + # Read the number of Surfaces and Surface IDs + n_surfaces = self._get_int( + path='{0}{1}/n_surfaces'.format(base, cell_key))[0] + + # If the Cell contains Surfaces, add them to the Cell + if n_surfaces > 0: + + surfaces = self._get_int( + n_surfaces, path='{0}{1}/surface_ids'.format(base, cell_key)) + + # Iterate over all Surfaces and add them to the Cell + for surface_key in surfaces: + + # Determine Surface halfspace + halfspace = np.sign(surface_key) + + surface_key = np.abs(surface_key) + + # Get the user-defined Sell ID and the corresponding Sell + surface_index = self._surface_ids.index(surface_key) + surface_id = self._surface_keys[surface_index] + surface = self._surfaces[surface_id] + + # Add this Surface, halsfpace to the Cell + cell.addSurface(surface, halfspace) + + # Add the Cell to the global dictionary of all Cells + self._cells[cell_key] = cell + + + + def _read_universes(self): + + # Initialize dictionary for each Universe + # Keys - Universe keys + # Values - Universe objects + self._universes = dict() + + self._n_universes = self._get_int(path='model/universes/n_universes')[0] + + # OpenMC Universe IDs (redefined internally from user definitions) + self._universe_ids = self._get_int( + self._n_universes, path='model/universes/ids') + + # User-defined Universe IDs + self._universe_keys = self._get_int( + self._n_universes, path='model/universes/keys') + + # Build dictionary of Universes + base = 'model/universes/universe ' + + # Iterate over all Universes + for universe_key in self._universe_keys: + + # Create this Universe + universe = openmc.Universe(universe_id=universe_key) + + # Read the number of Cells in this Universe + n_cells = self._get_int( + path='{0}{1}/n_cells'.format(base, universe_key))[0] + + # If the Universe contains Cells, add them to the Universe + if n_cells > 0: + + # Read Cell keys (internal OpenMC IDs) + cell_keys = self._get_int( + n_cells, path='{0}{1}/cell_ids'.format(base, universe_key)) + + # Add each Cell to the Universe + for cell_key in cell_keys: + + # Get the user-defined Cell ID and the corresponding Cell + cell_index = self._cell_ids.index(cell_key) + cell_id = self._cell_keys[cell_index] + cell = self._cells[cell_id] + + # Add the Cell to the Universe + universe.addCell(cell) + + # Add the Universe to the global list of Universes + self._universes[universe_key] = universe + + + def _read_lattices(self): + + # Initialize lattices for each Lattice + # Keys - Lattice keys + # Values - Lattice objects + self._lattices = dict() + + self._n_lattices = self._get_int(path='model/lattices/n_lattices')[0] + + # Read a list of the IDs for each object in the geometry + if self._n_lattices > 0: + + # OpenMC Lattice IDs (redefined internally from user definitions) + self._lattice_ids = self._get_int( + self._n_lattices, path='model/lattices/ids') + + # User-defined Lattice IDs + self._lattice_keys = self._get_int( + self._n_lattices, path='model/lattices/keys') + + else: + self._lattice_keys = dict() + self._lattice_ids = dict() + + # Build dictionary of Lattices + base = 'model/lattices/lattice ' + + # Iterate over all Lattices + for lattice_key in self._lattice_keys: + + # Read the type of Lattice (e.g., rectangular or hexagonal) + type = self._get_int(path='{0}{1}/type'.format(base, lattice_key))[0] + + # Read the Lattice dimensions + dimension = self._get_int( + 3, path='{0}{1}/dimension'.format(base, lattice_key)) + + # Read the Lattice cell widths + width = self._get_double( + 3, path='{0}{1}/width'.format(base, lattice_key)) + + # Read the lower left coordinates + lower_left = self._get_double( + 3, path='{0}{1}/lower_left'.format(base, lattice_key)) + + # Read the ID of the Material outside the lattice (-1 for void) + outside = self._get_int( + path='{0}{1}/outside'.format(base, lattice_key))[0] + + # Read the offset(s) + offset_size = self._get_int( + path='{0}{1}/offset_size'.format(base, lattice_key))[0] + + # Read the lattice maps + maps = self._get_int(path='{0}{1}/maps'.format(base, lattice_key))[0] + + if offset_size > 0: + + path = '{0}{1}/offsets'.format(base, lattice_key) + + if self._hdf5: + offsets = self._f[path].value + offsets = np.swapaxes(offsets, 0, 1) + offsets = np.swapaxes(offsets, 1, 2) + else: + offsets = np.array(self._get_int(offset_size, path=path)) + offsets.shape = (dimension[2], dimension[1], dimension[0], maps) + + # Read the Universes filling the Lattice + path = '{0}{1}/universes'.format(base, lattice_key) + + # Import the Universe IDs and convert to 3D array in (x,y,z) order + if self._hdf5: + universe_ids = self._f[path][...] + universe_ids = np.swapaxes(universe_ids, 0, 1) + universe_ids = np.swapaxes(universe_ids, 1, 2) + else: + num_values = dimension[0] * dimension[1] * dimension[2] + universe_ids = np.array(self._get_int(num_values, path=path)) + universe_ids.shape = dimension + + # Create the Lattice + lattice = openmc.Lattice(lattice_id=lattice_key, type=LATTICE_TYPES[type]) + lattice.setDimension(tuple(dimension)) + lattice.setLowerLeft(lower_left) + lattice.setWidth(width) + + # Build array of Universe pointers for the Lattice + universes = np.ndarray(tuple(universe_ids.shape), + dtype=openmc.Universe) + + 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._universes[universe_ids[x,y,z]] + + lattice.setUniverses(universes) + + # If the Material specified outside the Lattice is not void (-1) + if outside != -1: + lattice.setOutside(self._materials[outside]) + + if offset_size > 0: + lattice.setOffsets(offsets) + + # Add the Lattice to the global dictionary of all Lattices + self._lattices[lattice_key] = lattice + + + def _finalize_geometry(self): + + # Initialize Geometry object + self._geometry = openmc.Geometry() + + # Iterate over all Cells and add fill Materials, Universes and Lattices + for cell_key in self._cell_fills.keys(): + + # Determine the fill type ('normal', 'universe', or 'lattice') and ID + fill_type = self._cell_fills[cell_key][0] + fill_id = self._cell_fills[cell_key][1] + + # Retrieve the object corresponding to the fill type and ID + if fill_type == 'normal': + if fill_id > 0: + fill = self._materials[fill_id] + else: + fill = 'void' + elif fill_type == 'fill': + fill = self._universes[fill_id] + else: + fill = self._lattices[fill_id] + + # Set the fill for the Cell + self._cells[cell_key].setFill(fill) + + # Set the root universe for the Geometry + root_universe = self._universes[0] + self._geometry.setRootUniverse(root_universe) + + + def _read_meshes(self): + + # Initialize dictionaries for the Meshes + # Keys - Mesh IDs + # Values - Mesh objects + self._meshes = dict() + + # Read the number of Meshes + self._n_meshes = self._get_int(path='meshes/n_meshes')[0] + + # Read a list of the IDs for each Mesh + if self._n_meshes > 0: + + # OpenMC Mesh IDs (redefined internally from user definitions) + self._mesh_ids = self._get_int( + self._n_meshes, path='meshes/ids') + + # User-defined Mesh IDs + self._mesh_keys = self._get_int( + self._n_meshes, path='meshes/keys') + + else: + self._mesh_keys = list() + self._mesh_ids = list() + + # Build dictionary of Meshes + base = 'meshes/mesh ' + + # Iterate over all Meshes + for mesh_key in self._mesh_keys: + + # Read the user-specified Mesh ID and type + mesh_id = self._get_int(path='{0}{1}/id'.format(base, mesh_key))[0] + mesh_type = self._get_int(path='{0}{1}/type'.format(base, mesh_key))[0] + + # Get the Mesh dimension + n_dimension = self._get_int( + path='{0}{1}/n_dimension'.format(base, mesh_key))[0] + + # Read the mesh dimensions, lower-left coordinates, + # upper-right coordinates, and width of each mesh cell + dimension = self._get_int( + n_dimension, path='{0}{1}/dimension'.format(base, mesh_key)) + lower_left = self._get_double( + n_dimension, path='{0}{1}/lower_left'.format(base, mesh_key)) + upper_right = self._get_double( + n_dimension, path='{0}{1}/upper_right'.format(base, mesh_key)) + width = self._get_double( + n_dimension, path='{0}{1}/width'.format(base, mesh_key)) + + # Create the Mesh and assign properties to it + mesh = openmc.Mesh(mesh_id) + + mesh.setDimension(dimension) + mesh.setWidth(width) + mesh.setLowerLeft(lower_left) + mesh.setUpperRight(upper_right) + + #FIXME: Set the mesh type to 'rectangular' by default + mesh.setType('rectangular') + + # 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 + self._tallies = dict() + + # Read the number of tallies + self._n_tallies = self._get_int(path='/tallies/n_tallies')[0] + + # Read a list of the IDs for each Tally + if self._n_tallies > 0: + + # OpenMC Tally IDs (redefined internally from user definitions) + self._tally_ids = self._get_int( + self._n_tallies, path='tallies/ids') + + # User-defined Tally IDs + self._tally_keys = self._get_int( + self._n_tallies, path='tallies/keys') + + else: + self._tally_keys = list() + self._tally_ids = list() + + base = 'tallies/tally ' + + # Iterate over all Tallies + for tally_key in self._tally_keys: + + # Read user-specified Tally label (if specified) + label_size = self._get_int( + path='{0}{1}/label_size'.format(base, tally_key))[0] + + if label_size > 0: + label = self._get_string( + label_size, path='{0}{1}/label'.format(base, tally_key)) + + # Remove leading and trailing characters from string label + label = label.lstrip('[\'') + label = label.rstrip('\']') + + # Read integer Tally estimator type code (analog or tracklength) + estimator_type = self._get_int( + path='{0}{1}/estimator'.format(base, tally_key))[0] + + # Read the Tally size specifications + n_realizations = self._get_int( + path='{0}{1}/n_realizations'.format(base, tally_key))[0] + + # Create Tally object and assign basic properties + tally = openmc.Tally(tally_key, label) + tally.setEstimator(ESTIMATOR_TYPES[estimator_type]) + tally.setNumRealizations(n_realizations) + + # Read the number of Filters + n_filters = self._get_int( + path='{0}{1}/n_filters'.format(base, tally_key))[0] + + subbase = '{0}{1}/filter '.format(base, tally_key) + + # Initialize the stride + stride = 1 + + # Initialize all Filters + for j in range(1, n_filters+1): + + # Read the integer Filter type code + filter_type = self._get_int( + path='{0}{1}/type'.format(subbase, j))[0] + + # Read the Filter offset + offset = self._get_int(path='{0}{1}/offset'.format(subbase, j))[0] + + n_bins = self._get_int(path='{0}{1}/n_bins'.format(subbase, j))[0] + + if n_bins <= 0: + msg = 'Unable to create Filter {0} for Tally ID={2} since ' \ + 'no bins were specified'.format(j, tally_key) + raise ValueError(msg) + + # Read the bin values + if FILTER_TYPES[filter_type] in ['energy', 'energyout']: + bins = self._get_double( + n_bins+1, path='{0}{1}/bins'.format(subbase, j)) + + elif FILTER_TYPES[filter_type] == 'mesh': + bins = self._get_int(path='{0}{1}/bins'.format(subbase, j))[0] + + elif FILTER_TYPES[filter_type] == 'surface': + + bins = self._get_int(n_bins, path='{0}{1}/bins'.format(subbase, j)) + + # Get the user-defined Surface IDs from the internal OpenMC + # Surface IDs for the bin from the Surface key mapping + for k, bin in enumerate(bins): + surface_index = self._surface_ids.index(bin) + surface_id = self._surface_keys[surface_index] + bins[k] = surface_id + + elif FILTER_TYPES[filter_type] == 'distribcell': + + bins = self._get_int(path='{0}{1}/bins'.format(subbase, j))[0] + + # Get the user-defined Cell ID from the internal OpenMC + # Cell ID for the bin from the cell key mapping + cell_index = self._cell_ids.index(bins) + cell_id = self._cell_keys[cell_index] + bins = cell_id + + elif FILTER_TYPES[filter_type] in ['cell', 'cellborn']: + + bins = self._get_int(n_bins, path='{0}{1}/bins'.format(subbase, j)) + + # Get the user-defined Cell IDs from the internal OpenMC + # Cell IDs for the bin from the cell key mapping + for k, bin in enumerate(bins): + cell_index = self._cell_ids.index(bin) + cell_id = self._cell_keys[cell_index] + bins[k] = cell_id + + elif FILTER_TYPES[filter_type] == 'universe': + + bins = self._get_int( + n_bins, path='{0}{1}/bins'.format(subbase, j)) + + # Get the user-defined Universe IDs from the internal OpenMC + # Surface IDs for the bin from the Surface key mapping + for k, bin in enumerate(bins): + universe_index = self._universe_ids.index(bin) + universe_id = self._universe_keys[universe_index] + bins[k] = universe_id + + elif FILTER_TYPES[filter_type] == 'material': + + bins = self._get_int(n_bins, path='{0}{1}/bins'.format(subbase, j)) + + # Get the user-defined Material IDs from the internal OpenMC + # Surface IDs for the bin from the Surface key mapping + for k, bin in enumerate(bins): + material_index = self._material_ids.index(bin) + material_id = self._material_keys[material_index] + bins[k] = material_id + + # Create Filter object + filter = openmc.Filter(FILTER_TYPES[filter_type], bins) + filter.setOffset(offset) + filter.setStride(stride) + filter.setNumBins(n_bins) + + if FILTER_TYPES[filter_type] == 'mesh': + filter.setMesh(self._meshes[bins]) + + # Add Filter to the Tally + tally.addFilter(filter) + + # Update the stride for the next Filter + stride *= n_bins + + # Read Nuclide bins + n_nuclides = self._get_int( + path='{0}{1}/n_nuclides'.format(base, tally_key))[0] + + nuclide_ids = self._get_int( + n_nuclides, path='{0}{1}/nuclides'.format(base, tally_key)) + + # Add all Nuclides to the Tally + for nuclide_id in nuclide_ids: + if nuclide_id == -1: + tally.addNuclide(openmc.Nuclide('total')) + else: + tally.addNuclide(self._nuclides[nuclide_id]) + + # Read score bins + n_score_bins = self._get_int( + path='{0}{1}/n_score_bins'.format(base, tally_key))[0] + + tally.setNumScoreBins(n_score_bins) + + scores = [SCORE_TYPES[j] for j in self._get_int( + n_score_bins, path='{0}{1}/score_bins'.format(base, tally_key))] + + # Read the scattering moment order for all scores + scatt_order = self._get_int( + n_score_bins, path='{0}{1}/moment_order'.format(base, tally_key)) + + # Add the scores to the Tally + for j, score in enumerate(scores): + + # If this is a scattering moment, insert the scattering order + if 'scatter-n' in score or 'scatter-yn' in score or 'scatter-pn' in score: + score = score.replace('n', str(scatt_order[j])) + + tally.addScore(score) + + # Add Tally to the global dictionary of all Tallies + self._tallies[tally_key] = tally + + + + def read_results(self): + + # Number of realizations for global Tallies + self._n_realizations = self._get_int(path='n_realizations')[0] + + #TODO: Make global tallies use the general Tally class + # Read global Tallies + n_global_tallies = self._get_int(path='n_global_tallies')[0] + + if self._hdf5: + data = self._f['global_tallies'].value + self._global_tallies = np.column_stack((data['sum'], data['sum_sq'])) + + else: + self._global_tallies = np.array(self._get_double(2*n_global_tallies)) + self._global_tallies.shape = (n_global_tallies, 2) + + # Flag indicating if Tallies are present + self._tallies_present = self._get_int(path='tallies/tallies_present')[0] + + base = 'tallies/tally ' + + # Read Tally results + if self._tallies_present: + + # Iterate over and extract the results for all Tallies + for tally_key in self._tally_keys: + + # Get this Tally + tally = self._tallies[tally_key] + + # Compute the total number of bins for this Tally + num_tot_bins = tally.getNumBins() + + # Extract Tally data from the file + if self._hdf5: + data = self._f['{0}{1}/results'.format(base, tally_key)].value + sum = data['sum'] + sum_sq = data['sum_sq'] +# sum = np.column_stack(data['sum']) +# sum_sq = np.column_stack(data['sum_sq']) + + else: + results = np.array(self._get_double(2*num_tot_bins)) + sum = results[0::2] + sum_sq = results[1::2] + + # Define a routine to convert 0 to 1 + nonzero = lambda val: 1 if not val else val + + # Reshape the results arrays + new_shape = (nonzero(tally.getNumFilterBins()), + nonzero(tally.getNumNuclides()), + nonzero(tally.getNumScoreBins())) + + sum = np.reshape(sum, new_shape) + sum_sq = np.reshape(sum_sq, new_shape) + + # Set the data for this Tally + tally.setResults(sum=sum, sum_sq=sum_sq) + + # Indicate that Tally results have been read + self._results = True + + + def read_source(self): + + # Check whether Tally results have been read + if not self._results: + self.read_results() + + # Check if source bank is in statepoint + if not self._source_present: + print('Unable to read the source since it is not in the statepoint file') + return + + # Initialize a NumPy array for the source sites + self._source = np.empty(self._n_particles, dtype=SourceSite) + + # For HDF5 state points, copy entire bank + if self._hdf5: + source_sites = self._f['source_bank'].value + + # Initialize SourceSite object for each particle + for i in range(self._n_particles): + + # Initialize new source site + site = SourceSite() + + # Read position, angle, and energy + if self._hdf5: + site._weight, site._xyz, site._uvw, site._E = source_sites[i] + else: + site._weight = self._get_double()[0] + site._xyz = self._get_double(3) + site._uvw = self._get_double(3) + site._E = self._get_double()[0] + + # 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.""" + + # Determine significance level and percentile for two-sided CI + alpha = 1 - confidence + percentile = 1 - alpha/2 + + # Calculate t-value + 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 + for each Tally bin. + """ + + # Determine number of realizations + n = self._n_realizations + + # Calculate the standard deviation for each global tally + for i in range(len(self._global_tallies)): + + # Get sum and sum of squares + s, s2 = self._global_tallies[i] + + # Calculate sample mean and replace value + s /= n + self._global_tallies[i, 0] = s + + # Calculate standard deviation + if s != 0.0: + self._global_tallies[i, 1] = t_value * np.sqrt((s2 / n - s**2) / (n-1)) + + + # Calculate the sample mean and standard deviation for user-defined Tallies + for tally_id, tally in self._tallies.items(): + tally.computeStdDev(t_value) + + + def get_tally(self, score, filters, nuclides, + label='', estimator='tracklength'): + """Finds and returns a Tally object with certain properties. + + Parameters + ---------- + score : str + The score string + + filters : list + A list of Filter objects + + nuclides : list + A list of Nuclide objects + + label : str + The label specified for the Tally (default is '') + + estimator: str + The type of estimator ('tracklength' (default) or 'analog') + """ + + # Loop over the domain-to-tallies mapping to find the Tally + tally = None + + # Iterate over all tallies to find the appropriate one + for tally_id, test_tally in self._tallies.items(): + + # Determine if the queried Tally label is the same as this Tally + if not label == test_tally._label: + continue + + # Determine if the queried Tally estimator is the same as this Tally + if not estimator == test_tally._estimator: + continue + + # Determine if the queried Tally scores are the same as this Tally + if not score in test_tally._scores: + continue + + # Determine if the queried Tally filters is the same length as this Tally + if len(filters) != len(test_tally._filters): + continue + + # Determine if the queried Tally filters are the same as this Tally + contains_filters = True + + # Iterate over the filters requested by the user + for filter in filters: + if not filter in test_tally._filters: + contains_filters = False + break + + # Determine if the queried Nuclide is in this Tally + contains_nuclides = True + + # Iterate over the Nuclides requested by the user + for nuclide in nuclides: + if not nuclide in test_tally._nuclides: + contains_nuclides = False + break + + # If the Tally contained all Filters and Nuclides, return the Tally + if contains_filters and contains_nuclides: + tally = test_tally + break + + # If we did not find the Tally, return an error message + if tally is None: + raise LookupError('Unable to get Tally') + + return tally + + + def get_tally_id(self, score, filters, label='', estimator='tracklength'): + """Retrieve the Tally ID for a given list of filters and score(s). + + Parameters + ---------- + score : str + The score string + + filters : list + A list of Filter objects + + label : str + The label specified for the Tally (default is '') + + estimator: str + The type of estimator ('tracklength' (default) or 'analog') + """ + + tally = self.get_tally(score, filters, label, estimator) + return tally._id + + + + def _get_data(self, n, typeCode, size): + return list(struct.unpack('={0}{1}'.format(n,typeCode), + self._f.read(n*size))) + + def _get_int(self, n=1, path=None): + if self._hdf5: + return [int(v) for v in self._f[path].value] + else: + return [int(v) for v in self._get_data(n, 'i', 4)] + + def _get_long(self, n=1, path=None): + if self._hdf5: + return [long(v) for v in self._f[path].value] + else: + return [long(v) for v in self._get_data(n, 'q', 8)] + + def _get_float(self, n=1, path=None): + if self._hdf5: + return [float(v) for v in self._f[path].value] + else: + return [float(v) for v in self._get_data(n, 'f', 4)] + + def _get_double(self, n=1, path=None): + if self._hdf5: + return [float(v) for v in self._f[path].value] + else: + return [float(v) for v in self._get_data(n, 'd', 8)] + + def _get_double_array(self, n=1, path=None): + if self._hdf5: + return self._f[path].value + else: + return self._get_data(n, 'd', 8) + + def _get_string(self, n=1, path=None): + if self._hdf5: + return str(self._f[path].value) + else: + return str(self._get_data(n, 's', 1)[0]) diff --git a/src/utils/openmc/surface.py b/src/utils/openmc/surface.py new file mode 100644 index 0000000000..898658002b --- /dev/null +++ b/src/utils/openmc/surface.py @@ -0,0 +1,593 @@ +#!/usr/bin/env python + +from openmc.checkvalue import * +from openmc.constants import BC_TYPES +from xml.etree import ElementTree as ET + + +# A static variable for auto-generated Surface IDs +AUTO_SURFACE_ID = 10000 + +def reset_auto_surface_id(): + global AUTO_SURFACE_ID + AUTO_SURFACE_ID = 10000 + + + +class Surface(object): + + def __init__(self, surface_id=None, bc_type='transmission', name=''): + + # Initialize class attributes + self._id = None + self._name = '' + self._type = '' + self._bc_type = '' + + # A dictionary of the quadratic surface coefficients + # Key - coefficeint name + # Value - coefficient value + self._coeffs = dict() + + # An ordered list of the coefficient names to export to XML in the + # proper order + self._coeff_keys = list() + + self.setId(surface_id) + self.setBoundaryType(bc_type) + self.setName(name) + + + def setId(self, surface_id=None): + + if surface_id is None: + global AUTO_SURFACE_ID + self._id = AUTO_SURFACE_ID + AUTO_SURFACE_ID += 1 + + # Check that the ID is an integer and wasn't already used + elif not is_integer(surface_id): + msg = 'Unable to set a non-integer Surface ID {0}'.format(surface_id) + raise ValueError(msg) + + elif surface_id < 0: + msg = 'Unable to set Surface ID to {0} since it must be a ' \ + 'non-negative integer'.format(surface_id) + raise ValueError(msg) + + else: + self._id = surface_id + + + def setName(self, name): + + if not is_string(name): + msg = 'Unable to set name for Surface ID={0} with a non-string ' \ + 'value {1}'.format(self._id, name) + raise ValueError(msg) + + else: + self._name = name + + + def setBoundaryType(self, bc_type): + + if not is_string(bc_type): + msg = 'Unable to set boundary type for Surface ID={0} with a ' \ + 'non-string value {1}'.format(self._id, bc_type) + raise ValueError(msg) + + elif not bc_type in BC_TYPES.values(): + msg = 'Unable to set boundary type for Surface ID={0} to {1} which ' \ + 'is not trasmission, vacuum or reflective'.format(bc_type) + raise ValueError(msg) + + else: + self._bc_type = bc_type + + + def __repr__(self): + + string = 'Surface\n' + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t{0}', self._id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) + string += '{0: <16}{1}{2}\n'.format('\tBoundary', '=\t', self._bc_type) + + coeffs = '{0: <16}'.format('\tCoefficients') + '\n' + + for coeff in self._coeffs.keys(): + coeffs += '{0: <16}{1}{2}\n'.format(coeff, '=\t', self._coeffs[coeff]) + + string += coeffs + + return string + + + def createXMLSubElement(self): + + element = ET.Element("surface") + element.set("id", str(self._id)) + element.set("type", self._type) + element.set("boundary", self._bc_type) + + coeffs = '' + + for coeff in self._coeff_keys: + coeffs += '{0} '.format(self._coeffs[coeff]) + + element.set("coeffs", coeffs.rstrip(' ')) + + return element + + + +class Plane(Surface): + + def __init__(self, surface_id=None, bc_type='transmission', + A=None, B=None, C=None, D=None, name='',): + + # Initialize Plane class attributes + super(Plane, self).__init__(surface_id, bc_type, name=name) + + self._A = None + self._B = None + self._C = None + self._D = None + self._type = 'plane' + self._coeff_keys = ['A', 'B', 'C', 'D'] + + if not A is None: + self.setA(A) + + if not B is None: + self.setB(B) + + if not C is None: + self.setC(C) + + if not D is None: + self.setD(D) + + + def setA(self, A): + + if not is_integer(A) and not is_float(A): + msg = 'Unable to set A coefficient for Plane ID={0} to a non-integer ' \ + 'value {1}'.format(self._id, A) + raise ValueError(msg) + + self._coeffs['A'] = A + + + def setB(self, B): + + if not is_integer(B) and not is_float(B): + msg = 'Unable to set B coefficient for Plane ID={0} to a non-integer ' \ + 'value {1}'.format(self._id, B) + raise ValueError(msg) + + self._coeffs['B'] = B + + + def setC(self, C): + + if not is_integer(C) and not is_float(C): + msg = 'Unable to set C coefficient for Plane ID={0} to a non-integer ' \ + 'value {1}'.format(self._id, C) + raise ValueError(msg) + + self._coeffs['C'] = C + + + def setD(self, D): + + if not is_integer(D) and not is_float(D): + msg = 'Unable to set D coefficient for Plane ID={0} to a non-integer ' \ + 'value {1}'.format(self._id, D) + raise ValueError(msg) + + self._coeffs['D'] = D + + + +class XPlane(Plane): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, name=''): + + # Initialize XPlane class attributes + super(XPlane, self).__init__(surface_id, bc_type, name=name) + + self._x0 = None + self._type = 'x-plane' + self._coeff_keys = ['x0'] + + if not x0 is None: + self.setX0(x0) + + + def setX0(self, x0): + + if not is_integer(x0) and not is_float(x0): + msg = 'Unable to set x0 coefficient for XPlane ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, x0) + raise ValueError(msg) + + self._coeffs['x0'] = x0 + + + +class YPlane(Plane): + + def __init__(self, surface_id=None, bc_type='transmission', + y0=None, name=''): + + # Initialize YPlane class attributes + super(YPlane, self).__init__(surface_id, bc_type, name=name) + + self._y0 = None + self._type = 'y-plane' + self._coeff_keys = ['y0'] + + if not y0 is None: + self.setY0(y0) + + + def setY0(self, y0): + + if not is_integer(y0) and not is_float(y0): + msg = 'Unable to set y0 coefficient for XPlane ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, y0) + raise ValueError(msg) + + self._coeffs['y0'] = y0 + + + +class ZPlane(Plane): + + def __init__(self, surface_id=None, bc_type='transmission', + z0=None, name=''): + + # Initialize ZPlane class attributes + super(ZPlane, self).__init__(surface_id, bc_type, name=name) + + self._z0 = None + self._type = 'z-plane' + self._coeff_keys = ['z0'] + + if not z0 is None: + self.setZ0(z0) + + + def setZ0(self, z0): + + if not is_integer(z0) and not is_float(z0): + msg = 'Unable to set z0 coefficient for ZPlane ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, z0) + raise ValueError(msg) + + self._coeffs['z0'] = z0 + + + +class Cylinder(Surface): + + def __init__(self, surface_id=None, bc_type='transmission', + R=None, name=''): + + # Initialize Cylinder class attributes + super(Cylinder, self).__init__(surface_id, bc_type, name=name) + + self._R = None + self._coeff_keys = ['R'] + + if not R is None: + self.setR(R) + + + def setR(self, R): + + if not is_integer(R) and not is_float(R): + msg = 'Unable to set R coefficient for Cylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, R) + raise ValueError(msg) + + self._coeffs['R'] = R + + + +class XCylinder(Cylinder): + + def __init__(self, surface_id=None, bc_type='transmission', + y0=None, z0=None, R=None, name=''): + + # Initialize XCylinder class attributes + super(XCylinder, self).__init__(surface_id, bc_type, R, name=name) + + self._type = 'x-cylinder' + self._y0 = None + self._z0 = None + self._coeff_keys = ['y0', 'z0', 'R'] + + if not y0 is None: + self.setY0(y0) + + if not z0 is None: + self.setZ0(z0) + + + def setY0(self, y0): + + if not is_integer(y0) and not is_float(y0): + msg = 'Unable to set y0 coefficient for XCylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, y0) + raise ValueError(msg) + + self._coeffs['y0'] = y0 + + + def setZ0(self, z0): + + if not is_integer(z0) and not is_float(z0): + msg = 'Unable to set z0 coefficient for XCylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, z0) + raise ValueError(msg) + + self._coeffs['z0'] = z0 + + + +class YCylinder(Cylinder): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, z0=None, R=None, name=''): + + # Initialize YCylinder class attributes + super(YCylinder, self).__init__(surface_id, bc_type, R, name=name) + + self._type = 'y-cylinder' + self._x0 = None + self._z0 = None + self._coeff_keys = ['x0', 'z0', 'R'] + + if not x0 is None: + self.setX0(x0) + + if not z0 is None: + self.setZ0(z0) + + + def setX0(self, x0): + + if not is_integer(x0) and not is_float(x0): + msg = 'Unable to set x0 coefficient for YCylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, x0) + raise ValueError(msg) + + self._coeffs['x0'] = x0 + + + def setZ0(self, z0): + + if not is_integer(z0) and not is_float(z0): + msg = 'Unable to set z0 coefficient for YCylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, z0) + raise ValueError(msg) + + self._coeffs['z0'] = z0 + + +class ZCylinder(Cylinder): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, y0=None, R=None, name=''): + + # Initialize ZCylinder class attributes + # Initialize YPlane class attributes + super(ZCylinder, self).__init__(surface_id, bc_type, R, name=name) + + self._type = 'z-cylinder' + self._x0 = None + self._y0 = None + self._coeff_keys = ['x0', 'y0', 'R'] + + if not x0 is None: + self.setX0(x0) + + if not y0 is None: + self.setY0(y0) + + + def setX0(self, x0): + + if not is_integer(x0) and not is_float(x0): + msg = 'Unable to set x0 coefficient for ZCylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, x0) + raise ValueError(msg) + + self._coeffs['x0'] = x0 + + + def setY0(self, y0): + + if not is_integer(y0) and not is_float(y0): + msg = 'Unable to set y0 coefficient for ZCylinder ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, y0) + raise ValueError(msg) + + self._coeffs['y0'] = y0 + + + +class Sphere(Surface): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, y0=None, z0=None, R=None, name=''): + + # Initialize Sphere class attributes + super(Sphere, self).__init__(surface_id, bc_type, name=name) + + self._type = 'sphere' + self._x0 = None + self._y0 = None + self._z0 = None + self._R = None + self._coeff_keys = ['x0', 'y0', 'z0', 'R'] + + if not x0 is None: + self.setX0(x0) + + if not y0 is None: + self.setY0(y0) + + if not z0 is None: + self.setZ0(z0) + + if not R is None: + self.setZ0(R) + + + def setX0(self, x0): + + if not is_integer(x0) and not is_float(x0): + msg = 'Unable to set x0 coefficient for Sphere ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, x0) + raise ValueError(msg) + + self._coeffs['x0'] = x0 + + + def setY0(self, y0): + + if not is_integer(y0) and not is_float(y0): + msg = 'Unable to set y0 coefficient for Sphere ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, y0) + raise ValueError(msg) + + self._coeffs['y0'] = y0 + + + def setZ0(self, z0): + + if not is_integer(z0) and not is_float(z0): + msg = 'Unable to set z0 coefficient for Sphere ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, z0) + raise ValueError(msg) + + self._coeffs['z0'] = z0 + + + def setR(self, R): + + if not is_integer(R) and not is_float(R): + msg = 'Unable to set R coefficient for Sphere ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, R) + raise ValueError(msg) + + self._coeffs['R'] = R + + + +class Cone(Surface): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, y0=None, z0=None, R2=None, name=''): + + # Initialize Cone class attributes + super(Cone, self).__init__(surface_id, bc_type, name=name) + + self._x0 = None + self._y0 = None + self._z0 = None + self._R2 = None + self._coeff_keys = ['x0', 'y0', 'z0', 'R2'] + + if not x0 is None: + self.setX0(x0) + + if not y0 is None: + self.setY0(y0) + + if not z0 is None: + self.setZ0(z0) + + if not R2 is None: + self.setZ0(R2) + + + def setX0(self, x0): + + if not is_integer(x0) and not is_float(x0): + msg = 'Unable to set x0 coefficient for Cone ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, x0) + raise ValueError(msg) + + self._coeffs['x0'] = x0 + + + def setY0(self, y0): + + if not is_integer(y0) and not is_float(y0): + msg = 'Unable to set y0 coefficient for Cone ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, y0) + raise ValueError(msg) + + self._coeffs['y0'] = y0 + + + def setZ0(self, z0): + + if not is_integer(z0) and not is_float(z0): + msg = 'Unable to set z0 coefficient for Cone ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, z0) + raise ValueError(msg) + + self._coeffs['z0'] = z0 + + + def setR2(self, R2): + + if not is_integer(R2) and not is_float(R2): + msg = 'Unable to set R^2 coefficient for Cone ID={0} to a ' \ + 'non-integer value {1}'.format(self._id, R2) + raise ValueError(msg) + + self._coeffs['R2'] = R2 + + + +class XCone(Cone): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, y0=None, z0=None, R2=None, name=''): + + # Initialize XCone class attributes + super(XCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) + + self._type = 'x-cone' + + + +class YCone(Cone): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, y0=None, z0=None, R2=None, name=''): + + # Initialize YCone class attributes + super(YCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) + + self._type = 'y-cone' + + + +class ZCone(Cone): + + def __init__(self, surface_id=None, bc_type='transmission', + x0=None, y0=None, z0=None, R2=None, name=''): + + # Initialize ZCone class attributes + super(ZCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) + + self._type = 'z-cone' diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py new file mode 100644 index 0000000000..ae88d25819 --- /dev/null +++ b/src/utils/openmc/tallies.py @@ -0,0 +1,1147 @@ +#!/usr/bin/env python + +from openmc import Nuclide +from openmc.clean_xml import * +from openmc.checkvalue import * +from openmc.constants import * +from xml.etree import ElementTree as ET +import numpy as np +import os + + +# "Static" variables for auto-generated Tally and Mesh IDs +AUTO_TALLY_ID = 10000 +AUTO_MESH_ID = 10000 + + +def reset_auto_tally_id(): + global AUTO_TALLY_ID + AUTO_TALLY_ID = 10000 + + +def reset_auto_mesh_id(): + global AUTO_MESH_ID + AUTO_MESH_ID = 10000 + + +class Filter(object): + + # Initialize Filter class attributes + def __init__(self, type=None, bins=None): + + self._type = None + self._bins = None + self._mesh = None + self._offset = -1 + self._stride = None + + # FIXME + self._num_bins = None + + if not type is None: + self.setType(type) + + if not bins is None: + self.setBinEdges(bins) + + + def __eq__(self, filter2): + + # Check type + if self._type != filter2._type: + return False + + # Check bin edges + elif list(self._bins) != list(filter2._bins): + return False + + else: + return True + + + def __hash__(self): + hashable = list() + hashable.append(self._type) + hashable.append(self._bins) + return hash(tuple(hashable)) + + + def setType(self, type): + + if not type 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 + + + def setBinEdges(self, bins): + + if self._type is None: + msg = 'Unable to set bin edges 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)): + bins = [bins] + + # If the bins are in a collection, convert it to a list + else: + bins = list(bins) + + if self._type in ['cell', 'cellborn', 'surface', 'material', + 'universe', 'distribcell']: + + for edge in bins: + + if not is_integer(edge): + msg = 'Unable to add bin edge {0} to a {1} Filter since it ' \ + 'is a non-integer'.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 integer'.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): + msg = 'Unable to add bin edge {0} to {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 {1} Filter since it is ' \ + 'a negative value'.format(edge, self._type) + raise ValueError(msg) + + # 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 {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 bin edges {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]): + msg = 'Unable to add bin edge {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 edge {0} to mesh Filter since it ' \ + 'is a negative integer'.format(bins[0]) + raise ValueError(msg) + + # If all error checks passed, add bin edges + self._bins = bins + + + # FIXME + def setNumBins(self, num_bins): + + if not is_integer(num_bins) 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) + raise ValueError(msg) + + self._num_bins = num_bins + + + def setMesh(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) + + self._mesh = mesh + self.setType('mesh') + self.setBinEdges(self._mesh._id) + + + def setOffset(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) + + self._offset = offset + + + def setStride(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) + + if stride < 0: + msg = 'Unable to set stride {0} for a {1} Filter since it is a ' \ + 'negative value'.format(stride, self._type) + raise ValueError(msg) + + self._stride = stride + + + def getNumBins(self): + + # FIXME + #if self._type == 'mesh': + # num_bins = self._mesh.getNumMeshCells() + #elif self._type == 'energy' or self._type == 'energyout': + # num_bins = len(self._bins) - 1 + #else: + # num_bins = len(self._bins) + + #return num_bins + + return self._num_bins + + + def getBinIndex(self, bin): + + try: + index = self._bins.index(bin) + + except ValueError: + msg = 'Unable to get the bin index for Filter since {0} ' \ + 'is not one of the bins'.format(bin) + raise ValueError(msg) + + return 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) + string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offset) + return string + + + +class Mesh(object): + + def __init__(self, mesh_id=None, name=''): + + # Initialize Mesh class attributes + self._id = None + self._name = '' + self._type = 'rectangular' + self._dimension = None + self._lower_left = None + self._upper_right = None + self._width = None + + self.setId(mesh_id) + self.setName(name) + + + def setId(self, mesh_id=None): + + 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: + self._id = mesh_id + + + def setName(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 + + + def setType(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 setDimension(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) + + self._dimension = dimension + + + def setLowerLeft(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) + + self._lower_left = lower_left + + + def setUpperRight(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) + + self._upper_right = upper_right + + + def setWidth(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) + + 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) + string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) + string += '{0: <16}{1}{2}\n'.format('\tBasis', '=\t', self._dimension) + string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._lower_left) + string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._upper_right) + string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._width) + return string + + + def getNumMeshCells(self): + return np.prod(self._dimension) + + + def getMeshXML(self): + + 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]) + + 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]) + + 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 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]) + + return element + + + +class Tally(object): + + def __init__(self, tally_id=None, label=''): + + # Initialize Tally class attributes + self._id = None + self._label = None + self._filters = list() + self._nuclides = list() + self._scores = list() + self._estimator = None + + self._num_score_bins = 0 + self._num_realizations = 0 + + self._sum = None + self._sum_sq = None + self._mean = None + self._std_dev = None + + self.setId(tally_id) + self.setLabel(label) + + + def __eq__(self, tally2): + + # Check all filters + for filter in self._filters: + if not filter in tally2._filters: + return False + + # Check all nuclides + for nuclide in self._nuclides: + if not nuclide in tally2._nuclides: + return False + + # Check all scores + for score in self._scores: + if not score in tally2._scores: + return False + + if self._estimator != tally2._estimator: + return False + + return True + + + def __hash__(self): + hashable = list() + + for filter in self._filters: + hashable.append((filter._type, tuple(filter._bins))) + + for nuclide in self._nuclides: + hashable.append(nuclide._name) + + for score in self._scores: + hashable.append(score) + + hashable.append(self._estimator) + hashable.append(self._label) + + return hash(tuple(hashable)) + + + def setEstimator(self, estimator): + if not estimator in ['analog', 'tracklength']: + msg = 'Unable to set the estimator for Tally ID={0} to {1} since it ' \ + 'is not a valid estimator type'.format(self._id, estimator) + raise ValueError(msg) + + self._estimator = estimator + + + def __add__(self, other): + + # FIXME: Error checking: must check that results has been + # set and that # bins is the same + + new_tally = Tally() + new_tally._mean = self._mean + other._mean + new_tally._std_dev = np.sqrt(self._std_dev**2 + other._std_dev**2) + + + def setId(self, tally_id=None): + + if tally_id is None: + global AUTO_TALLY_ID + self._id = AUTO_TALLY_ID + AUTO_TALLY_ID += 1 + + # Check that the ID is an integer and wasn't already used + elif not is_integer(tally_id): + msg = 'Unable to set a non-integer Tally ID {0}'.format(tally_id) + raise ValueError(msg) + + elif tally_id < 0: + msg = 'Unable to set Tally ID to {0} since it must be a ' \ + 'non-negative integer'.format(tally_id) + raise ValueError(msg) + + else: + self._id = tally_id + + + def setLabel(self, label=None): + + if not is_string(label): + msg = 'Unable to set name for Tally ID={0} with a non-string ' \ + 'value {1}'.format(self._id, label) + raise ValueError(msg) + + else: + self._label = label + + + def addFilter(self, filter): + + global filters + + if not isinstance(filter, Filter): + msg = 'Unable to add Filter {0} to Tally ID={1} since it is not ' \ + 'a Filter object'.format(filter, self._id) + raise ValueError(msg) + + self._filters.append(filter) + + + def addNuclide(self, nuclide): + self._nuclides.append(nuclide) + + + def addScore(self, score): + + if not is_string(score): + msg = 'Unable to add score {0} to Tally ID={1} since it is not a ' \ + 'string'.format(score, self._id) + raise ValueError(msg) + + elif 'scatter-' in score: + + moment = score.split('-')[-1] + + if 'p' in moment.lower() or 'y' in moment.lower(): + moment = moment[1:] + + if int(moment) < 0 or int(moment) > 10: + msg = 'Unable to add score {0} to Tally ID={1} since OpenMC ' \ + 'can only tally the scattering moments between 0 and ' \ + '10'.format(score, self._id) + raise ValueError(msg) + + elif not score in SCORE_TYPES.values(): + msg = 'Unable to add score {0} to Tally ID={1} since it is not a ' \ + 'supported score in OpenMC'.format(score, self._id) + raise ValueError(msg) + + # If the score is already in the Tally, don't add it again + if score in self._scores: + return + else: + self._scores.append(score) + + + def setNumScoreBins(self, num_score_bins): + self._num_score_bins = num_score_bins + + + def setNumRealizations(self, num_realizations): + + if not is_integer(num_realizations): + msg = 'Unable to set the number of realizations to {0} for Tally ID={1} ' \ + 'since it is not an integer'.format(num_realizations) + raise ValueError(msg) + + elif num_realizations < 0: + msg = 'Unable to set the number of realizations to {0} for Tally ID={1} ' \ + 'since it is a negative value'.format(num_realizations) + raise ValueError(msg) + + self._num_realizations = num_realizations + + + def setResults(self, sum, sum_sq): + + if not isinstance(sum, (tuple, list, np.ndarray)): + msg = 'Unable to set the sum to {0}for Tally ID={1} since it is not ' \ + 'a Python tuple/list or NumPy array'.format(sum, self._id) + raise ValueError(msg) + + if not isinstance(sum_sq, (tuple, list, np.ndarray)): + msg = 'Unable to set the sum to {0}for Tally ID={1} since it is not ' \ + 'a Python tuple/list or NumPy array'.format(sum_sq, self._id) + raise ValueError(msg) + + self._sum = sum + self._sum_sq = sum_sq + + + def removeScore(self, score): + + if not score in self._scores: + msg = 'Unable to remove score {0} from Tally ID={1} since the Tally ' \ + 'does not contain this score'.format(score, self._id) + ValueError(msg) + + self._scores.remove(score) + + + def removeFilter(self, filter): + + if not filter in self._filters: + msg = 'Unable to remove filter {0} from Tally ID={1} since the Tally ' \ + 'does not contain this filter'.format(filter, self._id) + ValueError(msg) + + self._filters.remove(filter) + + + def removeNuclide(self, nuclide): + + if not nuclide in self._nuclides: + msg = 'Unable to remove nuclide {0} from Tally ID={1} since the Tally ' \ + 'does not contain this nuclide'.format(nuclide, self._id) + ValueError(msg) + + self._nuclides.remove(nuclide) + + + def computeStdDev(self, t_value=1.0): + + # Calculate sample mean and standard deviation + self._mean = self._sum / self._num_realizations + self._std_dev = t_value * np.sqrt((self._sum_sq / self._num_realizations - \ + self._mean**2) / (self._num_realizations -1)) + + + def __repr__(self): + + string = 'Tally\n' + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._label) + + string += '{0: <16}\n'.format('\tFilters') + + for filter in self._filters: + string += '{0: <16}\t\t{1}\t{2}\n'.format('', filter._type, + filter._bins) + + string += '{0: <16}{1}'.format('\tNuclides', '=\t') + + for nuclide in self._nuclides: + string += '{0} '.format(nuclide._name) + + string += '\n' + + string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self._scores) + string += '{0: <16}{1}{2}\n'.format('\tEstimator', '=\t', self._estimator) + + return string + + + def getNumFilters(self): + return len(self._filters) + + + def getNumFilterBins(self): + + num_bins = 1 + + for filter in self._filters: + num_bins *= filter.getNumBins() + + return num_bins + + + def getNumNuclides(self): + return len(self._nuclides) + + + def getNumScores(self): + return len(self._scores) + + + def getNumScoreBins(self): + return self._num_score_bins + + + def getNumBins(self): + num_bins = self.getNumFilterBins() + num_bins *= self.getNumNuclides() + num_bins *= self.getNumScoreBins() + return num_bins + + + def getTallyXML(self): + + element = ET.Element("tally") + + # Tally ID + element.set("id", str(self._id)) + + # Optional Tally label + if self._label != '': + element.set("label", self._label) + + # Optional Tally filters + for filter in self._filters: + + subelement = ET.SubElement(element, "filter") + subelement.set("type", str(filter._type)) + + if not filter._bins is None: + + bins = '' + for bin in filter._bins: + bins += '{0} '.format(bin) + + subelement.set("bins", bins.rstrip(' ')) + + # Optional Nuclides + if len(self._nuclides) > 0: + + nuclides = '' + for nuclide in self._nuclides: + nuclides += '{0} '.format(nuclide._name) + + subelement = ET.SubElement(element, "nuclides") + subelement.text = nuclides.rstrip(' ') + + # Scores + if len(self._scores) == 0: + msg = 'Unable to get XML for Tally ID={0} since it does not contain ' \ + 'any scores'.format(self._id) + raise ValueError(msg) + + else: + + scores = '' + for score in self._scores: + scores += '{0} '.format(score) + + subelement = ET.SubElement(element, "scores") + subelement.text = scores.rstrip(' ') + + # Tally estimator type + if not self._estimator is None: + subelement = ET.SubElement(element, "estimator") + subelement.text = self._estimator + + return element + + + def findFilter(self, filter_type, bins): + + filter = None + + for test_filter in self._filters: + + # Determine if the Filter has the same type as the one requested + if test_filter._type != filter_type: + continue + + # Determine if the Filter has the same bin edges as the one requested + elif test_filter._bins != bins: + continue + + else: + filter = test_filter + break + + # If we found the Filter, return it + if not filter is None: + return filter + + # Otherwise, throw an Exception + else: + msg = 'Unable to find filter type {0} with bin edges {1} in Tally ' \ + 'ID={2}'.format(filter_type, bins, self._id) + raise ValueError(msg) + + + def getScoreIndex(self, score): + + try: + index = self._scores.index(score) + + except ValueError: + msg = 'Unable to get the score index for Tally since {0} ' \ + 'is not one of the bins'.format(bin) + raise ValueError(msg) + + return index + + + def getValue(self, score, filters, filter_bins, nuclide, value='mean'): + """Returns a tally score value given a list of filters to satisfy. + + Parameters + ---------- + score : str + The score string of interest + + filters : list + A list of the filters of interest + + filter_bins : list + A list of the filter bins of interest. These are integers for + material, surface, cell, cellborn, distribcell, universe filters, + and floats for energy or energyout filters. The bins are tuples + of three integers (x,y,z) for mesh filters. The order of the bins + in the list is assumed to correspond to the order of the filters. + + nuclide : Nuclide + The Nuclide of interest + + value : str + A string for the type of value to return ('mean' (default), 'std_dev', + 'sum', or 'sum_sq' are accepted) + """ + + # Determine the score index from the score string + score_index = self._scores.index(score) + + # Determine the nuclide index from the nuclide string/object + nuclide_index = self._nuclides.index(nuclide) + + # Initialize index for Filter in Tally.results[:,:,:] + filter_index = 0 + + # Iterate over specified Filters to compute filter index + for i, filter in enumerate(filters): + + # Find the equivalent Filter in this Tally's list of Filters + test_filter = self.findFilter(filter._type, filter._bins) + + # Filter bins for a mesh are an (x,y,z) tuple + if filter._type == 'mesh': + + # Get the dimensions of the corresponding mesh + nx, ny, nz = test_filter._mesh._dimension + + # Convert (x,y,z) to a single bin -- this is similar to + # subroutine mesh_indices_to_bin in openmc/src/mesh.F90. + value = ((filter_bins[i][0] - 1) * ny * nz + + (filter_bins[i][1] - 1) * nz + + (filter_bins[i][2] - 1)) + filter_index += value * test_filter._stride + + # Filter bins for a distribcell are the "IDs" of each unique placement + # of the Cell in the Geometry (integers starting at 0) + elif filter._type == 'distribcell': + bin = filter_bins[i] + filter_index += bin * test_filter._stride + + else: + bin = filter_bins[i] + bin_index = test_filter.getBinIndex(bin) + filter_index += bin_index * test_filter._stride + + # Return the desired result from Tally + if value == 'mean': + return self._mean[filter_index, nuclide_index, score_index] + elif value == 'std_dev': + return self._std_dev[filter_index, nuclide_index, score_index] + elif value == 'sum': + return self._sum[filter_index, nuclide_index, score_index] + elif value == 'sum_sq': + return self._sum_sq[filter_index, nuclide_index, score_index] + else: + msg = 'Unable to return results from Tally ID={0} for score {1} ' \ + 'since the value {2} is not \'mean\', \'std_dev\', \'sum\', ' \ + 'or \'sum_sq\''.format(self._id, score, value) + raise LookupError(msg) + + + def exportResults(self, filename='tally-results', directory='.', + format='hdf5', append=True): + """Returns a tally score value given a list of filters to satisfy. + + Parameters + ---------- + filename : str + The name of the file for the results (default is 'tally-results') + + directory : str + The name of the directory for the results (default is '.') + + format : str + The format for the exported file - HDF5 ('hdf5', default), Python + pickle ('pkl'), comma-separated values ('csv') files are supported. + + append : bool + Whether or not to append the results to the file (default is True) + """ + + if not is_string(filename): + msg = 'Unable to export the results for Tally ID={0} to filename={1} ' \ + 'since it is not a string'.format(self._id, filename) + raise ValueError(msg) + + elif not is_string(directory): + msg = 'Unable to export the results for Tally ID={0} to directory={1} ' \ + 'since it is not a string'.format(self._id, directory) + raise ValueError(msg) + + elif not format in ['hdf5', 'pkl', 'csv']: + msg = 'Unable to export the results for Tally ID={0} to format {1} ' \ + 'since it is not supported'.format(self._id, format) + raise ValueError(msg) + + elif not isinstance(append, (bool, np.bool)): + msg = 'Unable to export the results for Tally ID={0} since the append ' \ + 'parameters is not True/False'.format(self._id, append) + raise ValueError(msg) + + # Make directory if it does not exist + if not os.path.exists(directory): + os.makedirs(directory) + + + # HDF5 binary file + if format == 'hdf5': + + import h5py + + filename = directory + '/' + filename + '.h5' + + if append: + tally_results = h5py.File(filename, 'a') + else: + tally_results = h5py.File(filename, 'w') + + # Create an HDF5 group within the file for this particular Tally + tally_group = tally_results.create_group('Tally-{0}'.format(self._id)) + + # Add basic Tally data to the HDF5 group + tally_group.create_dataset('id', data=self._id) + tally_group.create_dataset('label', data=self._label) + tally_group.create_dataset('estimator', data=self._estimator) + tally_group.create_dataset('scores', data=np.array(self._scores)) + + # Add a string array of the nuclides to the HDF5 group + nuclides = list() + + for nuclide in self._nuclides: + nuclides.append(nuclide._name) + + + tally_group.create_dataset('nuclides', data=np.array(nuclides)) + + # Create an HDF5 sub-group for the Filters + filter_group = tally_group.create_group('filters') + + for filter in self._filters: + filter_group.create_dataset(filter._type, data=filter._bins) + + # Add all results to the main HDF5 group for the Tally + tally_group.create_dataset('sum', data=self._sum) + tally_group.create_dataset('sum_sq', data=self._sum_sq) + tally_group.create_dataset('mean', data=self._mean) + tally_group.create_dataset('std_dev', data=self._std_dev) + + # Close the Tally results HDF5 file + tally_results.close() + + + # Python pickle binary file + elif format == 'pkl': + + import pickle + + # Load the dictionary from the Pickle file + filename = directory + '/' + filename + '.pkl' + + if os.path.exists(filename) and append: + tally_results = pickle.load(file(filename, 'rb')) + else: + tally_results = dict() + + # Create a nested dictionary within the file for this particular Tally + tally_results['Tally-{0}'.format(self._id)] = dict() + tally_group = tally_results['Tally-{0}'.format(self._id)] + + # Add basic Tally data to the nested dictionary + tally_group['id'] = self._id + tally_group['label'] = self._label + tally_group['estimator'] = self._estimator + tally_group['scores'] = np.array(self._scores) + + # Add a string array of the nuclides to the HDF5 group + nuclides = list() + + for nuclide in self._nuclides: + nuclides.append(nuclide._name) + + tally_group['nuclides']= np.array(nuclides) + + # Create a nested dictionary for the Filters + tally_group['filters'] = dict() + filter_group = tally_group['filters'] + + for filter in self._filters: + filter_group[filter._type] = filter._bins + + # Add all results to the main sub-dictionary for the Tally + tally_group['sum'] = self._sum + tally_group['sum_sq'] = self._sum_sq + tally_group['mean'] = self._mean + tally_group['std_dev'] = self._std_dev + + # Pickle the Tally results to a file + pickle.dump(tally_results, open(filename, 'wb')) + + +class TalliesFile(object): + + def __init__(self): + + # Initialize TalliesFile class attributes + self._tallies = list() + self._meshes = list() + self._tallies_file = ET.Element("tallies") + + + def addTally(self, tally): + + if not isinstance(tally, Tally): + msg = 'Unable to add a non-Tally {0} to the TalliesFile'.format(tally) + raise ValueError(msg) + + self._tallies.append(tally) + + + def removeTally(self, tally): + self._tallies.remove(tally) + + + def addMesh(self, mesh): + + if not isinstance(mesh, Mesh): + msg = 'Unable to add a non-Mesh {0} to the TalliesFile'.format(mesh) + raise ValueError(msg) + + self._meshes.append(mesh) + + + def removeMesh(self, mesh): + self._meshes.remove(mesh) + + + def createTallySubelements(self): + + for tally in self._tallies: + xml_element = tally.getTallyXML() + self._tallies_file.append(xml_element) + + + def createMeshSubelements(self): + + for mesh in self._meshes: + + if len(mesh._name) > 0: + self._tallies_file.append(ET.Comment(mesh._name)) + + xml_element = mesh.getMeshXML() + self._tallies_file.append(xml_element) + + + def exportToXML(self): + + self.createMeshSubelements() + self.createTallySubelements() + + # Clean the indentation in the file to be user-readable + clean_xml_indentation(self._tallies_file) + + # Write the XML Tree to the tallies.xml file + tree = ET.ElementTree(self._tallies_file) + tree.write("tallies.xml", xml_declaration=True, + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py new file mode 100644 index 0000000000..c6c8c5c777 --- /dev/null +++ b/src/utils/openmc/universe.py @@ -0,0 +1,932 @@ +#!/usr/bin/env python + +import openmc +from openmc.checkvalue import * +from xml.etree import ElementTree as ET +from collections import OrderedDict +import operator +import numpy as np + + +################################################################################ +#################################### Cell #################################### +################################################################################ + +# A static variable for auto-generated Cell IDs +AUTO_CELL_ID = 10000 + +def reset_auto_cell_id(): + global AUTO_CELL_ID + AUTO_CELL_ID = 10000 + + +class Cell(object): + + def __init__(self, cell_id=None, name=''): + + # Initialize Cell class attributes + self._id = None + self._name = None + self._fill = None + self._type = None + self._surfaces = dict() + self._rotation = None + self._translation = None + self._offset = None + + self.setId(cell_id) + self.setName(name) + + + def getOffset(self, path, filter_offset): + + # Get the current element and remove it from the list + cell_id = path[0] + path = path[1:] + + # If the Cell is filled by a Material + if self._type == 'normal': + offset = 0 + + # If the Cell is filled by a Universe + elif self._type == 'fill': + offset = self._offset[filter_offset-1] + offset += self._fill.getOffset(path, filter_offset) + + # If the Cell is filled by a Lattice + else: + offset = self._fill.getOffset(path, filter_offset) + + return offset + + # Make a recursive call to the Universe filling this Cell + offset = self._cells[cell_id].getOffset(path, filter_offset) + + # Return the offset computed at all nested Universe levels below this one + return offset + + + def setId(self, cell_id=None): + + if cell_id is None: + global AUTO_CELL_ID + self._id = AUTO_CELL_ID + AUTO_CELL_ID += 1 + + # Check that the ID is an integer and wasn't already used + elif not is_integer(cell_id): + msg = 'Unable to set a non-integer Cell ID {0}'.format(cell_id) + raise ValueError(msg) + + elif cell_id < 0: + msg = 'Unable to set Cell ID to {0} since it must be a ' \ + 'non-negative integer'.format(cell_id) + raise ValueError(msg) + + else: + self._id = cell_id + + + def setName(self, name): + + if not isinstance(name, str): + msg = 'Unable to set name for Cell ID={0} with a non-string ' \ + 'value {1}'.format(self._id, name) + raise ValueError(msg) + + else: + self._name = name + + + def setFill(self, fill): + + if not isinstance(fill, (openmc.Material, Universe, Lattice)) \ + and fill != 'void': + msg = 'Unable to set Cell ID={0} to use a a non-Material or ' \ + 'Universe fill {1}'.format(self._id, fill) + raise ValueError(msg) + + self._fill = fill + + if isinstance(fill, Lattice): + self._type = 'lattice' + elif isinstance(fill, Universe): + self._type = 'fill' + elif fill == 'void': + self._type = 'normal' + else: + self._type = 'normal' + + + def addSurface(self, surface, halfspace): + + if not isinstance(surface, openmc.Surface): + msg = 'Unable to add Surface {0} to Cell ID={1} since it is ' \ + 'not a Surface object'.format(surface, self._id) + raise ValueError(msg) + + if not halfspace in [-1, +1]: + msg = 'Unable to add Surface {0} to Cell ID={1} with halfspace {2} ' \ + 'since it is not +/-1'.format(surface, self._id, halfspace) + raise ValueError(msg) + + # If the Cell does not already contain the Surface, add it + if not surface._id in self._surfaces.keys(): + self._surfaces[surface._id] = (surface, halfspace) + + + def removeSurface(self, surface): + + if not isinstance(surface, openmc.Surface): + msg = 'Unable to remove Surface {0} from Cell ID={1} since it is ' \ + 'not a Surface object'.format(surface, self._id) + raise ValueError(msg) + + # If the Cell contains the Surface, delete it + if surface._id in self._surfaces: + del self._surfaces[surface._id] + + + def setSurfaces(self, surfaces, halfspaces): + + if not isinstance(surfaces, (tuple, list, np.ndarray)): + msg = 'Unable to set Cell ID={0} with Surfaces {1} since it is not a ' \ + 'Python tuple/list or NumPy array'.format(self._id, surfaces) + raise ValueError(msg) + + if not isinstance(halfspaces, (tuple, list, np.ndarray)): + msg = 'Unable to set Cell ID={0} with Surface halfspaces {1} ' \ + 'since it is not a Python tuple/list or NumPy ' \ + 'array'.format(self._id, halfspaces) + raise ValueError(msg) + + for surface in surfaces: + self.addSurface(surface) + + + def setRotation(self, rotation): + + if not isinstance(rotation, (tuple, list, np.ndarray)): + msg = 'Unable to add rotation {0} to Cell ID={1} since it is not a ' \ + 'Python tuple/list or NumPy array'.format(rotation, self._id) + raise ValueError(msg) + + elif len(rotation) != 3: + msg = 'Unable to add rotation {0} to Cell ID={1} since it does not ' \ + 'contain 3 values'.format(rotation, self._id) + raise ValueError(msg) + + for axis in rotation: + + if not is_integer(axis) and not is_float(axis): + msg = 'Unable to add rotation {0} to Cell ID={1} since it is not ' \ + 'an integer or floating point value'.format(axis, self._id) + raise ValueError(msg) + + self._rotation = rotation + + + def setTranslation(self, translation): + + if not isinstance(translation, (tuple, list, np.ndarray)): + msg = 'Unable to add translation {0} to Cell ID={1} since it is not a ' \ + 'Python tuple/list or NumPy array'.format(translation, self._id) + raise ValueError(msg) + + elif len(translation) != 3: + msg = 'Unable to add translation {0} to Cell ID={1} since it does not ' \ + 'contain 3 values'.format(translation, self._id) + raise ValueError(msg) + + for axis in translation: + if not is_integer(axis) and not is_float(axis): + msg = 'Unable to add translation {0} to Cell ID={1} since it is ' \ + 'not an integer or floating point value'.format(axis, self._id) + raise ValueError(msg) + + self._translation = translation + + def setOffset(self, offset): + + if not isinstance(offset, (tuple, list, np.ndarray)): + msg = 'Unable to set offset {0} to Cell ID={1} since it is not a ' \ + 'Python tuple/list or NumPy array'.format(offset, self._id) + raise ValueError(msg) + + self._offset = offset + + + def getAllNuclides(self): + + nuclides = dict() + + if self._type != 'void': + nuclides.update(self._fill.getAllNuclides()) + + return nuclides + + + def getAllCells(self): + + cells = dict() + + if self._type == 'fill' or self._type == 'lattice': + cells.update(self._fill.getAllCells()) + + return cells + + + def getAllUniverses(self): + + universes = dict() + + if self._type == 'fill': + universes[self._fill._id] = self._fill + universes.update(self._fill.getAllUniverses()) + elif self._type == 'lattice': + universes.update(self._fill.getAllUniverses()) + + return universes + + + def __repr__(self): + + string = 'Cell\n' + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + + if isinstance(self._fill, openmc.Material): + string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', self._fill._id) + elif isinstance(self._fill, (Universe, Lattice)): + string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill._id) + else: + string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill) + + string += '{0: <16}{1}\n'.format('\tSurfaces', '=\t') + + for surface_id in self._surfaces: + halfspace = self._surfaces[surface_id][1] + string += '{0} '.format(halfspace * surface_id) + + string = string.rstrip(' ') + '\n' + + string += '{0: <16}{1}{2}\n'.format('\tRotation', '=\t', self._rotation) + string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t', self._translation) + string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offset) + + return string + + + def createXMLSubElement(self, xml_element): + + element = ET.Element("cell") + element.set("id", str(self._id)) + + if isinstance(self._fill, openmc.Material): + element.set("material", str(self._fill._id)) + + elif isinstance(self._fill, (Universe, Lattice)): + element.set("fill", str(self._fill._id)) + self._fill.createXMLSubElement(xml_element) + + else: + element.set("fill", str(self._fill)) + self._fill.createXMLSubElement(xml_element) + + + if not self._surfaces is None: + + surfaces = '' + + for surface_id in self._surfaces.keys(): + + # Determine if XML element already includes subelement for this Surface + path = './surface[@id=\'{0}\']'.format(surface_id) + test = xml_element.find(path) + + # If the element does not contain the Surface subelement, then add it + if test is None: + + # Create the XML subelement for this Surface and add it to the element + surface = self._surfaces[surface_id][0] + surface_subelement = surface.createXMLSubElement() + + if len(surface._name) > 0: + xml_element.append(ET.Comment(surface._name)) + + xml_element.append(surface_subelement) + + # Append to the Surfaces's XML attribute the halfspace and Surface ID + halfspace = self._surfaces[surface_id][1] + surfaces += '{0} '.format(halfspace * surface_id) + + element.set("surfaces", surfaces.rstrip(' ')) + + + if not self._translation is None: + + translation = '' + + for axis in self._translation: + translation += '{0} '.format(axis) + + element.set("translation", translation.rstrip(' ')) + + + if not self._rotation is None: + + rotation = '' + + for axis in self._rotation: + rotation += '{0} '.format(axis) + + element.set("rotation", rotation.rstrip(' ')) + + return element + + + +################################################################################ +################################### Universe ################################# +################################################################################ + +# A static variable for auto-generated Lattice (Universe) IDs +AUTO_UNIVERSE_ID = 10000 + +def reset_auto_universe_id(): + global AUTO_UNIVERSE_ID + AUTO_UNIVERSE_ID = 10000 + + +class Universe(object): + + def __init__(self, universe_id=None, name=''): + + # Initialize Cell class attributes + self._id = None + self._name = None + + # Keys - Cell IDs + # Values - Cells + self._cells = dict() + + # Keys - Cell IDs + # Values - Offsets + self._cell_offsets = OrderedDict() + self._num_regions = 0 + + self.setId(universe_id) + self.setName(name) + + + def setId(self, universe_id=None): + + if universe_id is None: + global AUTO_UNIVERSE_ID + self._id = AUTO_UNIVERSE_ID + AUTO_UNIVERSE_ID += 1 + + # Check that the ID is an integer and wasn't already used + elif not is_integer(universe_id): + msg = 'Unable to set Universe ID to a non-integer {0}'.format(universe_id) + raise ValueError(msg) + + elif universe_id < 0: + msg = 'Unable to set Universe ID to {0} since it must be a ' \ + 'non-negative integer'.format(universe_id) + raise ValueError(msg) + + else: + self._id = universe_id + + + def setName(self, name): + + if not is_string(name): + msg = 'Unable to set name for Universe ID={0} with a non-string ' \ + 'value {1}'.format(self._id, name) + raise ValueError(msg) + + else: + self._name = name + + + def addCell(self, cell): + + if not isinstance(cell, Cell): + msg = 'Unable to add a Cell to Universe ID={0} since {1} is not ' \ + 'a Cell'.format(self._id, cell) + raise ValueError(msg) + + cell_id = cell._id + + if not cell_id in self._cells.keys(): + self._cells[cell_id] = cell + + + def addCells(self, cells): + + if not isinstance(cells, (list, tuple, np.ndarray)): + msg = 'Unable to add Cells to Universe ID={0} since {1} is not a ' \ + 'Python tuple/list or NumPy array'.format(self._id, cells) + raise ValueError(msg) + + for i in range(len(cells)): + self.addCell(cells[i]) + + + def removeCell(self, cell): + + if not isinstance(cell, Cell): + msg = 'Unable to remove a Cell from Universe ID={0} since {1} is ' \ + 'not a Cell'.format(self._id, cell) + raise ValueError(msg) + + cell_id = cell.getId() + + # If the Cell is in the Universe's list of Cells, delete it + if cell_id in self._cells.keys(): + del self._cells[cell_id] + + + def clearCells(self): + self._cells.clear() + + + def getOffset(self, path, filter_offset): + + # Get the current element and remove it from the list + path = path[1:] + + # Get the Cell ID + cell_id = path[0] + + # Make a recursive call to the Cell within this Universe + offset = self._cells[cell_id].getOffset(path, filter_offset) + + # Return the offset computed at all nested Universe levels below this one + return offset + + + def getAllNuclides(self): + + nuclides = dict() + + # Append all Nuclides in each Cell in the Universe to the dictionary + for cell_id, cell in self._cells.items(): + nuclides.update(cell.getAllNuclides()) + + return nuclides + + + def getAllCells(self): + + cells = dict() + + # Add this Universe's cells to the dictionary + cells.update(self._cells) + + # Append all Cells in each Cell in the Universe to the dictionary + for cell_id, cell in self._cells.items(): + cells.update(cell.getAllCells()) + + return cells + + + def getAllUniverses(self): + + # Get all Cells in this Universe + cells = self.getAllCells() + + universes = dict() + + # Append all Universes containing each Cell to the dictionary + for cell_id, cell in cells.items(): + universes.update(cell.getAllUniverses()) + + return universes + + + def __repr__(self): + + string = 'Universe\n' + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <16}{1}{2}\n'.format('\tCells', '=\t', self._cells.keys()) + string += '{0: <16}{1}{2}\n'.format('\t# Regions', '=\t', self._num_regions) + return string + + + def createXMLSubElement(self, xml_element): + + # Iterate over all Cells + for cell_id, cell in self._cells.items(): + + # Determine if XML element already contains subelement for this Cell + path = './cell[@id=\'{0}\']'.format(cell_id) + test = xml_element.find(path) + + # If the element does not contain the Cell subelement, then add it + if test is None: + + # Create the XML subelement for this Cell and everything beneath it + cell_subelement = cell.createXMLSubElement(xml_element) + + # Append the Universe ID to the subelement and add it to the element + cell_subelement.set("universe", str(self._id)) + + if len(cell._name) > 0: + xml_element.append(ET.Comment(cell._name)) + + xml_element.append(cell_subelement) + + + +################################################################################ +################################### Lattice ################################## +################################################################################ + + +class Lattice(object): + + def __init__(self, lattice_id=None, name='', type='rectangular'): + + # Initialize Lattice class attributes + self._id = None + self._name = None + self._type = '' + self._dimension = None + self._lower_left = None + self._width = None + self._outside = None + self._universes = None + self._offsets = None + + + self.setId(lattice_id) + self.setName(name) + self.setType(type) + + + def setId(self, lattice_id=None): + + if lattice_id is None: + global AUTO_UNIVERSE_ID + self._id = AUTO_UNIVERSE_ID + AUTO_UNIVERSE_ID += 1 + + # Check that the ID is an integer and wasn't already used + elif not is_integer(lattice_id): + msg = 'Unable to set a non-integer Lattice ID {0}'.format(lattice_id) + raise ValueError(msg) + + elif lattice_id < 0: + msg = 'Unable to set Lattice ID to {0} since it must be a ' \ + 'non-negative integer'.format(lattice_id) + raise ValueError(msg) + + else: + self._id = lattice_id + + + def setName(self, name): + + if not is_string(name): + msg = 'Unable to set name for Lattice ID={0} with a non-string ' \ + 'value {1}'.format(self._id, name) + raise ValueError(msg) + + else: + self._name = name + + + def setType(self, type): + + if not is_string(type): + msg = 'Unable to set type for Lattice ID={0} with a non-string ' \ + 'value {1}'.format(self._id, type) + raise ValueError(msg) + + elif not type in ['rectangular', 'hexagonal']: + msg = 'Unable to set type for Lattice ID={0} as {1} since it is not ' \ + 'rectangular or hexagonal'.format(self._id, type) + raise ValueError(msg) + + self._type = type + + + def setDimension(self, dimension): + + if not isinstance(dimension, (tuple, list, np.ndarray)): + msg = 'Unable to set Lattice ID={0} dimension to {1} since it is not ' \ + 'a Python tuple/list or NumPy array'.format(self._id, dimension) + raise ValueError(msg) + + elif len(dimension) != 2 and len(dimension) != 3: + msg = 'Unable to set Lattice ID={0} dimension to {1} since it does ' \ + 'not contain 2 or 3 coordinates'.format(self._id, dimension) + raise ValueError(msg) + + for dim in dimension: + + if not is_integer(dim) and not is_float(dim): + msg = 'Unable to set Lattice ID={0} dimension to {1} since it is ' \ + 'not an integer or floating point value'.format(self._id, dim) + raise ValueError(msg) + + + elif dim < 0: + msg = 'Unable to set Lattice ID={0} dimension to {1} since it ' \ + 'is a negative value'.format(self._id, dim) + raise ValueError(msg) + + self._dimension = dimension + + + def setLowerLeft(self, lower_left): + + if not isinstance(lower_left, (tuple, list, np.ndarray)): + msg = 'Unable to set Lattice ID={0} lower_left to {1} since it is not ' \ + 'a Python tuple/list 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 Lattice ID={0} lower_left to {1} since it does ' \ + 'not contain 2 or 3 coordinates'.format(self._id, lower_left) + raise ValueError(msg) + + + for dim in lower_left: + + if not is_integer(dim) and not is_float(dim): + msg = 'Unable to set Lattice ID={0} lower_left to {1} since it is ' \ + 'is not an integer or floating point value'.format(self._id, dim) + raise ValueError(msg) + + self._lower_left = lower_left + + + def setWidth(self, width): + + if not isinstance(width, (tuple, list, np.ndarray)): + msg = 'Unable to set Lattice ID={0} width to {1} since it 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 set Lattice ID={0} width to {1} since it does ' \ + 'not contain 2 or 3 coordinates'.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 Lattice ID={0} width to {1} since it is not an ' \ + 'an integer or floating point value'.format(self._id, dim) + raise ValueError(msg) + + elif dim < 0: + msg = 'Unable to set Lattice ID={0} width to {1} since it ' \ + 'is a negative value'.format(self._id, dim) + raise ValueError(msg) + + self._width = width + + + def setOutside(self, outside): + + if not isinstance(outside, Universe): + msg = 'Unable to set Lattice ID={0} outside universe to {1} ' \ + 'since it is not a Universe object'.format(self._id, outside) + raise ValueError(msg) + + self._outside = outside + + + def setUniverses(self, universes): + + if not isinstance(universes, (tuple, list, np.ndarray)): + msg = 'Unable to set Lattice ID={0} universes to {1} since it is not ' \ + 'a Python tuple/list or NumPy array'.format(self._id, universes) + raise ValueError(msg) + + self._universes = np.asarray(universes, dtype=Universe) + + + def setOffsets(self, offsets): + + if not isinstance(offsets, (tuple, list, np.ndarray)): + msg = 'Unable to set Lattice ID={0} offsets to {1} since it is not ' \ + 'a Python tuple/list or NumPy array'.format(self._id, offsets) + raise ValueError(msg) + + self._offsets = offsets + + + def getOffset(self, path, filter_offset): + + # Get the current element and remove it from the list + i = path[0] + path = path[1:] + + # For 2D Lattices + if len(self._dimension) == 2: + offset = self._offsets[i[1]-1, i[2]-1, 0, filter_offset-1] + offset += self._universes[i[1]][i[2]].getOffset(path, filter_offset) + + # For 3D Lattices + else: + offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, filter_offset-1] + offset += self._universes[i[1]-1][i[2]-1][i[3]-1].getOffset(path, \ + filter_offset) + + return offset + + + def getUniqueUniverses(self): + + unique_universes = np.unique(self._universes.ravel()) + universes = dict() + + for universe in unique_universes: + universes[universe._id] = universe + + return universes + + + def getAllNuclides(self): + + nuclides = dict() + + # Get all unique Universes contained in each of the lattice cells + unique_universes = self.getUniqueUniverses() + + # Append all Universes containing each cell to the dictionary + for universe_id, universe in unique_universes.items(): + nuclides.update(universe.getAllNuclides()) + + return nuclides + + + def getAllCells(self): + + cells = dict() + unique_universes = self.getUniqueUniverses() + + for universe_id, universe in unique_universes.items(): + cells.update(universe.getAllCells()) + + return cells + + + def getAllUniverses(self): + + # Initialize a dictionary of all Universes contained by the Lattice + # in each nested Universe level + all_universes = dict() + + # Get all unique Universes contained in each of the lattice cells + unique_universes = self.getUniqueUniverses() + + # Add the unique Universes filling each Lattice cell + all_universes.update(unique_universes) + + # Append all Universes containing each cell to the dictionary + for universe_id, universe in unique_universes.items(): + all_universes.update(universe.getAllUniverses()) + + return all_universes + + + def __repr__(self): + + string = 'Lattice\n' + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) + string += '{0: <16}{1}{2}\n'.format('\tDimension', '=\t', self._dimension) + string += '{0: <16}{1}{2}\n'.format('\tLower Left', '=\t', self._lower_left) + string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._width) + + if self._outside is not None: + string += '{0: <16}{1}{2}\n'.format('\tOutside', '=\t', self._outside._id) + else: + string += '{0: <16}{1}{2}\n'.format('\tOutside', '=\t', self._outside) + + string += '{0: <16}\n'.format('\tUniverses') + + # Lattice nested Universe IDs - column major for Fortran + for i, universe in enumerate(np.ravel(self._universes)): + string += '{0} '.format(universe._id) + + # Add a newline character every time we reach the end of a row of cells + if (i+1) % self._dimension[-1] == 0: + string += '\n' + + string = string.rstrip('\n') + + if self._offsets is not None: + string += '{0: <16}\n'.format('\tOffsets') + + # Lattice cell offsets + for i, offset in enumerate(np.ravel(self._offsets)): + + string += '{0} '.format(offset) + + # Add a newline character every time we reach the end of a row of cells + if (i+1) % self._dimension[-1] == 0: + string += '\n' + + string = string.rstrip('\n') + + return string + + + def createXMLSubElement(self, xml_element): + + # Determine if XML element already contains subelement for this Lattice + path = './lattice[@id=\'{0}\']'.format(self._id) + test = xml_element.find(path) + + # If the element does contain the Lattice subelement, then return + if not test is None: + return + + lattice_subelement = ET.Element("lattice") + lattice_subelement.set("id", str(self._id)) + lattice_subelement.set("type", self._type) + + # Export Lattice cell dimensions + dimension = ET.SubElement(lattice_subelement, "dimension") + dimension.text = '{0} {1} {2}'.format(self._dimension[0], \ + self._dimension[1], \ + self._dimension[2]) + + # Export Lattice lower left + lower_left = ET.SubElement(lattice_subelement, "lower_left") + lower_left.text = '{0} {1} {2}'.format(self._lower_left[0], \ + self._lower_left[1], \ + self._lower_left[2]) + + # Export the Lattice cell width/height + width = ET.SubElement(lattice_subelement, "width") + width.text = '{0} {1} {2}'.format(self._width[0], \ + self._width[1], \ + self._width[2]) + + # Export the Lattice outside Universe (if specified) + if self._outside is not None: + outside = ET.SubElement(lattice_subelement, "outside") + outside.text = '{0}'.format(self._outside._id) + + # Export the Lattice nested Universe IDs - column major for Fortran + universe_ids = '\n' + + # 3D Lattices + if len(self._dimension) == 3: + for z in range(self._dimension[2]): + for y in range(self._dimension[1]): + for x in range(self._dimension[0]): + + universe = self._universes[x][y][z] + + # Append Universe ID to the string for the Lattice XML subelement + universe_ids += '{0} '.format(universe._id) + + # Create XML subelement for this Universe and everything beneath it + universe.createXMLSubElement(xml_element) + + # Add newline character every time we reach the end of row of cells + universe_ids += '\n' + + # Add newline character every time we reach the end of row of cells + universe_ids += '\n' + + # 2D Lattices + else: + for y in range(self._dimension[1]): + for x in range(self._dimension[0]): + + universe = self._universes[x][y] + + # Append Universe ID to the string for the Lattice XML subelement + universe_ids += '{0} '.format(universe._id) + + # Create XML subelement for this Universe and everything beneath it + universe.createXMLSubElement(xml_element) + + # Add newline character every time we reach the end of row of cells + universe_ids += '\n' + + # Remove trailing newline character from Universe IDs string + universe_ids = universe_ids.rstrip('\n') + + universes = ET.SubElement(lattice_subelement, "universes") + universes.text = universe_ids + + if len(self._name) > 0: + xml_element.append(ET.Comment(self._name)) + + # Append the XML subelement for this Lattice to the XML element + xml_element.append(lattice_subelement) diff --git a/src/utils/setup.py b/src/utils/setup.py index 1000e509b3..f706ce9688 100644 --- a/src/utils/setup.py +++ b/src/utils/setup.py @@ -9,3 +9,11 @@ setup(name='statepoint', author_email='paul.k.romano@gmail.com', url='https://github.com/mit-crpg/openmc', py_modules=['statepoint']) + +setup(name='openmc', + version='0.6.0', + description='OpenMC Python API', + author='Will Boyd', + author_email='wbinventor@gmail.com', + url='https://github.com/mit-crpg/openmc', + packages=['openmc'])