This commit is contained in:
Adam Nelson 2016-08-11 19:18:16 -04:00
commit daa8fbc3f6
46 changed files with 5105 additions and 247 deletions

View file

@ -30,3 +30,4 @@ Output Files
particle_restart
track
voxel
volume

View file

@ -1,4 +1,4 @@
.. _usersguide_nuclear_data:
.. _io_nuclear_data:
========================
Nuclear Data File Format

View file

@ -0,0 +1,22 @@
.. _io_volume:
==================
Volume File Format
==================
**/**
:Attributes: - **samples** (*int*) -- Number of samples
- **lower_left** (*double[3]*) -- Lower-left coordinates of
bounding box
- **upper_right** (*double[3]*) -- Upper-right coordinates of
bounding box
**/cell_<id>/**
:Datasets: - **volume** (*double[2]*) -- Calculated volume and its uncertainty
in cubic centimeters
- **nuclides** (*char[][]*) -- Names of nuclides identified in the
cell
- **atoms** (*double[][2]*) -- Total number of atoms of each nuclide
and its uncertainty

View file

@ -324,6 +324,16 @@ Functions
:mod:`openmc.data` -- Nuclear Data Interface
--------------------------------------------
Physical Data
-------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.data.atomic_mass
Core Classes
------------

View file

@ -836,6 +836,35 @@ displayed. This element takes the following attributes:
*Default*: 5
``<volume_calc>`` Element
-------------------------
The ``<volume_calc>`` element indicates that a stochastic volume calculation
should be run at the beginning of the simulation. This element has the following
sub-elements/attributes:
:cells:
The unique IDs of cells for which the volume should be estimated.
*Default*: None
:samples:
The number of samples used to estimate volumes.
*Default*: None
:lower_left:
The lower-left Cartesian coordinates of a bounding box that is used to
sample points within.
*Default*: None
:upper_right:
The upper-right Cartesian coordinates of a bounding box that is used to
sample points within.
*Default*: None
--------------------------------------
Geometry Specification -- geometry.xml
--------------------------------------

View file

@ -3,7 +3,7 @@
openmc \- Executes the OpenMC Monte Carlo code
.SH DESCRIPTION
This command is used to execute the OpenMC Monte Carlo code. It is assumed that
a set of XML input files has already been created and that ACE format cross
a set of XML input files has already been created and that HDF5 format cross
sections are available.
.SH SYNOPSIS
\fBopenmc\fR [\fIoptions\fR] [\fIpath\fR]
@ -40,11 +40,21 @@ The behavior of
.B openmc
is affected by the following environment variables.
.TP
.B CROSS_SECTIONS
.B OPENMC_CROSS_SECTIONS
Indicates the default path to the cross_sections.xml summary file that is used
to locate ACE format cross section libraries if the user has not specified the
to locate HDF5 format cross section libraries if the user has not specified the
<cross_sections> tag in
.I settings.xml\fP.
.TP
.B OPENMC_MG_CROSS_SECTIONS
Indicates the default path to the mgxs.xml file that contains multi-group cross
section libraries if the user has not specified the <cross_sections> tag in
.I settings.xml\fP.
.TP
.B OPENMC_MULTIPOLE_LIBRARY
Indicates the default path to a directory containing windowed multipole data if
the user has not specified the <multipole_library> tag in
.I settings.xml\fP.
.SH LICENSE
Copyright \(co 2011-2016 Massachusetts Institute of Technology.
.PP

View file

@ -6,6 +6,9 @@ from openmc.nuclide import *
from openmc.macroscopic import *
from openmc.material import *
from openmc.plots import *
from openmc.region import *
from openmc.volume import *
from openmc.source import *
from openmc.settings import *
from openmc.surface import *
from openmc.universe import *
@ -18,8 +21,6 @@ from openmc.cmfd import *
from openmc.executor import *
from openmc.statepoint import *
from openmc.summary import *
from openmc.region import *
from openmc.source import *
from openmc.particle_restart import *
try:

View file

@ -85,6 +85,10 @@ class Cell(object):
Array of offsets used for distributed cell searches
distribcell_index : int
Index of this cell in distribcell arrays
volume_information : dict
Estimate of the volume and total number of atoms of each nuclide from a
stochastic volume calculation. This information is set with the
:meth:`Cell.add_volume_information` method.
"""
@ -100,6 +104,7 @@ class Cell(object):
self._translation = None
self._offsets = None
self._distribcell_index = None
self._volume_information = None
def __contains__(self, point):
if self.region is None:
@ -212,6 +217,10 @@ class Cell(object):
def distribcell_index(self):
return self._distribcell_index
@property
def volume_information(self):
return self._volume_information
@id.setter
def id(self, cell_id):
if cell_id is None:
@ -351,6 +360,25 @@ class Cell(object):
else:
self.region = Intersection(self.region, region)
def add_volume_information(self, volume_calc):
"""Add volume information to a cell.
Parameters
----------
volume_calc : openmc.VolumeCalculation
Results from a stochastic volume calculation
"""
if volume_calc.domain_type == 'cell':
for cell_id in volume_calc.results:
if cell_id == self.id:
self._volume_information = volume_calc.results[cell_id]
break
else:
raise ValueError('No volume information found for this cell.')
else:
raise ValueError('No volume information found for this cell.')
def get_cell_instance(self, path, distribcell_index):
# If the Cell is filled by a Material
@ -368,8 +396,19 @@ class Cell(object):
return offset
def get_all_nuclides(self):
"""Return all nuclides contained in the cell
def get_nuclides(self):
"""Returns all nuclides in the cell
Returns
-------
nuclides : list of str
List of nuclide names
"""
return self.fill.get_nuclides() if self.fill_type != 'void' else []
def get_nuclide_densities(self):
"""Return all nuclides contained in the cell and their densities
Returns
-------
@ -381,8 +420,24 @@ class Cell(object):
nuclides = OrderedDict()
if self.fill_type != 'void':
nuclides.update(self.fill.get_all_nuclides())
if self.fill_type == 'material':
nuclides.update(self.fill.get_nuclide_densities())
elif self.fill_type == 'void':
pass
else:
if self.volume_information is not None:
volume = self.volume_information['volume'][0]
for full_name, atoms in self.volume_information['atoms']:
name, xs = full_name.split('.')
nuclide = openmc.Nuclide(name, xs)
density = 1.0e-24 * atoms[0]/volume # density in atoms/b-cm
nuclides[name] = (nuclide, density)
else:
raise RuntimeError(
'Volume information is needed to calculate microscopic cross '
'sections for cell {}. This can be done by running a '
'stochastic volume calculation via the '
'openmc.VolumeCalculation object'.format(self.id))
return nuclides

View file

@ -282,7 +282,7 @@ class Library(object):
lines = [ace_file.readline() for i in range(13)]
while len(lines) != 0 and lines[0] != '':
while len(lines) != 0 and lines[0].strip() != '':
# Read name of table, atomic mass ratio, and temperature. If first
# line is empty, we are at end of file

View file

@ -1,3 +1,7 @@
import itertools
import os
# Isotopic abundances from M. Berglund and M. E. Wieser, "Isotopic compositions
# of the elements 2009 (IUPAC Technical Report)", Pure. Appl. Chem. 83 (2),
# pp. 397--410 (2011).
@ -119,9 +123,12 @@ ATOMIC_SYMBOL = {1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C', 7: 'N',
98: 'Cf', 99: 'Es', 100: 'Fm', 101: 'Md', 102: 'No',
103: 'Lr', 104: 'Rf', 105: 'Db', 106: 'Sg', 107: 'Bh',
108: 'Hs', 109: 'Mt', 110: 'Ds', 111: 'Rg', 112: 'Cn',
114: 'Fl', 116: 'Lv'}
113: 'Nh', 114: 'Fl', 115: 'Mc', 116: 'Lv', 117: 'Ts',
118: 'Og'}
ATOMIC_NUMBER = {value: key for key, value in ATOMIC_SYMBOL.items()}
_ATOMIC_MASS = {}
REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
5: '(n,misc)', 11: '(n,2nd)', 16: '(n,2n)', 17: '(n,3n)',
18: '(n,fission)', 19: '(n,f)', 20: '(n,nf)', 21: '(n,2nf)',
@ -175,3 +182,39 @@ SUM_RULES = {1: [2, 3],
105: list(range(700, 750)),
106: list(range(750, 800)),
107: list(range(800, 850))}
def atomic_mass(isotope):
"""Return atomic mass of isotope in atomic mass units.
Atomic mass data comes from the Atomic Mass Evaluation 2012, published in
Chinese Physics C 36 (2012), 1287--1602.
Parameters
----------
isotope : str
Name of isotope, e.g. 'Pu239'
Returns
-------
float or None
Atomic mass of isotope in atomic mass units. If the isotope listed does
not have a known atomic mass, None is returned.
"""
if not _ATOMIC_MASS:
# Load data from AME2012 file
mass_file = os.path.join(os.path.dirname(__file__), 'mass.mas12')
with open(mass_file, 'r') as ame:
# Read lines in file starting at line 40
for line in itertools.islice(ame, 40, None):
name = '{}{}'.format(line[20:22].strip(), int(line[16:19]))
mass = float(line[96:99]) + 1e-6*float(
line[100:106] + '.' + line[107:112])
_ATOMIC_MASS[name.lower()] = mass
# Get rid of metastable information
if '_' in isotope:
isotope = isotope[:isotope.find('_')]
return _ATOMIC_MASS.get(isotope.lower())

View file

@ -97,7 +97,7 @@ class Tabulated1D(object):
if self.interpolation[k] == 1:
# Histogram
y[contined] = yi
y[contained] = yi
elif self.interpolation[k] == 2:
# Linear-linear

3392
openmc/data/mass.mas12 Normal file

File diff suppressed because it is too large Load diff

View file

@ -178,7 +178,7 @@ class IncidentNeutron(object):
@temperature.setter
def temperature(self, temperature):
cv.check_type('temperature', temperature, Real)
cv.check_greater_than('temperature', temperature, 0.0)
cv.check_greater_than('temperature', temperature, 0.0, True)
self._temperature = temperature
@energy.setter
@ -448,6 +448,10 @@ class IncidentNeutron(object):
# Create summed reaction with appropriate cross section
rx = Reaction(mt)
mts = data.get_reaction_components(mt)
if len(mts) == 0:
warn('Photon production is present for MT={} but no '
'reaction components exist.'.format(mt))
continue
rx.xs = Sum([data.reactions[mt_i].xs for mt_i in mts])
# Determine summed cross section

View file

@ -16,7 +16,8 @@ from openmc.stats import Discrete, Tabular
_THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27',
'be': 'c_Be',
'bebeo': 'c_Be_in_BeO', 'be-o': 'c_Be_in_BeO',
'beo': 'c_BeO',
'bebeo': 'c_Be_in_BeO', 'be-o': 'c_Be_in_BeO', 'be/o': 'c_Be_in_BeO',
'benz': 'c_Benzine',
'cah': 'c_Ca_in_CaH2',
'dd2o': 'c_D_in_D2O', 'hwtr': 'c_D_in_D2O',
@ -25,18 +26,18 @@ _THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27',
'hca': 'c_H_in_CaH2',
'hch2': 'c_H_in_CH2', 'poly': 'c_H_in_CH2',
'hh2o': 'c_H_in_H2O', 'lwtr': 'c_H_in_H2O',
'hzrh': 'c_H_in_ZrH', 'h-zr': 'c_H_in_ZrH',
'hzrh': 'c_H_in_ZrH', 'h-zr': 'c_H_in_ZrH', 'h/zr': 'c_H_in_ZrH',
'lch4': 'c_liquid_CH4', 'lmeth': 'c_liquid_CH4',
'mg': 'c_Mg24',
'obeo': 'c_O_in_BeO', 'o-be': 'c_O_in_BeO',
'obeo': 'c_O_in_BeO', 'o-be': 'c_O_in_BeO', 'o/be': 'c_O_in_BeO',
'orthod': 'c_ortho_D', 'dortho': 'c_ortho_D',
'orthoh': 'c_ortho_H', 'hortho': 'c_ortho_H',
'ouo2': 'c_O_in_UO2', 'o2-u': 'c_O_in_UO2',
'ouo2': 'c_O_in_UO2', 'o2-u': 'c_O_in_UO2', 'o2/u': 'c_O_in_UO2',
'parad': 'c_para_D', 'dpara': 'c_para_D',
'parah': 'c_para_H', 'hpara': 'c_para_H',
'sch4': 'c_solid_CH4', 'smeth': 'c_solid_CH4',
'uuo2': 'c_U_in_UO2', 'u-o2': 'c_U_in_UO2',
'zrzrh': 'c_Zr_in_ZrH', 'zr-h': 'c_Zr_in_ZrH'}
'uuo2': 'c_U_in_UO2', 'u-o2': 'c_U_in_UO2', 'u/o2': 'c_U_in_UO2',
'zrzrh': 'c_Zr_in_ZrH', 'zr-h': 'c_Zr_in_ZrH', 'zr/h': 'c_Zr_in_ZrH'}
class CoherentElastic(object):

View file

@ -50,6 +50,20 @@ class Geometry(object):
self._root_universe = root_universe
def add_volume_information(self, volume_calc):
"""Add volume information from a stochastic volume calculation.
Parameters
----------
volume_calc : openmc.VolumeCalculation
Results from a stochastic volume calculation
"""
if volume_calc.domain_type == 'cell':
for cell in self.get_all_cells():
if cell.id in volume_calc.results:
cell.add_volume_information(volume_calc)
def export_to_xml(self):
"""Create a geometry.xml file that can be used for a simulation.
@ -166,24 +180,6 @@ class Geometry(object):
universes.sort(key=lambda x: x.id)
return universes
def get_all_nuclides(self):
"""Return all nuclides assigned to a material in the geometry
Returns
-------
list of openmc.Nuclide
Nuclides in the geometry
"""
nuclides = OrderedDict()
materials = self.get_all_materials()
for material in materials:
nuclides.update(material.get_all_nuclides())
return nuclides
def get_all_materials(self):
"""Return all materials assigned to a cell

View file

@ -144,25 +144,26 @@ class Lattice(object):
return univs
def get_all_nuclides(self):
"""Return all nuclides contained in the lattice
def get_nuclides(self):
"""Returns all nuclides in the lattice
Returns
-------
nuclides : collections.OrderedDict
Dictionary whose keys are nuclide names and values are 2-tuples of
(nuclide, density)
nuclides : list of str
List of nuclide names
"""
nuclides = OrderedDict()
nuclides = []
# Get all unique Universes contained in each of the lattice cells
unique_universes = self.get_unique_universes()
# Append all Universes containing each cell to the dictionary
for universe_id, universe in unique_universes.items():
nuclides.update(universe.get_all_nuclides())
for universe in unique_universes.values():
for nuclide in universe.get_nuclides():
if nuclide not in nuclides:
nuclides.append(nuclide)
return nuclides

View file

@ -492,9 +492,31 @@ class Material(object):
for element, percent, percent_type in self._elements:
element.scattering = 'iso-in-lab'
def get_all_nuclides(self):
def get_nuclides(self):
"""Returns all nuclides in the material
Returns
-------
nuclides : list of str
List of nuclide names
"""
nuclides = []
for nuclide, density, density_type in self._nuclides:
nuclides.append(nuclide.name)
for element, density, density_type in self._elements:
# Expand natural element into isotopes
for isotope, abundance in element.expand():
nuclides.append(isotope.name)
return nuclides
def get_nuclide_densities(self):
"""Returns all nuclides in the material and their densities
Returns
-------
nuclides : dict

View file

@ -1069,7 +1069,7 @@ class Library(object):
# Create the xsdata object and add it to the mgxs_file
for i, domain in enumerate(self.domains):
if self.by_nuclide:
nuclides = list(domain.get_all_nuclides().keys())
nuclides = domain.get_nuclides()
else:
nuclides = ['total']
for nuclide in nuclides:

View file

@ -288,9 +288,7 @@ class MGXS(object):
# If this is a by-nuclide cross-section, add nuclides to Tally
if self.by_nuclide and score != 'flux':
all_nuclides = self.get_all_nuclides()
for nuclide in all_nuclides:
self._tallies[key].nuclides.append(nuclide)
self._tallies[key].nuclides += self.get_nuclides()
else:
self._tallies[key].nuclides.append('total')
@ -329,14 +327,14 @@ class MGXS(object):
@property
def num_nuclides(self):
if self.by_nuclide:
return len(self.get_all_nuclides())
return len(self.get_nuclides())
else:
return 1
@property
def nuclides(self):
if self.by_nuclide:
return self.get_all_nuclides()
return self.get_nuclides()
else:
return 'sum'
@ -503,7 +501,7 @@ class MGXS(object):
mgxs.name = name
return mgxs
def get_all_nuclides(self):
def get_nuclides(self):
"""Get all nuclides in the cross section's spatial domain.
Returns
@ -528,8 +526,7 @@ class MGXS(object):
# Otherwise, return all nuclides in the spatial domain
else:
nuclides = self.domain.get_all_nuclides()
return list(nuclides.keys())
return self.domain.get_nuclides()
def get_nuclide_density(self, nuclide):
"""Get the atomic number density in units of atoms/b-cm for a nuclide
@ -556,7 +553,7 @@ class MGXS(object):
cv.check_type('nuclide', nuclide, basestring)
# Get list of all nuclides in the spatial domain
nuclides = self.domain.get_all_nuclides()
nuclides = self.domain.get_nuclide_densities()
if nuclide not in nuclides:
msg = 'Unable to get density for nuclide "{0}" which is not in ' \
@ -597,14 +594,14 @@ class MGXS(object):
# Sum the atomic number densities for all nuclides
if nuclides == 'sum':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
densities = np.zeros(1, dtype=np.float)
for nuclide in nuclides:
densities[0] += self.get_nuclide_density(nuclide)
# Tabulate the atomic number densities for all nuclides
elif nuclides == 'all':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
densities = np.zeros(self.num_nuclides, dtype=np.float)
for i, nuclide in enumerate(nuclides):
densities[i] += self.get_nuclide_density(nuclide)
@ -635,7 +632,7 @@ class MGXS(object):
# If computing xs for each nuclide, replace CrossNuclides with originals
if self.by_nuclide:
self.xs_tally._nuclides = []
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
for nuclide in nuclides:
self.xs_tally.nuclides.append(openmc.Nuclide(nuclide))
@ -797,7 +794,7 @@ class MGXS(object):
# Construct a collection of the nuclides to retrieve from the xs tally
if self.by_nuclide:
if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']:
query_nuclides = self.get_all_nuclides()
query_nuclides = self.get_nuclides()
else:
query_nuclides = nuclides
else:
@ -1165,7 +1162,7 @@ class MGXS(object):
# Construct a collection of the nuclides to report
if self.by_nuclide:
if nuclides == 'all':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
elif nuclides == 'sum':
nuclides = ['sum']
else:
@ -1303,7 +1300,7 @@ class MGXS(object):
# Construct a collection of the nuclides to report
if self.by_nuclide:
if nuclides == 'all':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
densities = np.zeros(len(nuclides), dtype=np.float)
elif nuclides == 'sum':
nuclides = ['sum']
@ -1485,7 +1482,7 @@ class MGXS(object):
if self.by_nuclide and nuclides == 'sum':
# Use tally summation to sum across all nuclides
query_nuclides = self.get_all_nuclides()
query_nuclides = self.get_nuclides()
xs_tally = self.xs_tally.summation(nuclides=query_nuclides)
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
@ -1791,7 +1788,7 @@ class MatrixMGXS(MGXS):
# Construct a collection of the nuclides to retrieve from the xs tally
if self.by_nuclide:
if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']:
query_nuclides = self.get_all_nuclides()
query_nuclides = self.get_nuclides()
else:
query_nuclides = nuclides
else:
@ -1938,7 +1935,7 @@ class MatrixMGXS(MGXS):
# Construct a collection of the nuclides to report
if self.by_nuclide:
if nuclides == 'all':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
if nuclides == 'sum':
nuclides = ['sum']
else:
@ -3623,7 +3620,7 @@ class ScatterMatrixXS(MatrixMGXS):
# Construct a collection of the nuclides to retrieve from the xs tally
if self.by_nuclide:
if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']:
query_nuclides = self.get_all_nuclides()
query_nuclides = self.get_nuclides()
else:
query_nuclides = nuclides
else:
@ -3790,7 +3787,7 @@ class ScatterMatrixXS(MatrixMGXS):
# Construct a collection of the nuclides to report
if self.by_nuclide:
if nuclides == 'all':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
if nuclides == 'sum':
nuclides = ['sum']
else:
@ -4596,7 +4593,7 @@ class Chi(MGXS):
nu_fission_out = self.tallies['nu-fission-out']
# Sum out all nuclides
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
nu_fission_in = nu_fission_in.summation(nuclides=nuclides)
nu_fission_out = nu_fission_out.summation(nuclides=nuclides)
@ -4616,7 +4613,7 @@ class Chi(MGXS):
# Get chi for all nuclides in the domain
elif nuclides == 'all':
nuclides = self.get_all_nuclides()
nuclides = self.get_nuclides()
xs = self.xs_tally.get_values(filters=filters,
filter_bins=filter_bins,
nuclides=nuclides, value=value)

View file

@ -689,13 +689,14 @@ def get_openmc_cell(opencg_cell):
translation = np.asarray(opencg_cell.translation, dtype=np.float64)
openmc_cell.translation = translation
surfaces = []
operators = []
for surface, halfspace in opencg_cell.surfaces.values():
surfaces.append(get_openmc_surface(surface))
operators.append(operator.neg if halfspace == -1 else operator.pos)
openmc_cell.region = openmc.Intersection(
*[op(s) for op, s in zip(operators, surfaces)])
if opencg_cell.surfaces:
surfaces = []
operators = []
for surface, halfspace in opencg_cell.surfaces.values():
surfaces.append(get_openmc_surface(surface))
operators.append(operator.neg if halfspace == -1 else operator.pos)
openmc_cell.region = openmc.Intersection(
*[op(s) for op, s in zip(operators, surfaces)])
# Add the OpenMC Cell to the global collection of all OpenMC Cells
OPENMC_CELLS[cell_id] = openmc_cell

View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections import Iterable, MutableSequence
from numbers import Real, Integral
import warnings
from xml.etree import ElementTree as ET
@ -7,10 +7,8 @@ import sys
import numpy as np
from openmc.clean_xml import clean_xml_indentation
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than, check_less_than)
from openmc import Nuclide
from openmc.source import Source
import openmc.checkvalue as cv
from openmc import Nuclide, VolumeCalculation, Source
if sys.version_info[0] >= 3:
basestring = str
@ -137,6 +135,8 @@ class Settings(object):
resonance cross sections.
resonance_scattering : ResonanceScattering or iterable of ResonanceScattering
The elastic scattering model to use for resonant isotopes
volume_calculations : VolumeCalculation or iterable of VolumeCalculation
Stochastic volume calculation specifications
"""
@ -155,7 +155,7 @@ class Settings(object):
self._max_order = None
# Source subelement
self._source = None
self._source = cv.CheckedList(Source, 'source distributions')
self._confidence_intervals = None
self._cross_sections = None
@ -216,10 +216,12 @@ class Settings(object):
self._settings_file = ET.Element("settings")
self._run_mode_subelement = None
self._source_element = None
self._multipole_active = None
self._resonance_scattering = None
self._resonance_scattering = cv.CheckedList(
ResonanceScattering, 'resonance scattering models')
self._volume_calculations = cv.CheckedList(
VolumeCalculation, 'volume calculations')
@property
def run_mode(self):
@ -421,6 +423,10 @@ class Settings(object):
def resonance_scattering(self):
return self._resonance_scattering
@property
def volume_calculations(self):
return self._volume_calculations
@run_mode.setter
def run_mode(self, run_mode):
if run_mode not in ['eigenvalue', 'fixed source']:
@ -431,26 +437,26 @@ class Settings(object):
@batches.setter
def batches(self, batches):
check_type('batches', batches, Integral)
check_greater_than('batches', batches, 0)
cv.check_type('batches', batches, Integral)
cv.check_greater_than('batches', batches, 0)
self._batches = batches
@generations_per_batch.setter
def generations_per_batch(self, generations_per_batch):
check_type('generations per patch', generations_per_batch, Integral)
check_greater_than('generations per batch', generations_per_batch, 0)
cv.check_type('generations per patch', generations_per_batch, Integral)
cv.check_greater_than('generations per batch', generations_per_batch, 0)
self._generations_per_batch = generations_per_batch
@inactive.setter
def inactive(self, inactive):
check_type('inactive batches', inactive, Integral)
check_greater_than('inactive batches', inactive, 0, True)
cv.check_type('inactive batches', inactive, Integral)
cv.check_greater_than('inactive batches', inactive, 0, True)
self._inactive = inactive
@particles.setter
def particles(self, particles):
check_type('particles', particles, Integral)
check_greater_than('particles', particles, 0)
cv.check_type('particles', particles, Integral)
cv.check_greater_than('particles', particles, 0)
self._particles = particles
@keff_trigger.setter
@ -484,23 +490,21 @@ class Settings(object):
@energy_mode.setter
def energy_mode(self, energy_mode):
check_value('energy mode', energy_mode,
cv.check_value('energy mode', energy_mode,
['continuous-energy', 'multi-group'])
self._energy_mode = energy_mode
@max_order.setter
def max_order(self, max_order):
check_type('maximum scattering order', max_order, Integral)
check_greater_than('maximum scattering order', max_order, 0, True)
cv.check_type('maximum scattering order', max_order, Integral)
cv.check_greater_than('maximum scattering order', max_order, 0, True)
self._max_order = max_order
@source.setter
def source(self, source):
if isinstance(source, Source):
self._source = [source,]
else:
check_type('source distribution', source, Iterable, Source)
self._source = source
if not isinstance(source, MutableSequence):
source = [source]
self._source = cv.CheckedList(Source, 'source distributions', source)
@output.setter
def output(self, output):
@ -525,197 +529,197 @@ class Settings(object):
@output_path.setter
def output_path(self, output_path):
check_type('output path', output_path, basestring)
cv.check_type('output path', output_path, basestring)
self._output_path = output_path
@verbosity.setter
def verbosity(self, verbosity):
check_type('verbosity', verbosity, Integral)
check_greater_than('verbosity', verbosity, 1, True)
check_less_than('verbosity', verbosity, 10, True)
cv.check_type('verbosity', verbosity, Integral)
cv.check_greater_than('verbosity', verbosity, 1, True)
cv.check_less_than('verbosity', verbosity, 10, True)
self._verbosity = verbosity
@statepoint_batches.setter
def statepoint_batches(self, batches):
check_type('statepoint batches', batches, Iterable, Integral)
cv.check_type('statepoint batches', batches, Iterable, Integral)
for batch in batches:
check_greater_than('statepoint batch', batch, 0)
cv.check_greater_than('statepoint batch', batch, 0)
self._statepoint_batches = batches
@statepoint_interval.setter
def statepoint_interval(self, interval):
check_type('statepoint interval', interval, Integral)
cv.check_type('statepoint interval', interval, Integral)
self._statepoint_interval = interval
@sourcepoint_batches.setter
def sourcepoint_batches(self, batches):
check_type('sourcepoint batches', batches, Iterable, Integral)
cv.check_type('sourcepoint batches', batches, Iterable, Integral)
for batch in batches:
check_greater_than('sourcepoint batch', batch, 0)
cv.check_greater_than('sourcepoint batch', batch, 0)
self._sourcepoint_batches = batches
@sourcepoint_interval.setter
def sourcepoint_interval(self, interval):
check_type('sourcepoint interval', interval, Integral)
cv.check_type('sourcepoint interval', interval, Integral)
self._sourcepoint_interval = interval
@sourcepoint_separate.setter
def sourcepoint_separate(self, source_separate):
check_type('sourcepoint separate', source_separate, bool)
cv.check_type('sourcepoint separate', source_separate, bool)
self._sourcepoint_separate = source_separate
@sourcepoint_write.setter
def sourcepoint_write(self, source_write):
check_type('sourcepoint write', source_write, bool)
cv.check_type('sourcepoint write', source_write, bool)
self._sourcepoint_write = source_write
@sourcepoint_overwrite.setter
def sourcepoint_overwrite(self, source_overwrite):
check_type('sourcepoint overwrite', source_overwrite, bool)
cv.check_type('sourcepoint overwrite', source_overwrite, bool)
self._sourcepoint_overwrite = source_overwrite
@confidence_intervals.setter
def confidence_intervals(self, confidence_intervals):
check_type('confidence interval', confidence_intervals, bool)
cv.check_type('confidence interval', confidence_intervals, bool)
self._confidence_intervals = confidence_intervals
@cross_sections.setter
def cross_sections(self, cross_sections):
check_type('cross sections', cross_sections, basestring)
cv.check_type('cross sections', cross_sections, basestring)
self._cross_sections = cross_sections
@multipole_library.setter
def multipole_library(self, multipole_library):
check_type('cross sections', multipole_library, basestring)
cv.check_type('cross sections', multipole_library, basestring)
self._multipole_library = multipole_library
@energy_grid.setter
def energy_grid(self, energy_grid):
check_value('energy grid', energy_grid,
cv.check_value('energy grid', energy_grid,
['nuclide', 'logarithm', 'material-union'])
self._energy_grid = energy_grid
@ptables.setter
def ptables(self, ptables):
check_type('probability tables', ptables, bool)
cv.check_type('probability tables', ptables, bool)
self._ptables = ptables
@run_cmfd.setter
def run_cmfd(self, run_cmfd):
check_type('run_cmfd', run_cmfd, bool)
cv.check_type('run_cmfd', run_cmfd, bool)
self._run_cmfd = run_cmfd
@seed.setter
def seed(self, seed):
check_type('random number generator seed', seed, Integral)
check_greater_than('random number generator seed', seed, 0)
cv.check_type('random number generator seed', seed, Integral)
cv.check_greater_than('random number generator seed', seed, 0)
self._seed = seed
@survival_biasing.setter
def survival_biasing(self, survival_biasing):
check_type('survival biasing', survival_biasing, bool)
cv.check_type('survival biasing', survival_biasing, bool)
self._survival_biasing = survival_biasing
@weight.setter
def weight(self, weight):
check_type('weight cutoff', weight, Real)
check_greater_than('weight cutoff', weight, 0.0)
cv.check_type('weight cutoff', weight, Real)
cv.check_greater_than('weight cutoff', weight, 0.0)
self._weight = weight
@weight_avg.setter
def weight_avg(self, weight_avg):
check_type('average survival weight', weight_avg, Real)
check_greater_than('average survival weight', weight_avg, 0.0)
cv.check_type('average survival weight', weight_avg, Real)
cv.check_greater_than('average survival weight', weight_avg, 0.0)
self._weight_avg = weight_avg
@entropy_dimension.setter
def entropy_dimension(self, dimension):
check_type('entropy mesh dimension', dimension, Iterable, Integral)
check_length('entropy mesh dimension', dimension, 3)
cv.check_type('entropy mesh dimension', dimension, Iterable, Integral)
cv.check_length('entropy mesh dimension', dimension, 3)
self._entropy_dimension = dimension
@entropy_lower_left.setter
def entropy_lower_left(self, lower_left):
check_type('entropy mesh lower left corner', lower_left,
cv.check_type('entropy mesh lower left corner', lower_left,
Iterable, Real)
check_length('entropy mesh lower left corner', lower_left, 3)
cv.check_length('entropy mesh lower left corner', lower_left, 3)
self._entropy_lower_left = lower_left
@entropy_upper_right.setter
def entropy_upper_right(self, upper_right):
check_type('entropy mesh upper right corner', upper_right,
cv.check_type('entropy mesh upper right corner', upper_right,
Iterable, Real)
check_length('entropy mesh upper right corner', upper_right, 3)
cv.check_length('entropy mesh upper right corner', upper_right, 3)
self._entropy_upper_right = upper_right
@trigger_active.setter
def trigger_active(self, trigger_active):
check_type('trigger active', trigger_active, bool)
cv.check_type('trigger active', trigger_active, bool)
self._trigger_active = trigger_active
@trigger_max_batches.setter
def trigger_max_batches(self, trigger_max_batches):
check_type('trigger maximum batches', trigger_max_batches, Integral)
check_greater_than('trigger maximum batches', trigger_max_batches, 0)
cv.check_type('trigger maximum batches', trigger_max_batches, Integral)
cv.check_greater_than('trigger maximum batches', trigger_max_batches, 0)
self._trigger_max_batches = trigger_max_batches
@trigger_batch_interval.setter
def trigger_batch_interval(self, trigger_batch_interval):
check_type('trigger batch interval', trigger_batch_interval, Integral)
check_greater_than('trigger batch interval', trigger_batch_interval, 0)
cv.check_type('trigger batch interval', trigger_batch_interval, Integral)
cv.check_greater_than('trigger batch interval', trigger_batch_interval, 0)
self._trigger_batch_interval = trigger_batch_interval
@no_reduce.setter
def no_reduce(self, no_reduce):
check_type('no reduction option', no_reduce, bool)
cv.check_type('no reduction option', no_reduce, bool)
self._no_reduce = no_reduce
@threads.setter
def threads(self, threads):
check_type('number of threads', threads, Integral)
check_greater_than('number of threads', threads, 0)
cv.check_type('number of threads', threads, Integral)
cv.check_greater_than('number of threads', threads, 0)
self._threads = threads
@trace.setter
def trace(self, trace):
check_type('trace', trace, Iterable, Integral)
check_length('trace', trace, 3)
check_greater_than('trace batch', trace[0], 0)
check_greater_than('trace generation', trace[1], 0)
check_greater_than('trace particle', trace[2], 0)
cv.check_type('trace', trace, Iterable, Integral)
cv.check_length('trace', trace, 3)
cv.check_greater_than('trace batch', trace[0], 0)
cv.check_greater_than('trace generation', trace[1], 0)
cv.check_greater_than('trace particle', trace[2], 0)
self._trace = trace
@track.setter
def track(self, track):
check_type('track', track, Iterable, Integral)
cv.check_type('track', track, Iterable, Integral)
if len(track) % 3 != 0:
msg = 'Unable to set the track to "{0}" since its length is ' \
'not a multiple of 3'.format(track)
raise ValueError(msg)
for t in zip(track[::3], track[1::3], track[2::3]):
check_greater_than('track batch', t[0], 0)
check_greater_than('track generation', t[0], 0)
check_greater_than('track particle', t[0], 0)
cv.check_greater_than('track batch', t[0], 0)
cv.check_greater_than('track generation', t[0], 0)
cv.check_greater_than('track particle', t[0], 0)
self._track = track
@ufs_dimension.setter
def ufs_dimension(self, dimension):
check_type('UFS mesh dimension', dimension, Iterable, Integral)
check_length('UFS mesh dimension', dimension, 3)
cv.check_type('UFS mesh dimension', dimension, Iterable, Integral)
cv.check_length('UFS mesh dimension', dimension, 3)
for dim in dimension:
check_greater_than('UFS mesh dimension', dim, 1, True)
cv.check_greater_than('UFS mesh dimension', dim, 1, True)
self._ufs_dimension = dimension
@ufs_lower_left.setter
def ufs_lower_left(self, lower_left):
check_type('UFS mesh lower left corner', lower_left, Iterable, Real)
check_length('UFS mesh lower left corner', lower_left, 3)
cv.check_type('UFS mesh lower left corner', lower_left, Iterable, Real)
cv.check_length('UFS mesh lower left corner', lower_left, 3)
self._ufs_lower_left = lower_left
@ufs_upper_right.setter
def ufs_upper_right(self, upper_right):
check_type('UFS mesh upper right corner', upper_right, Iterable, Real)
check_length('UFS mesh upper right corner', upper_right, 3)
cv.check_type('UFS mesh upper right corner', upper_right, Iterable, Real)
cv.check_length('UFS mesh upper right corner', upper_right, 3)
self._ufs_upper_right = upper_right
@dd_mesh_dimension.setter
@ -724,8 +728,8 @@ class Settings(object):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
check_type('DD mesh dimension', dimension, Iterable, Integral)
check_length('DD mesh dimension', dimension, 3)
cv.check_type('DD mesh dimension', dimension, Iterable, Integral)
cv.check_length('DD mesh dimension', dimension, 3)
self._dd_mesh_dimension = dimension
@ -735,8 +739,8 @@ class Settings(object):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
check_type('DD mesh lower left corner', lower_left, Iterable, Real)
check_length('DD mesh lower left corner', lower_left, 3)
cv.check_type('DD mesh lower left corner', lower_left, Iterable, Real)
cv.check_length('DD mesh lower left corner', lower_left, 3)
self._dd_mesh_lower_left = lower_left
@ -746,8 +750,8 @@ class Settings(object):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
check_type('DD mesh upper right corner', upper_right, Iterable, Real)
check_length('DD mesh upper right corner', upper_right, 3)
cv.check_type('DD mesh upper right corner', upper_right, Iterable, Real)
cv.check_length('DD mesh upper right corner', upper_right, 3)
self._dd_mesh_upper_right = upper_right
@ -757,7 +761,7 @@ class Settings(object):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
check_type('DD nodemap', nodemap, Iterable)
cv.check_type('DD nodemap', nodemap, Iterable)
nodemap = np.array(nodemap).flatten()
@ -782,7 +786,7 @@ class Settings(object):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
check_type('DD allow leakage', allow, bool)
cv.check_type('DD allow leakage', allow, bool)
self._dd_allow_leakage = allow
@ -793,24 +797,28 @@ class Settings(object):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
check_type('DD count interactions', interactions, bool)
cv.check_type('DD count interactions', interactions, bool)
self._dd_count_interactions = interactions
@use_windowed_multipole.setter
def use_windowed_multipole(self, active):
check_type('use_windowed_multipole', active, bool)
cv.check_type('use_windowed_multipole', active, bool)
self._multipole_active = active
@resonance_scattering.setter
def resonance_scattering(self, res):
if isinstance(res, Iterable):
check_type('resonance_scattering', res, Iterable,
ResonanceScattering)
self._resonance_scattering = res
else:
check_type('resonance_scattering', res, ResonanceScattering)
self._resonance_scattering = [res]
if not isinstance(res, MutableSequence):
res = [res]
self._resonance_scattering = cv.CheckedList(
ResonanceScattering, 'resonance scattering models', res)
@volume_calculations.setter
def volume_calculations(self, vol_calcs):
if not isinstance(vol_calcs, MutableSequence):
vol_calcs = [vol_calcs]
self._volume_calculations = cv.CheckedList(
VolumeCalculation, 'stochastic volume calculations', vol_calcs)
def _create_run_mode_subelement(self):
@ -869,9 +877,12 @@ class Settings(object):
element.text = str(self._max_order)
def _create_source_subelement(self):
if self.source is not None:
for source in self.source:
self._settings_file.append(source.to_xml())
for source in self.source:
self._settings_file.append(source.to_xml_element())
def _create_volume_calcs_subelement(self):
for calc in self.volume_calculations:
self._settings_file.append(calc.to_xml_element())
def _create_output_subelement(self):
if self._output is not None:
@ -1102,18 +1113,15 @@ class Settings(object):
"use_windowed_multipole")
element.text = str(self._multipole_active)
def _create_resonance_scattering_element(self):
if self.resonance_scattering is None:
return
element = ET.SubElement(self._settings_file, "resonance_scattering")
for r in self.resonance_scattering:
if r.nuclide.name != r.nuclide_0K.name:
raise ValueError("The nuclide and nuclide_0K attributes of "
"a ResonantScattering object must have "
"identical names.")
r.create_xml_subelement(element)
def _create_resonance_scattering_subelement(self):
if len(self.resonance_scattering) > 0:
elem = ET.SubElement(self._settings_file, 'resonance_scattering')
for r in self.resonance_scattering:
if r.nuclide.name != r.nuclide_0K.name:
raise ValueError("The nuclide and nuclide_0K attributes of "
"a ResonantScattering object must have "
"identical names.")
elem.append(r.to_xml_element())
def export_to_xml(self):
"""Create a settings.xml file that can be used for a simulation.
@ -1125,7 +1133,6 @@ class Settings(object):
self._source_subelement = None
self._trigger_subelement = None
self._run_mode_subelement = None
self._source_element = None
self._create_run_mode_subelement()
self._create_source_subelement()
@ -1153,7 +1160,8 @@ class Settings(object):
self._create_ufs_subelement()
self._create_dd_subelement()
self._create_use_multipole_subelement()
self._create_resonance_scattering_element()
self._create_resonance_scattering_subelement()
self._create_volume_calcs_subelement()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._settings_file)
@ -1216,33 +1224,41 @@ class ResonanceScattering(object):
@nuclide.setter
def nuclide(self, nuc):
check_type('nuclide', nuc, Nuclide)
cv.check_type('nuclide', nuc, Nuclide)
self._nuclide = nuc
@nuclide_0K.setter
def nuclide_0K(self, nuc):
check_type('nuclide_0K', nuc, Nuclide)
cv.check_type('nuclide_0K', nuc, Nuclide)
self._nuclide_0K = nuc
@method.setter
def method(self, m):
check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM'))
cv.check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM'))
self._method = m
@E_min.setter
def E_min(self, E):
check_type('E_min', E, Real)
check_greater_than('E_min', E, 0, True)
cv.check_type('E_min', E, Real)
cv.check_greater_than('E_min', E, 0, True)
self._E_min = E
@E_max.setter
def E_max(self, E):
check_type('E_max', E, Real)
check_greater_than('E_max', E, 0, True)
cv.check_type('E_max', E, Real)
cv.check_greater_than('E_max', E, 0, True)
self._E_max = E
def create_xml_subelement(self, xml_element):
scatterer = ET.SubElement(xml_element, "scatterer")
def to_xml_element(self):
"""Return XML representation of the resonance scattering model
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing resonance scattering model
"""
scatterer = ET.Element("scatterer")
subelement = ET.SubElement(scatterer, 'nuclide')
subelement.text = self.nuclide.name
if self.method is not None:
@ -1258,3 +1274,4 @@ class ResonanceScattering(object):
if self.E_max is not None:
subelement = ET.SubElement(scatterer, 'E_max')
subelement.text = str(self.E_max)
return scatterer

View file

@ -103,7 +103,7 @@ class Source(object):
cv.check_greater_than('source strength', strength, 0.0, True)
self._strength = strength
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the source
Returns
@ -117,9 +117,9 @@ class Source(object):
if self.file is not None:
element.set("file", self.file)
if self.space is not None:
element.append(self.space.to_xml())
element.append(self.space.to_xml_element())
if self.angle is not None:
element.append(self.angle.to_xml())
element.append(self.angle.to_xml_element())
if self.energy is not None:
element.append(self.energy.to_xml('energy'))
element.append(self.energy.to_xml_element('energy'))
return element

View file

@ -2,6 +2,7 @@ import sys
import re
import os
import warnings
import glob
import numpy as np
@ -22,8 +23,9 @@ class StatePoint(object):
filename : str
Path to file to load
autolink : bool, optional
Whether to automatically link in metadata from a summary.h5
file. Defaults to True.
Whether to automatically link in metadata from a summary.h5 file and
stochastic volume calculation results from volume_*.h5 files. Defaults
to True.
Attributes
----------
@ -143,6 +145,12 @@ class StatePoint(object):
su = openmc.Summary(path_summary)
self.link_with_summary(su)
path_volume = os.path.join(os.path.dirname(filename), 'volume_*.h5')
for path_i in glob.glob(path_volume):
if re.search(r'volume_\d+\.h5', path_i):
vol = openmc.VolumeCalculation.from_hdf5(path_i)
self.add_volume_information(vol)
def close(self):
self._f.close()
@ -501,6 +509,18 @@ class StatePoint(object):
for tally_id in self.tallies:
self.tallies[tally_id].sparse = self.sparse
def add_volume_information(self, volume_calc):
"""Add volume information to the geometry within the file
Parameters
----------
volume_calc : openmc.VolumeCalculation
Results from a stochastic volume calculation
"""
if self.summary is not None:
self.summary.add_volume_information(volume_calc)
def get_tally(self, scores=[], filters=[], nuclides=[],
name=None, id=None, estimator=None, exact_filters=False,
exact_nuclides=False, exact_scores=False):

View file

@ -50,7 +50,7 @@ class UnitSphere(object):
self._reference_uvw = uvw/np.linalg.norm(uvw)
@abstractmethod
def to_xml(self):
def to_xml_element(self):
return ''
@ -109,13 +109,21 @@ class PolarAzimuthal(UnitSphere):
cv.check_type('azimuthal angle', phi, Univariate)
self._phi = phi
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the angular distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing angular distribution data
"""
element = ET.Element('angle')
element.set("type", "mu-phi")
if self.reference_uvw is not None:
element.set("reference_uvw", ' '.join(map(str, self.reference_uvw)))
element.append(self.mu.to_xml('mu'))
element.append(self.phi.to_xml('phi'))
element.append(self.mu.to_xml_element('mu'))
element.append(self.phi.to_xml_element('phi'))
return element
@ -127,7 +135,15 @@ class Isotropic(UnitSphere):
def __init__(self):
super(Isotropic, self).__init__()
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the isotropic distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing isotropic distribution data
"""
element = ET.Element('angle')
element.set("type", "isotropic")
return element
@ -152,7 +168,15 @@ class Monodirectional(UnitSphere):
def __init__(self, reference_uvw=[1., 0., 0.]):
super(Monodirectional, self).__init__(reference_uvw)
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the monodirectional distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing monodirectional distribution data
"""
element = ET.Element('angle')
element.set("type", "monodirectional")
if self.reference_uvw is not None:
@ -174,7 +198,7 @@ class Spatial(object):
pass
@abstractmethod
def to_xml(self):
def to_xml_element(self):
return ''
@ -238,12 +262,20 @@ class CartesianIndependent(Spatial):
cv.check_type('z coordinate', z, Univariate)
self._z = z
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the spatial distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing spatial distribution data
"""
element = ET.Element('space')
element.set('type', 'cartesian')
element.append(self.x.to_xml('x'))
element.append(self.y.to_xml('y'))
element.append(self.z.to_xml('z'))
element.append(self.x.to_xml_element('x'))
element.append(self.y.to_xml_element('y'))
element.append(self.z.to_xml_element('z'))
return element
@ -308,7 +340,15 @@ class Box(Spatial):
cv.check_type('only fissionable', only_fissionable, bool)
self._only_fissionable = only_fissionable
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the box distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing box distribution data
"""
element = ET.Element('space')
if self.only_fissionable:
element.set("type", "fission")
@ -352,7 +392,15 @@ class Point(Spatial):
cv.check_length('coordinate', xyz, 3)
self._xyz = xyz
def to_xml(self):
def to_xml_element(self):
"""Return XML representation of the point distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing point distribution location
"""
element = ET.Element('space')
element.set("type", "point")
params = ET.SubElement(element, "parameters")

View file

@ -29,7 +29,7 @@ class Univariate(object):
pass
@abstractmethod
def to_xml(self, element_name):
def to_xml_element(self, element_name):
return ''
@abstractmethod
@ -92,7 +92,20 @@ class Discrete(Univariate):
cv.check_greater_than('discrete probability', pk, 0.0, True)
self._p = p
def to_xml(self, element_name):
def to_xml_element(self, element_name):
"""Return XML representation of the discrete distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing discrete distribution data
"""
element = ET.Element(element_name)
element.set("type", "discrete")
@ -153,7 +166,20 @@ class Uniform(Univariate):
t.c = [0., 1.]
return t
def to_xml(self, element_name):
def to_xml_element(self, element_name):
"""Return XML representation of the uniform distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing uniform distribution data
"""
element = ET.Element(element_name)
element.set("type", "uniform")
element.set("parameters", '{} {}'.format(self.a, self.b))
@ -196,7 +222,20 @@ class Maxwell(Univariate):
cv.check_greater_than('Maxwell temperature', theta, 0.0)
self._theta = theta
def to_xml(self, element_name):
def to_xml_element(self, element_name):
"""Return XML representation of the Maxwellian distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing Maxwellian distribution data
"""
element = ET.Element(element_name)
element.set("type", "maxwell")
element.set("parameters", str(self.theta))
@ -254,7 +293,20 @@ class Watt(Univariate):
cv.check_greater_than('Watt b', b, 0.0)
self._b = b
def to_xml(self, element_name):
def to_xml_element(self, element_name):
"""Return XML representation of the Watt distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing Watt distribution data
"""
element = ET.Element(element_name)
element.set("type", "watt")
element.set("parameters", '{} {}'.format(self.a, self.b))
@ -333,7 +385,20 @@ class Tabular(Univariate):
cv.check_value('interpolation', interpolation, _INTERPOLATION_SCHEMES)
self._interpolation = interpolation
def to_xml(self, element_name):
def to_xml_element(self, element_name):
"""Return XML representation of the tabular distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing tabular distribution data
"""
element = ET.Element(element_name)
element.set("type", "tabular")
element.set("interpolation", self.interpolation)
@ -386,7 +451,7 @@ class Legendre(Univariate):
self._legendre_polynomial = np.polynomial.legendre.Legendre(
coefficients)
def to_xml(self, element_name):
def to_xml_element(self, element_name):
raise NotImplementedError
@ -440,5 +505,5 @@ class Mixture(Univariate):
Iterable, Univariate)
self._distribution = distribution
def to_xml(self, element_name):
def to_xml_element(self, element_name):
raise NotImplementedError

View file

@ -594,6 +594,17 @@ class Summary(object):
# Add Tally to the global dictionary of all Tallies
self.tallies[tally_id] = tally
def add_volume_information(self, volume_calc):
"""Add volume information to the geometry within the summary file
Parameters
----------
volume_calc : openmc.VolumeCalculation
Results from a stochastic volume calculation
"""
self.openmc_geometry.add_volume_information(volume_calc)
def get_material_by_id(self, material_id):
"""Return a Material object given the material id

View file

@ -349,7 +349,27 @@ class Universe(object):
# Return the offset computed at all nested Universe levels
return offset
def get_all_nuclides(self):
def get_nuclides(self):
"""Returns all nuclides in the universe
Returns
-------
nuclides : list of str
List of nuclide names
"""
nuclides = []
# Append all Nuclides in each Cell in the Universe to the dictionary
for cell in self.cells.values():
for nuclide in cell.get_nuclides():
if nuclide not in nuclides:
nuclides.append(nuclide)
return nuclides
def get_nuclide_densities(self):
"""Return all nuclides contained in the universe
Returns
@ -360,13 +380,8 @@ class Universe(object):
"""
nuclides = OrderedDict()
# Append all Nuclides in each Cell in the Universe to the dictionary
for cell in self._cells.values():
nuclides.update(cell.get_all_nuclides())
return nuclides
raise NotImplementedError('Determining average nuclide densities over '
'an entire universe not yet supported.')
def get_all_cells(self):
"""Return all cells that are contained within the universe

231
openmc/volume.py Normal file
View file

@ -0,0 +1,231 @@
from collections import Iterable, Mapping
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import numpy as np
import pandas as pd
import openmc
import openmc.checkvalue as cv
class VolumeCalculation(object):
"""Stochastic volume calculation specifications and results.
Parameters
----------
domains : Iterable of openmc.Cell, openmc.Material, or openmc.Universe
Domains to find volumes of
samples : int
Number of samples used to generate volume estimates
lower_left : Iterable of float
Lower-left coordinates of bounding box used to sample points. If this
argument is not supplied, an attempt is made to automatically determine
a bounding box.
upper_right : Iterable of float
Upper-right coordinates of bounding box used to sample points. If this
argument is not supplied, an attempt is made to automatically determine
a bounding box.
Attributes
----------
ids : Iterable of int
IDs of domains to find volumes of
domain_type : {'cell', 'material', 'universe'}
Type of each domain
samples : int
Number of samples used to generate volume estimates
lower_left : Iterable of float
Lower-left coordinates of bounding box used to sample points
upper_right : Iterable of float
Upper-right coordinates of bounding box used to sample points
results : dict
Dictionary whose keys are unique IDs of domains and values are
dictionaries with calculated volumes and total number of atoms for each
nuclide present in the domain.
volumes : dict
Dictionary whose keys are unique IDs of domains and values are the
estimated volumes
atoms_dataframe : pandas.DataFrame
DataFrame showing the estimated number of atoms for each nuclide present
in each domain specified.
"""
def __init__(self, domains, samples, lower_left=None,
upper_right=None):
self._results = None
cv.check_type('domains', domains, Iterable,
(openmc.Cell, openmc.Material, openmc.Universe))
if isinstance(domains[0], openmc.Cell):
self._domain_type = 'cell'
elif isinstance(domains[0], openmc.Material):
self._domain_type = 'material'
elif isinstance(domains[0], openmc.Universe):
self._domain_type = 'universe'
self.ids = [d.id for d in domains]
self.samples = samples
if lower_left is not None:
self.lower_left = lower_left
if upper_right is None:
raise ValueError('Both lower-left and upper-right coordinates '
'should be specified')
self.upper_right = upper_right
else:
if self.domain_type == 'cell':
ll, ur = openmc.Union(*[c.region for c in domains]).bounding_box
if np.any(np.isinf(ll)) or np.any(np.isinf(ur)):
raise ValueError('Could not automatically determine bounding '
'box for stochastic volume calculation.')
else:
self.lower_left = ll
self.upper_right = ur
else:
raise ValueError('Could not automatically determine bounding box '
'for stochastic volume calculation.')
@property
def ids(self):
return self._ids
@property
def samples(self):
return self._samples
@property
def lower_left(self):
return self._lower_left
@property
def upper_right(self):
return self._upper_right
@property
def results(self):
return self._results
@property
def domain_type(self):
return self._domain_type
@property
def volumes(self):
return {uid: results['volume'] for uid, results in self.results.items()}
@property
def atoms_dataframe(self):
items = []
columns = [self.domain_type.capitalize(), 'Nuclide', 'Atoms',
'Uncertainty']
for uid, results in self.results.items():
for name, atoms in results['atoms']:
items.append((uid, name, atoms[0], atoms[1]))
return pd.DataFrame.from_records(items, columns=columns)
@ids.setter
def ids(self, ids):
cv.check_type('domain IDs', ids, Iterable, Real)
self._ids = ids
@samples.setter
def samples(self, samples):
cv.check_type('number of samples', samples, Integral)
cv.check_greater_than('number of samples', samples, 0)
self._samples = samples
@lower_left.setter
def lower_left(self, lower_left):
name = 'lower-left bounding box coordinates',
cv.check_type(name, lower_left, Iterable, Real)
cv.check_length(name, lower_left, 3)
self._lower_left = lower_left
@upper_right.setter
def upper_right(self, upper_right):
name = 'upper-right bounding box coordinates'
cv.check_type(name, upper_right, Iterable, Real)
cv.check_length(name, upper_right, 3)
self._upper_right = upper_right
@results.setter
def results(self, results):
cv.check_type('results', results, Mapping)
self._results = results
@classmethod
def from_hdf5(cls, filename):
"""Load stochastic volume calculation results from HDF5 file.
Parameters
----------
filename : str
Path to volume.h5 file
Returns
-------
openmc.VolumeCalculation
Results of the stochastic volume calculation
"""
import h5py
with h5py.File(filename, 'r') as f:
domain_type = f.attrs['domain_type'].decode()
samples = f.attrs['samples']
lower_left = f.attrs['lower_left']
upper_right = f.attrs['upper_right']
results = {}
ids = []
for obj_name in f:
if obj_name.startswith('domain_'):
domain_id = int(obj_name[7:])
ids.append(domain_id)
group = f[obj_name]
volume = tuple(group['volume'].value)
nucnames = group['nuclides'].value
atoms = group['atoms'].value
atom_list = []
for name_i, atoms_i in zip(nucnames, atoms):
atom_list.append((name_i.decode(), tuple(atoms_i)))
results[domain_id] = {'volume': volume, 'atoms': atom_list}
# Instantiate some throw-away domains that are used by the constructor
# to assign IDs
if domain_type == 'cell':
domains = [openmc.Cell(uid) for uid in ids]
elif domain_type == 'material':
domains = [openmc.Material(uid) for uid in ids]
elif domain_type == 'universe':
domains = [openmc.Universe(uid) for uid in ids]
# Instantiate the class and assign results
vol = cls(domains, samples, lower_left, upper_right)
vol.results = results
return vol
def to_xml_element(self):
"""Return XML representation of the volume calculation
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing volume calculation data
"""
element = ET.Element("volume_calc")
dt_elem = ET.SubElement(element, "domain_type")
dt_elem.text = self.domain_type
id_elem = ET.SubElement(element, "domain_ids")
id_elem.text = ' '.join(str(uid) for uid in self.ids)
samples_elem = ET.SubElement(element, "samples")
samples_elem.text = str(self.samples)
ll_elem = ET.SubElement(element, "lower_left")
ll_elem.text = ' '.join(str(x) for x in self.lower_left)
ur_elem = ET.SubElement(element, "upper_right")
ur_elem.text = ' '.join(str(x) for x in self.upper_right)
return element

View file

@ -7,7 +7,7 @@ OpenMC Monte Carlo Particle Transport Code
The OpenMC project aims to provide a fully-featured Monte Carlo particle
transport code based on modern methods. It is a constructive solid geometry,
continuous-energy transport code that uses ACE format cross sections. The
continuous-energy transport code that uses HDF5 format cross sections. The
project started under the Computational Reactor Physics Group at MIT.
Complete documentation on the usage of OpenMC is hosted on Read the Docs at

View file

@ -64,7 +64,9 @@ if args.xml is not None:
directory = os.path.dirname(args.xml)
for ace_table in root.findall('ace_table'):
ace_libraries.append(os.path.join(directory, ace_table.attrib['path']))
path = os.path.join(directory, ace_table.attrib['path'])
if path not in ace_libraries:
ace_libraries.append(path)
elif args.xsdir is not None:
# Find 'directory' section
@ -75,8 +77,15 @@ elif args.xsdir is not None:
else:
raise IOError("Could not find 'directory' section in MCNP xsdir file")
# Handle continuation lines indicated by '+' at end of line
lines = lines[index + 1:]
continue_lines = [i for i, line in enumerate(lines)
if line.strip().endswith('+')]
for i in reversed(continue_lines):
lines[i] += lines[i].strip()[:-1] + lines.pop(i + 1)
# Create list of ACE libraries
for line in lines[index + 1:]:
for line in lines:
words = line.split()
if len(words) < 3:
continue
@ -109,14 +118,18 @@ for filename in ace_libraries:
for table in lib.tables:
if table.name.endswith('c'):
# Continuous-energy neutron data
neutron = openmc.data.IncidentNeutron.from_ace(
table, args.metastable)
print(neutron.name)
try:
neutron = openmc.data.IncidentNeutron.from_ace(
table, args.metastable)
except Exception as e:
print('Failed to convert {}: {}'.format(table.name, e))
continue
print('Converting {} (ACE) to {} (HDF5)'.format(table.name, neutron.name))
# Determine filename
outfile = os.path.join(args.destination,
neutron.name.replace('.', '_') + '.h5')
neutron.export_to_hdf5(outfile)
neutron.export_to_hdf5(outfile, 'w')
# Register with library
library.register_file(outfile)
@ -124,12 +137,12 @@ for filename in ace_libraries:
elif table.name.endswith('t'):
# Thermal scattering data
thermal = openmc.data.ThermalScattering.from_ace(table)
print(thermal.name)
print('Converting {} (ACE) to {} (HDF5)'.format(table.name, thermal.name))
# Determine filename
outfile = os.path.join(args.destination,
thermal.name.replace('.', '_') + '.h5')
thermal.export_to_hdf5(outfile)
thermal.export_to_hdf5(outfile, 'w')
# Register with library
library.register_file(outfile, 'thermal')

View file

@ -40,6 +40,12 @@ if have_setuptools:
'sparse' : ['scipy'],
'vtk': ['vtk', 'silomesh'],
'validate': ['lxml']
}})
},
# Data files
'package_data': {
'openmc.data': ['mass.mas12']
},
})
setup(**kwargs)

View file

@ -380,11 +380,12 @@ module constants
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS
integer, parameter :: N_STREAMS = 4
integer, parameter :: N_STREAMS = 5
integer, parameter :: STREAM_TRACKING = 1
integer, parameter :: STREAM_TALLIES = 2
integer, parameter :: STREAM_SOURCE = 3
integer, parameter :: STREAM_URR_PTABLE = 4
integer, parameter :: STREAM_VOLUME = 5
! ============================================================================
! MISCELLANEOUS CONSTANTS

View file

@ -17,6 +17,7 @@ module global
use tally_header, only: TallyObject, TallyResult
use trigger_header, only: KTrigger
use timer_header, only: Timer
use volume_header, only: VolumeCalculation
#ifdef MPIF08
use mpi_f08
@ -35,6 +36,8 @@ module global
type(Material), allocatable, target :: materials(:)
type(ObjectPlot), allocatable, target :: plots(:)
type(VolumeCalculation), allocatable :: volume_calcs(:)
! Size of main arrays
integer :: n_cells ! # of cells
integer :: n_universes ! # of universes

View file

@ -75,8 +75,15 @@ module hdf5_interface
module procedure read_attribute_string
end interface read_attribute
interface write_attribute
module procedure write_attribute_double
module procedure write_attribute_double_1D
module procedure write_attribute_integer
end interface write_attribute
public :: write_dataset
public :: read_dataset
public :: write_attribute
public :: read_attribute
public :: file_create
public :: file_open
@ -2059,6 +2066,25 @@ contains
call h5aclose_f(attr_id, hdf5_err)
end subroutine read_attribute_double
subroutine write_attribute_double(obj_id, name, buffer)
integer(HID_T), intent(in) :: obj_id
character(*), intent(in) :: name
real(8), intent(in), target :: buffer
integer :: hdf5_err
integer(HID_T) :: dspace_id
integer(HID_T) :: attr_id
type(C_PTR) :: f_ptr
call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err)
call h5acreate_f(obj_id, trim(name), H5T_NATIVE_DOUBLE, dspace_id, &
attr_id, hdf5_err)
f_ptr = c_loc(buffer)
call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err)
call h5aclose_f(attr_id, hdf5_err)
call h5sclose_f(dspace_id, hdf5_err)
end subroutine write_attribute_double
subroutine read_attribute_double_1D(buffer, obj_id, name)
real(8), target, allocatable, intent(inout) :: buffer(:)
integer(HID_T), intent(in) :: obj_id
@ -2097,6 +2123,37 @@ contains
call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err)
end subroutine read_attribute_double_1D_explicit
subroutine write_attribute_double_1D(obj_id, name, buffer)
integer(HID_T), intent(in) :: obj_id
character(*), intent(in) :: name
real(8), target, intent(in) :: buffer(:)
integer(HSIZE_T) :: dims(1)
dims(:) = shape(buffer)
call write_attribute_double_1D_explicit(obj_id, dims, name, buffer)
end subroutine write_attribute_double_1D
subroutine write_attribute_double_1D_explicit(obj_id, dims, name, buffer)
integer(HID_T), intent(in) :: obj_id
integer(HSIZE_T), intent(in) :: dims(1)
character(*), intent(in) :: name
real(8), target, intent(in) :: buffer(dims(1))
integer :: hdf5_err
integer(HID_T) :: dspace_id
integer(HID_T) :: attr_id
type(C_PTR) :: f_ptr
call h5screate_simple_f(1, dims, dspace_id, hdf5_err)
call h5acreate_f(obj_id, trim(name), H5T_NATIVE_DOUBLE, dspace_id, &
attr_id, hdf5_err)
f_ptr = c_loc(buffer)
call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err)
call h5aclose_f(attr_id, hdf5_err)
call h5sclose_f(dspace_id, hdf5_err)
end subroutine write_attribute_double_1D_explicit
subroutine read_attribute_double_2D(buffer, obj_id, name)
real(8), target, allocatable, intent(inout) :: buffer(:,:)
integer(HID_T), intent(in) :: obj_id
@ -2150,6 +2207,25 @@ contains
call h5aclose_f(attr_id, hdf5_err)
end subroutine read_attribute_integer
subroutine write_attribute_integer(obj_id, name, buffer)
integer(HID_T), intent(in) :: obj_id
character(*), intent(in) :: name
integer, intent(in), target :: buffer
integer :: hdf5_err
integer(HID_T) :: dspace_id
integer(HID_T) :: attr_id
type(C_PTR) :: f_ptr
call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err)
call h5acreate_f(obj_id, trim(name), H5T_NATIVE_INTEGER, dspace_id, &
attr_id, hdf5_err)
f_ptr = c_loc(buffer)
call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err)
call h5aclose_f(attr_id, hdf5_err)
call h5sclose_f(dspace_id, hdf5_err)
end subroutine write_attribute_integer
subroutine read_attribute_integer_1D(buffer, obj_id, name)
integer, target, allocatable, intent(inout) :: buffer(:)
integer(HID_T), intent(in) :: obj_id

View file

@ -88,8 +88,10 @@ contains
type(Node), pointer :: node_scatterer => null()
type(Node), pointer :: node_trigger => null()
type(Node), pointer :: node_keff_trigger => null()
type(Node), pointer :: node_vol => null()
type(NodeList), pointer :: node_scat_list => null()
type(NodeList), pointer :: node_source_list => null()
type(NodeList), pointer :: node_vol_list => null()
! Check if settings.xml exists
filename = trim(path_input) // "settings.xml"
@ -1120,6 +1122,14 @@ contains
end select
end if
call get_node_list(doc, "volume_calc", node_vol_list)
n = get_list_size(node_vol_list)
allocate(volume_calcs(n))
do i = 1, n
call get_list_item(node_vol_list, i, node_vol)
call volume_calcs(i) % from_xml(node_vol)
end do
! Close settings XML file
call close_xmldoc(doc)

View file

@ -29,7 +29,6 @@ module random_lcg
public :: set_particle_seed
public :: advance_prn_seed
public :: prn_set_stream
public :: STREAM_TRACKING, STREAM_TALLIES
contains

View file

@ -142,6 +142,19 @@ element settings {
element verbosity { xsd:positiveInteger }? &
element volume_calc {
(element domain_type { xsd:string } |
attribute domain_type { xsd:string }) &
(element domain_ids { list { xsd:integer+ } } |
attribute domain_ids { list { xsd:integer+ } }) &
(element samples { xsd:positiveInteger } |
attribute samples { xsd:positiveInteger }) &
(element lower_left { list { xsd:double+ } } |
attribute lower_left { list { xsd:double+ } }) &
(element upper_right { list { xsd:double+ } } |
attribute upper_right { list { xsd:double+ } })
}+ &
element uniform_fs{
(element dimension { list { xsd:positiveInteger+ } } |
attribute dimension { list { xsd:positiveInteger+ } }) &

View file

@ -625,6 +625,76 @@
<data type="positiveInteger"/>
</element>
</optional>
<oneOrMore>
<element name="volume_calc">
<interleave>
<choice>
<element name="domain_type">
<data type="string"/>
</element>
<attribute name="domain_type">
<data type="string"/>
</attribute>
</choice>
<choice>
<element name="domain_ids">
<list>
<oneOrMore>
<data type="integer"/>
</oneOrMore>
</list>
</element>
<attribute name="domain_ids">
<list>
<oneOrMore>
<data type="integer"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="samples">
<data type="positiveInteger"/>
</element>
<attribute name="samples">
<data type="positiveInteger"/>
</attribute>
</choice>
<choice>
<element name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</interleave>
</element>
</oneOrMore>
<optional>
<element name="uniform_fs">
<interleave>

View file

@ -24,6 +24,7 @@ module simulation
reset_result
use trigger, only: check_triggers
use tracking, only: transport
use volume_calc, only: run_volume_calculations
implicit none
private
@ -42,6 +43,9 @@ contains
type(Particle) :: p
integer(8) :: i_work
! Volume calculations
if (size(volume_calcs) > 0) call run_volume_calculations()
if (.not. restart_run) call initialize_source()
! Display header

View file

@ -114,6 +114,11 @@ contains
! Since integer is trivially destructible, we only need to set size to zero
! and can leave capacity as is
this%size_ = 0
if (allocated(this % data)) then
this%capacity_ = size(this % data)
else
this%capacity_ = 0
end if
end subroutine clear_int
subroutine initialize_fill_int(this, n, val)
@ -249,6 +254,11 @@ contains
! Since real is trivially destructible, we only need to set size to zero and
! can leave capacity as is
this%size_ = 0
if (allocated(this % data)) then
this%capacity_ = size(this % data)
else
this%capacity_ = 0
end if
end subroutine clear_real
subroutine initialize_fill_real(this, n, val)
@ -384,6 +394,11 @@ contains
! Since char is trivially destructible, we only need to set size to zero and
! can leave capacity as is
this%size_ = 0
if (allocated(this % data)) then
this%capacity_ = size(this % data)
else
this%capacity_ = 0
end if
end subroutine clear_char
subroutine initialize_fill_char(this, n, val)

475
src/volume_calc.F90 Normal file
View file

@ -0,0 +1,475 @@
module volume_calc
use hdf5, only: HID_T
#ifdef _OPENMP
use omp_lib
#endif
use constants
use geometry, only: find_cell
use global
use hdf5_interface, only: file_create, file_close, write_attribute, &
create_group, close_group, write_dataset, write_attribute_string
use output, only: write_message, header
use message_passing
use particle_header, only: Particle
use random_lcg, only: prn, prn_set_stream, set_particle_seed
use stl_vector, only: VectorInt, VectorReal
use timer_header, only: Timer
use volume_header
implicit none
private
public :: run_volume_calculations
contains
!===============================================================================
! RUN_VOLUME_CALCULATIONS runs each of the stochastic volume calculations that
! the user has specified and writes results to HDF5 files
!===============================================================================
subroutine run_volume_calculations()
integer :: i, j
integer :: n
real(8), allocatable :: volume(:,:) ! volume mean/stdev in each domain
character(10) :: domain_type
character(MAX_FILE_LEN) :: filename ! filename for HDF5 file
type(Timer) :: time_volume ! timer for volume calculation
type(VectorInt), allocatable :: nuclide_vec(:) ! indices in nuclides array
type(VectorReal), allocatable :: atoms_vec(:) ! total # of atoms of each nuclide
type(VectorReal), allocatable :: uncertainty_vec(:) ! uncertainty of total # of atoms
if (master) then
call header("STOCHASTIC VOLUME CALCULATION", level=1)
call time_volume % start()
end if
do i = 1, size(volume_calcs)
n = size(volume_calcs(i) % domain_id)
allocate(nuclide_vec(n))
allocate(atoms_vec(n), uncertainty_vec(n))
allocate(volume(2,n))
if (master) then
call write_message("Running volume calculation " // trim(to_str(i)) &
// "...")
end if
call get_volume(volume_calcs(i), volume, nuclide_vec, atoms_vec, &
uncertainty_vec)
if (master) then
select case (volume_calcs(i) % domain_type)
case (FILTER_CELL)
domain_type = ' Cell'
case (FILTER_MATERIAL)
domain_type = ' Material'
case (FILTER_UNIVERSE)
domain_type = ' Universe'
end select
! Display domain volumes
do j = 1, size(volume_calcs(i) % domain_id)
call write_message(trim(domain_type) // " " // trim(to_str(&
volume_calcs(i) % domain_id(j))) // ": " // trim(to_str(&
volume(1,j))) // " +/- " // trim(to_str(volume(2,j))) // " cm^3")
end do
call write_message("")
filename = trim(path_output) // 'volume_' // trim(to_str(i)) // '.h5'
call write_volume(volume_calcs(i), filename, volume, nuclide_vec, &
atoms_vec, uncertainty_vec)
end if
deallocate(nuclide_vec, atoms_vec, uncertainty_vec, volume)
end do
! Show elapsed time
if (master) then
call time_volume % stop()
call write_message("Elapsed time: " // trim(to_str(time_volume % &
get_value())) // " s")
end if
end subroutine run_volume_calculations
!===============================================================================
! GET_VOLUME stochastically determines the volume of a set of domains along with
! the average number densities of nuclides within the domain
!===============================================================================
subroutine get_volume(this, volume, nuclide_vec, atoms_vec, uncertainty_vec)
type(VolumeCalculation), intent(in) :: this
real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each domain
type(VectorInt), intent(out) :: nuclide_vec(:) ! indices in nuclides array
type(VectorReal), intent(out) :: atoms_vec(:) ! total # of atoms of each nuclide
type(VectorReal), intent(out) :: uncertainty_vec(:) ! uncertainty of total # of atoms
! Variables that are private to each thread
integer(8) :: i
integer :: j, k
integer :: i_domain ! index in domain_id array
integer :: i_material ! index in materials array
integer :: level ! local coordinate level
integer :: n_mat(size(this % domain_id)) ! Number of materials for each domain
integer, allocatable :: indices(:,:) ! List of material indices for each domain
integer, allocatable :: hits(:,:) ! Number of hits for each material in each domain
logical :: found_cell
type(Particle) :: p
! Shared variables
integer :: i_start, i_end ! Starting/ending sample for each process
type(VectorInt) :: master_indices(size(this % domain_id))
type(VectorInt) :: master_hits(size(this % domain_id))
! Variables used outside of parallel region
integer :: i_nuclide ! index in nuclides array
integer :: total_hits ! total hits for a single domain (summed over materials)
integer :: min_samples ! minimum number of samples per process
integer :: remainder ! leftover samples from uneven divide
#ifdef MPI
integer :: m ! index over materials
integer :: n ! number of materials
integer, allocatable :: data(:) ! array used to send number of hits
#endif
real(8) :: f ! fraction of hits
real(8) :: var_f ! variance of fraction of hits
real(8) :: volume_sample ! total volume of sampled region
real(8) :: atoms(2, size(nuclides))
! Divide work over MPI processes
min_samples = this % samples / n_procs
remainder = mod(this % samples, n_procs)
if (rank < remainder) then
i_start = (min_samples + 1)*rank
i_end = i_start + min_samples
else
i_start = (min_samples + 1)*remainder + (rank - remainder)*min_samples
i_end = i_start + min_samples - 1
end if
call p % initialize()
!$omp parallel private(i, j, k, i_domain, i_material, level, found_cell, &
!$omp& indices, hits, n_mat) firstprivate(p)
! Create space for material indices and number of hits for each
allocate(indices(size(this % domain_id), 8))
allocate(hits(size(this % domain_id), 8))
n_mat(:) = 0
call prn_set_stream(STREAM_VOLUME)
! ==========================================================================
! SAMPLES LOCATIONS AND COUNT HITS
!$omp do
SAMPLE_LOOP: do i = i_start, i_end
call set_particle_seed(i)
p % n_coord = 1
p % coord(1) % xyz(1) = this % lower_left(1) + prn()*(&
this % upper_right(1) - this % lower_left(1))
p % coord(1) % xyz(2) = this % lower_left(2) + prn()*(&
this % upper_right(2) - this % lower_left(2))
p % coord(1) % xyz(3) = this % lower_left(3) + prn()*(&
this % upper_right(3) - this % lower_left(3))
p % coord(1) % uvw(:) = [HALF, HALF, HALF]
! If this location is not in the geometry at all, move on to the next
! block
call find_cell(p, found_cell)
if (.not. found_cell) cycle
if (this % domain_type == FILTER_MATERIAL) then
i_material = p % material
do i_domain = 1, size(this % domain_id)
if (i_material == materials(i_domain) % id) then
call check_hit(i_domain, i_material, indices, hits, n_mat)
end if
end do
elseif (this % domain_type == FILTER_CELL) THEN
do level = 1, p % n_coord
do i_domain = 1, size(this % domain_id)
if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then
i_material = p % material
call check_hit(i_domain, i_material, indices, hits, n_mat)
end if
end do
end do
elseif (this % domain_type == FILTER_UNIVERSE) then
do level = 1, p % n_coord
do i_domain = 1, size(this % domain_id)
if (universes(p % coord(level) % universe) % id == &
this % domain_id(i_domain)) then
i_material = p % material
call check_hit(i_domain, i_material, indices, hits, n_mat)
end if
end do
end do
end if
end do SAMPLE_LOOP
!$omp end do
! ==========================================================================
! REDUCE HITS ONTO MASTER THREAD
! At this point, each thread has its own pair of index/hits lists and we now
! need to reduce them. OpenMP is not nearly smart enough to do this on its
! own, so we have to manually reduce them.
#ifdef _OPENMP
!$omp do ordered schedule(static)
THREAD_LOOP: do i = 1, omp_get_num_threads()
!$omp ordered
do i_domain = 1, size(this % domain_id)
INDEX_LOOP: do j = 1, n_mat(i_domain)
! Check if this material has been added to the master list and if so,
! accumulate the number of hits
do k = 1, master_indices(i_domain) % size()
if (indices(i_domain, j) == master_indices(i_domain) % data(k)) then
master_hits(i_domain) % data(k) = &
master_hits(i_domain) % data(k) + hits(i_domain, j)
cycle INDEX_LOOP
end if
end do
! If we made it here, this means the material hasn't yet been added to
! the master list, so add an entry to both the master indices and master
! hits lists
call master_indices(i_domain) % push_back(indices(i_domain, j))
call master_hits(i_domain) % push_back(hits(i_domain, j))
end do INDEX_LOOP
end do
!$omp end ordered
end do THREAD_LOOP
!$omp end do
#else
do i_domain = 1, size(this % domain_id)
do j = 1, n_mat(i_domain)
call master_indices(i_domain) % push_back(indices(i_domain, j))
call master_hits(i_domain) % push_back(hits(i_domain, j))
end do
end do
#endif
call prn_set_stream(STREAM_TRACKING)
!$omp end parallel
! ==========================================================================
! REDUCE HITS ONTO MASTER PROCESS
volume_sample = product(this % upper_right - this % lower_left)
do i_domain = 1, size(this % domain_id)
atoms(:, :) = ZERO
total_hits = 0
if (master) then
#ifdef MPI
do j = 1, n_procs - 1
call MPI_RECV(n, 1, MPI_INTEGER, j, 0, MPI_COMM_WORLD, &
MPI_STATUS_IGNORE, mpi_err)
allocate(data(2*n))
call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, MPI_COMM_WORLD, &
MPI_STATUS_IGNORE, mpi_err)
do k = 0, n - 1
do m = 1, master_indices(i_domain) % size()
if (data(2*k + 1) == master_indices(i_domain) % data(m)) then
master_hits(i_domain) % data(m) = master_hits(i_domain) % data(m) + &
data(2*k + 2)
end if
end do
end do
deallocate(data)
end do
#endif
do j = 1, master_indices(i_domain) % size()
total_hits = total_hits + master_hits(i_domain) % data(j)
f = real(master_hits(i_domain) % data(j), 8) / this % samples
var_f = f*(ONE - f) / this % samples
i_material = master_indices(i_domain) % data(j)
if (i_material == MATERIAL_VOID) cycle
associate (mat => materials(i_material))
do k = 1, size(mat % nuclide)
! Accumulate nuclide density
i_nuclide = mat % nuclide(k)
atoms(1, i_nuclide) = atoms(1, i_nuclide) + &
mat % atom_density(k) * f
atoms(2, i_nuclide) = atoms(2, i_nuclide) + &
mat % atom_density(k)**2 * var_f
end do
end associate
end do
! Determine volume
volume(1, i_domain) = real(total_hits, 8) / this % samples * volume_sample
volume(2, i_domain) = sqrt(volume(1, i_domain) * (volume_sample - &
volume(1, i_domain)) / this % samples)
! Determine total number of atoms. At this point, we have values in
! atoms/b-cm. To get to atoms we multiple by 10^24 V.
do j = 1, size(atoms, 2)
atoms(1, j) = 1.0e24_8 * volume_sample * atoms(1, j)
atoms(2, j) = 1.0e24_8 * volume_sample * sqrt(atoms(2, j))
end do
! Convert full arrays to vectors
do j = 1, size(nuclides)
if (atoms(1, j) > ZERO) then
call nuclide_vec(i_domain) % push_back(j)
call atoms_vec(i_domain) % push_back(atoms(1, j))
call uncertainty_vec(i_domain) % push_back(atoms(2, j))
end if
end do
else
#ifdef MPI
n = master_indices(i_domain) % size()
allocate(data(2*n))
do k = 0, n - 1
data(2*k + 1) = master_indices(i_domain) % data(k + 1)
data(2*k + 2) = master_hits(i_domain) % data(k + 1)
end do
call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, mpi_err)
call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, MPI_COMM_WORLD, mpi_err)
deallocate(data)
#endif
end if
end do
contains
!===========================================================================
! CHECK_HIT is an internal subroutine that checks for whether a material has
! already been hit for a given domain. If not, it increases the list size by
! one (taking care of re-allocation if needed).
!===========================================================================
subroutine check_hit(i_domain, i_material, indices, hits, n_mat)
integer :: i_domain
integer :: i_material
integer, allocatable :: indices(:,:)
integer, allocatable :: hits(:,:)
integer :: n_mat(:)
integer, allocatable :: temp(:,:)
logical :: already_hit
integer :: j, k, nm
! Check if we've already had a hit in this material and if so,
! simply add one
already_hit = .false.
nm = n_mat(i_domain)
do j = 1, nm
if (indices(i_domain, j) == i_material) then
hits(i_domain, j) = hits(i_domain, j) + 1
already_hit = .true.
end if
end do
if (.not. already_hit) then
! If we make it here, that means we haven't yet had a hit in this
! material. First check if the indices and hits arrays are large enough
! and if not, double them.
if (nm == size(indices, 2)) then
k = 2*size(indices, 2)
allocate(temp(size(this % domain_id), k))
temp(:, 1:nm) = indices(:, 1:nm)
call move_alloc(FROM=temp, TO=indices)
allocate(temp(size(this % domain_id), k))
temp(:, 1:nm) = hits(:, 1:nm)
call move_alloc(FROM=temp, TO=indices)
end if
! Add an entry to both the indices list and the hits list
n_mat(i_domain) = n_mat(i_domain) + 1
indices(i_domain, n_mat(i_domain)) = i_material
hits(i_domain, n_mat(i_domain)) = 1
end if
end subroutine check_hit
end subroutine get_volume
!===============================================================================
! WRITE_VOLUME writes the results of a single stochastic volume calculation to
! an HDF5 file
!===============================================================================
subroutine write_volume(this, filename, volume, nuclide_vec, atoms_vec, &
uncertainty_vec)
type(VolumeCalculation), intent(in) :: this
character(*), intent(in) :: filename ! filename for HDF5 file
real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each domain
type(VectorInt), intent(in) :: nuclide_vec(:) ! indices in nuclides array
type(VectorReal), intent(in) :: atoms_vec(:) ! total # of atoms of each nuclide
type(VectorReal), intent(in) :: uncertainty_vec(:) ! uncertainty of total # of atoms
integer :: i, j
integer :: n
integer(HID_T) :: file_id
integer(HID_T) :: group_id
real(8), allocatable :: atom_data(:,:) ! mean/stdev of total # of atoms for
! each nuclide
character(MAX_WORD_LEN), allocatable :: nucnames(:) ! names of nuclides
! Create HDF5 file
file_id = file_create(filename)
! Write basic metadata
select case (this % domain_type)
case (FILTER_CELL)
call write_attribute_string(file_id, ".", "domain_type", "cell")
case (FILTER_MATERIAL)
call write_attribute_string(file_id, ".", "domain_type", "material")
case (FILTER_UNIVERSE)
call write_attribute_string(file_id, ".", "domain_type", "universe")
end select
call write_attribute(file_id, "samples", this % samples)
call write_attribute(file_id, "lower_left", this % lower_left)
call write_attribute(file_id, "upper_right", this % upper_right)
do i = 1, size(this % domain_id)
group_id = create_group(file_id, "domain_" // trim(to_str(&
this % domain_id(i))))
! Write volume for domain
call write_dataset(group_id, "volume", volume(:, i))
! Create array of nuclide names from the vector
n = nuclide_vec(i) % size()
if (n > 0) then
allocate(nucnames(n))
do j = 1, n
nucnames(j) = nuclides(nuclide_vec(i) % data(j)) % name
end do
! Create array of total # of atoms with uncertainty for each nuclide
allocate(atom_data(2, n))
atom_data(1, :) = atoms_vec(i) % data(1:n)
atom_data(2, :) = uncertainty_vec(i) % data(1:n)
! Write results
call write_dataset(group_id, "nuclides", nucnames)
call write_dataset(group_id, "atoms", atom_data)
deallocate(nucnames)
deallocate(atom_data)
end if
call close_group(group_id)
end do
call file_close(file_id)
end subroutine write_volume
end module volume_calc

59
src/volume_header.F90 Normal file
View file

@ -0,0 +1,59 @@
module volume_header
use constants, only: FILTER_CELL, FILTER_MATERIAL, FILTER_UNIVERSE
use error, only: fatal_error
use xml_interface
implicit none
type VolumeCalculation
integer :: domain_type
integer, allocatable :: domain_id(:)
real(8) :: lower_left(3)
real(8) :: upper_right(3)
integer :: samples
contains
procedure :: from_xml => volume_from_xml
end type VolumeCalculation
contains
subroutine volume_from_xml(this, node_vol)
class(VolumeCalculation), intent(out) :: this
type(Node), pointer :: node_vol
integer :: num_domains
character(10) :: temp_str
! Check domain type
call get_node_value(node_vol, "domain_type", temp_str)
select case (temp_str)
case ('cell')
this % domain_type = FILTER_CELL
case ('material')
this % domain_type = FILTER_MATERIAL
case ('universe')
this % domain_type = FILTER_UNIVERSE
case default
call fatal_error("Unrecognized domain type for stochastic volume &
&calculation: " // trim(temp_str))
end select
! Read cell IDs
if (check_for_node(node_vol, "domain_ids")) then
num_domains = get_arraysize_integer(node_vol, "domain_ids")
else
call fatal_error("Must specify at least one cell for a volume calculation")
end if
allocate(this % domain_id(num_domains))
call get_node_array(node_vol, "domain_ids", this % domain_id)
! Read lower-left and upper-right bounding coordinates
call get_node_array(node_vol, "lower_left", this % lower_left)
call get_node_array(node_vol, "upper_right", this % upper_right)
! Read number of samples
call get_node_value(node_vol, "samples", this % samples)
end subroutine volume_from_xml
end module volume_header

View file

@ -0,0 +1 @@
102569289552d021b6803f404a0c17a9c17a40578fdba43a6ba08b77a731e0368fffa6a8a7abd48555167cb9997c6dba9ec5044c8593b12056957b7e3ec44ed0

View file

@ -0,0 +1,31 @@
k-combined: 4.165451e-02 3.582531e-04
Volume calculation 0
Domain 1: 31.4693 +/- 0.0721 cm^3
Domain 2: 2.0933 +/- 0.0310 cm^3
Domain 3: 2.0486 +/- 0.0307 cm^3
Cell Nuclide Atoms Uncertainty
0 1 U235.71c 3.481769e+23 7.979991e+20
1 1 Mo99.71c 3.481769e+22 7.979991e+19
2 2 H1.71c 1.399770e+23 2.072914e+21
3 2 O16.71c 6.998852e+22 1.036457e+21
4 2 B10.71c 6.998852e+18 1.036457e+17
5 3 H1.71c 1.369920e+23 2.051689e+21
6 3 O16.71c 6.849599e+22 1.025844e+21
7 3 B10.71c 6.849599e+18 1.025844e+17
Volume calculation 1
Domain 1: 4.1419 +/- 0.0426 cm^3
Domain 2: 31.4693 +/- 0.0721 cm^3
Material Nuclide Atoms Uncertainty
0 1 H1.71c 2.769690e+23 2.850068e+21
1 1 O16.71c 1.384845e+23 1.425034e+21
2 1 B10.71c 1.384845e+19 1.425034e+17
3 2 U235.71c 3.481769e+23 7.979991e+20
4 2 Mo99.71c 3.481769e+22 7.979991e+19
Volume calculation 2
Domain 0: 35.6112 +/- 0.0664 cm^3
Universe Nuclide Atoms Uncertainty
0 0 H1.71c 2.769690e+23 2.850068e+21
1 0 O16.71c 1.384845e+23 1.425034e+21
2 0 B10.71c 1.384845e+19 1.425034e+17
3 0 U235.71c 3.481769e+23 7.979991e+20
4 0 Mo99.71c 3.481769e+22 7.979991e+19

View file

@ -0,0 +1,89 @@
#!/usr/bin/env python
import os
import glob
import sys
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
import openmc
class VolumeTest(PyAPITestHarness):
def _build_inputs(self):
# Define materials
water = openmc.Material(1)
water.add_nuclide('H1', 2.0)
water.add_nuclide('O16', 1.0)
water.add_nuclide('B10', 0.0001)
water.add_s_alpha_beta('c_H_in_H2O', '71t')
water.set_density('g/cc', 1.0)
fuel = openmc.Material(2)
fuel.add_nuclide('U235', 1.0)
fuel.add_nuclide('Mo99', 0.1)
fuel.set_density('g/cc', 4.5)
materials = openmc.Materials((water, fuel))
materials.default_xs = '71c'
materials.export_to_xml()
cyl = openmc.ZCylinder(1, R=1.0, boundary_type='vacuum')
top_sphere = openmc.Sphere(2, z0=5., R=1., boundary_type='vacuum')
top_plane = openmc.ZPlane(3, z0=5.)
bottom_sphere = openmc.Sphere(4, z0=-5., R=1., boundary_type='vacuum')
bottom_plane = openmc.ZPlane(5, z0=-5.)
# Define geometry
inside_cyl = openmc.Cell(1, fill=fuel, region=-cyl & -top_plane & +bottom_plane)
top_hemisphere = openmc.Cell(2, fill=water, region=-top_sphere & +top_plane)
bottom_hemisphere = openmc.Cell(3, fill=water, region=-bottom_sphere & -top_plane)
root = openmc.Universe(0, cells=(inside_cyl, top_hemisphere, bottom_hemisphere))
geometry = openmc.Geometry()
geometry.root_universe = root
geometry.export_to_xml()
# Set up stochastic volume calculation
ll, ur = openmc.Union(*[c.region for c in root.cells.values()]).bounding_box
vol_calcs = [
openmc.VolumeCalculation(list(root.cells.values()), 100000),
openmc.VolumeCalculation([water, fuel], 100000, ll, ur),
openmc.VolumeCalculation([root], 100000, ll, ur)
]
# Define settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 4
settings.inactive = 0
settings.source = openmc.Source(space=openmc.stats.Box(
[-1., -1., -5.], [1., 1., 5.]))
settings.volume_calculations = vol_calcs
settings.export_to_xml()
def _get_results(self):
# Read the statepoint file.
statepoint = os.path.join(os.getcwd(), self._sp_name)
sp = openmc.StatePoint(statepoint)
# Write out k-combined.
outstr = 'k-combined: {:12.6e} {:12.6e}\n'.format(*sp.k_combined)
for i, filename in enumerate(sorted(glob.glob(os.path.join(
os.getcwd(), 'volume_*.h5')))):
outstr += 'Volume calculation {}\n'.format(i)
# Read volume calculation results
vol = openmc.VolumeCalculation.from_hdf5(filename)
# Write cell volumes and total # of atoms for each nuclide
for uid, results in sorted(vol.results.items()):
outstr += 'Domain {0}: {1[0]:.4f} +/- {1[1]:.4f} cm^3\n'.format(
uid, results['volume'])
outstr += str(vol.atoms_dataframe) + '\n'
return outstr
if __name__ == '__main__':
harness = VolumeTest('statepoint.4.h5')
harness.main()

View file

@ -137,6 +137,7 @@ class TestHarness(object):
output.append(os.path.join(os.getcwd(), 'tallies.out'))
output.append(os.path.join(os.getcwd(), 'results_test.dat'))
output.append(os.path.join(os.getcwd(), 'summary.h5'))
output += glob.glob(os.path.join(os.getcwd(), 'volume_*.h5'))
for f in output:
if os.path.exists(f):
os.remove(f)