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

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

Browse source
This commit is contained in:
Sterling Harper 2016-12-02 14:37:27 -05:00
commit 52bd883e58
30 changed files with 1642 additions and 520 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/conf.py
View file

@ -25,7 +25,7 @@ except ImportError:
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'h5py', 'pandas', 'opencg']
'h5py', 'pandas', 'uncertainties', 'opencg']
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
import numpy as np

3
docs/source/developers.rst
View file

@ -13,6 +13,7 @@ Active development of the OpenMC Monte Carlo code is currently led by:
* `Jon Walsh <https://github.com/walshjon>`_
* `Sterling Harper <https://github.com/smharper>`_
* `Will Boyd <https://github.com/wbinventor>`_
* `Samuel Shaner <https://github.com/samuelshaner>`_
* `Benoit Forget <http://web.mit.edu/nse/people/faculty/forget.html>`_
* `Kord Smith <http://web.mit.edu/nse/people/faculty/smith.html>`_
* `Andrew Siegel <http://www.mcs.anl.gov/about/people_detail.php?id=404>`_
* `Andrew Siegel <http://www.mcs.anl.gov/person/andrew-siegel>`_

2
docs/source/devguide/workflow.rst
View file

@ -238,4 +238,4 @@ from your private repository into a public fork.
.. _paid plan: https://github.com/plans
.. _Bitbucket: https://bitbucket.org
.. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html
.. _NNDC: http://http://www.nndc.bnl.gov/endf/b7.1/acefiles.html
.. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html

6
docs/source/methods/geometry.rst
View file

@ -384,7 +384,7 @@ x - x_0`, :math:`\bar{y} = y - y_0`, and :math:`\bar{z} = z - z_0`. We then have
Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain
.. math::
:label: dist-xcylinder-2
:label: dist-xcone-2
(v^2 + w^2 - R^2u^2) d^2 + 2 (\bar{y}v + \bar{z}w - R^2\bar{x}u) d +
(\bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2) = 0
@ -392,7 +392,7 @@ Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain
Defining the terms
.. math::
:label: dist-quadric-terms
:label: dist-xcone-terms
a = v^2 + w^2 - R^2u^2
@ -896,4 +896,4 @@ Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, the gradient to the surface is
.. _surfaces: http://en.wikipedia.org/wiki/Surface
.. _MCNP: http://mcnp.lanl.gov
.. _Serpent: http://montecarlo.vtt.fi
.. _Monte Carlo Performance benchmark: https://github.com/paulromano/benchmarks/tree/master/mc-performance/openmc
.. _Monte Carlo Performance benchmark: https://github.com/mit-crpg/benchmarks/tree/master/mc-performance/openmc

28
docs/source/publications.rst
View file

@ -82,6 +82,10 @@ Coupling and Multi-physics
Geometry
--------
- Derek M. Lax, "Memory efficient indexing algorithm for physical properties in
OpenMC," S. M. Thesis, Massachusetts Institute of Technology
(2015). `<http://hdl.handle.net/1721.1/97862>`_
- Derek Lax, William Boyd, Nicholas Horelik, Benoit Forget, and Kord Smith, "A
memory efficient algorithm for classifying unique regions in constructive
solid geometries," *Proc. PHYSOR*, Kyoto, Japan, Sep. 28--Oct. 3 (2014).
@ -106,6 +110,10 @@ Miscellaneous
triggers for the OpenMC Monte Carlo code," *Trans. Am. Nucl. Soc.*, **112**,
637-640 (2015).
- Kyungkwan Noh and Deokjung Lee, "Whole Core Analysis using OpenMC Monte Carlo
Code," *Trans. Kor. Nucl. Soc. Autumn Meeting*, Gyeongju, Korea,
Oct. 24-25, 2013.
- Timothy P. Burke, Brian C. Kiedrowski, and William R. Martin, "Flux and
Reaction Rate Kernel Density Estimators in OpenMC," *Trans. Am. Nucl. Soc.*,
**109**, 683-686 (2013).
@ -139,7 +147,7 @@ Doppler Broadening
- Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed multipole
for cross section Doppler broadening," *J. Comput. Phys.*, **307**, 715-727
(2016). `<http://dx.doi.org/10.1016/jcp.2015.08.013>`_
(2016). `<http://dx.doi.org/10.1016/j.jcp.2015.08.013>`_
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown,
"On-the-fly Doppler Broadening of Unresolved Resonance Region Cross Sections
@ -169,6 +177,11 @@ Nuclear Data
- Paul K. Romano and Sterling M. Harper, "Nuclear data processing capabilities
in OpenMC", *Proc. Nuclear Data*, Sep. 11-16, 2016.
- Jonathan A. Walsh, Benoit Froget, Kord S. Smith, and Forrest B. Brown,
"Neutron Cross Section Processing Methods for Improved Integral Benchmarking
of Unresolved Resonance Region Evaluations," *Eur. Phys. J. Web Conf.*
**111**, 06001 (2016). `<http://dx.doi.org/10.1051/epjconf/201611106001>`_
- Jonathan A. Walsh, Paul K. Romano, Benoit Forget, and Kord S. Smith,
"Optimizations of the energy grid search algorithm in continuous-energy Monte
Carlo particle transport codes", *Comput. Phys. Commun.*, **196**, 134-142
@ -213,10 +226,19 @@ Parallelism
memory subsystem on Monte Carlo code performance," *Proc. Joint
Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015).
- Hajime Fujita, Nan Dun, Aiman Fang, Zachary A. Rubinstein, Ziming Zheng, Kamil
Iskra, Jeff Hammonds, Anshu Dubey, Pavan Balaji, and Andrew A. Chien, "Using
Global View Resilience (GVR) to add Resilience to Exascale Applications,"
*Proc. Supercomputing*, New Orleans, Louisiana, Nov. 16--21, 2014.
- Nicholas Horelik, Benoit Forget, Kord Smith, and Andrew Siegel, "Domain
decomposition and terabyte tallies with the OpenMC Monte Carlo neutron
transport code," *Proc. PHYSOR*, Kyoto Japan, Sep. 28--Oct. 3 (2014).
- John R. Tramm, Andrew R. Siegel, Tanzima Islam, and Martin Schulz, "XSBench --
the development and verification of a performance abstraction for Monte Carlo
reactor analysis," *Proc. PHYSOR*, Kyoto, Japan, Sep 28--Oct. 3, 2014.
- Nicholas Horelik, Andrew Siegel, Benoit Forget, and Kord Smith, "Monte Carlo
domain decomposition for robust nuclear reactor analysis," *Parallel Comput.*,
**40**, 646--660 (2014). `<http://dx.doi.org/10.1016/j.parco.2014.10.001>`_
@ -270,6 +292,10 @@ Parallelism
Depletion
---------
- Anas Gul, K. S. Chaudri, R. Khan, and M. Azeen, "Development and verification
of LOOP: A Linkage of ORIGEN2.2 and OpenMC," *Ann. Nucl. Energy*, **99**,
321--327 (2017). `<http://dx.doi.org/10.1016/j.anucene.2016.09.016>`_
- Kai Huang, Hongchun Wu, Yunzhao Li, and Liangzhi Cao, "Generalized depletion
chain simplification based of significance analysis," *Proc. PHYSOR*, Sun
Valley, Idaho, May 1-5, 2016.

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

File diff suppressed because one or more lines are too long

3
docs/source/pythonapi/index.rst
View file

@ -373,6 +373,8 @@ Core Classes
openmc.data.CoherentElastic
openmc.data.FissionEnergyRelease
openmc.data.DataLibrary
openmc.data.Decay
openmc.data.FissionProductYields
openmc.data.WindowedMultipole
Core Functions
@ -477,6 +479,7 @@ Functions
openmc.data.endf.float_endf
openmc.data.endf.get_cont_record
openmc.data.endf.get_evaluations
openmc.data.endf.get_head_record
openmc.data.endf.get_tab1_record
openmc.data.endf.get_tab2_record

10
examples/python/pincell_multigroup/build-xml.py
View file

@ -1,3 +1,5 @@
import numpy as np
import openmc
import openmc.mgxs
@ -26,7 +28,7 @@ uo2_xsdata.set_total(
0.5644058])
uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02,
3.0020E-02, 1.1126E-01, 2.8278E-01])
uo2_xsdata.set_scatter_matrix(
scatter_matrix = np.array(
[[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
@ -34,6 +36,8 @@ uo2_xsdata.set_scatter_matrix(
[0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]])
scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
uo2_xsdata.set_scatter_matrix(scatter_matrix)
uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03,
1.85648E-02, 1.78084E-02, 8.30348E-02,
2.16004E-01])
@ -50,7 +54,7 @@ h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04,
1.9406E-03, 5.7416E-03, 1.5001E-02,
3.7239E-02])
h2o_xsdata.set_scatter_matrix(
scatter_matrix = np.array(
[[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
[0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
[0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
@ -58,6 +62,8 @@ h2o_xsdata.set_scatter_matrix(
[0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
h2o_xsdata.set_scatter_matrix(scatter_matrix)
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])

BIN
examples/xml/pincell_multigroup/mgxs.h5
View file

Binary file not shown.

8
openmc/clean_xml.py
View file

@ -65,25 +65,25 @@ def sort_xml_elements(tree):
tree.extend(sorted_elements)
def clean_xml_indentation(element, level=0):
def clean_xml_indentation(element, level=0, spaces_per_level=4):
"""
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*" "
i = "\n" + level*spaces_per_level*" "
if len(element):
if not element.text or not element.text.strip():
element.text = i + " "
element.text = i + spaces_per_level*" "
if not element.tail or not element.tail.strip():
element.tail = i
for sub_element in element:
clean_xml_indentation(sub_element, level+1)
clean_xml_indentation(sub_element, level+1, spaces_per_level)
if not sub_element.tail or not sub_element.tail.strip():
sub_element.tail = i

1
openmc/data/__init__.py
View file

@ -9,6 +9,7 @@ WMP_VERSION = b'v0.2'
from .data import *
from .neutron import *
from .decay import *
from .reaction import *
from .ace import *
from .angle_distribution import *

485
openmc/data/decay.py Normal file
View file

@ -0,0 +1,485 @@
from collections import Iterable, namedtuple
from io import StringIO
from math import log
from numbers import Real
import re
from warnings import warn
from six import string_types
import numpy as np
try:
from uncertainties import ufloat, unumpy, UFloat
except ImportError:
ufloat = UFloat = namedtuple('UFloat', ['nominal_value', 'std_dev'])
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
# Gives name and (change in A, change in Z) resulting from decay
_DECAY_MODES = {
0: ('gamma', (0, 0)),
1: ('beta-', (0, 1)),
2: ('ec/beta+', (0, -1)),
3: ('IT', (0, 0)),
4: ('alpha', (-4, -2)),
5: ('n', (-1, 0)),
6: ('sf', None),
7: ('p', (-1, -1)),
8: ('e-', (0, 0)),
9: ('xray', (0, 0)),
10: ('unknown', None)
}
_RADIATION_TYPES = {
0: 'gamma',
1: 'beta-',
2: 'ec/beta+',
4: 'alpha',
5: 'n',
6: 'sf',
7: 'p',
8: 'e-',
9: 'xray',
10: 'anti-neutrino',
11: 'neutrino'
}
def get_decay_modes(value):
"""Return sequence of decay modes given an ENDF RTYP value.
Parameters
----------
value : float
ENDF definition of sequence of decay modes
Returns
-------
list of str
List of successive decays, e.g. ('beta-', 'neutron')
"""
return [_DECAY_MODES[int(x)][0] for x in
str(value).strip('0').replace('.', '')]
class FissionProductYields(EqualityMixin):
"""Independent and cumulative fission product yields.
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Attributes
----------
cumulative : list of dict
Cumulative yields for each tabulated energy. Each item in the list is a
dictionary whose keys are nuclide names and values are cumulative
yields. The i-th dictionary corresponds to the i-th incident neutron
energy.
energies : Iterable of float or None
Energies at which fission product yields are tabulated.
independent : list of dict
Independent yields for each tabulated energy. Each item in the list is a
dictionary whose keys are nuclide names and values are independent
yields. The i-th dictionary corresponds to the i-th incident neutron
energy.
nuclide : dict
Properties of the fissioning nuclide.
Notes
-----
Neutron fission yields are typically not measured with a monoenergetic
source of neutrons. As such, if the fission yields are given at, e.g.,
0.0253 eV, one should interpret this as meaning that they are derived from a
typical thermal reactor flux spectrum as opposed to a monoenergetic source
at 0.0253 eV.
"""
def __init__(self, ev_or_filename):
# Define function that can be used to read both independent and
# cumulative yields
def get_yields(file_obj):
# Determine number of energies
n_energy = get_head_record(file_obj)[2]
energies = np.zeros(n_energy)
data = []
for i in range(n_energy):
# Determine i-th energy and number of products
items, values = get_list_record(file_obj)
energies[i] = items[0]
n_products = items[5]
# Get yields for i-th energy
yields = {}
for j in range(n_products):
Z, A = divmod(int(values[4*j]), 1000)
isomeric_state = int(values[4*j + 1])
name = ATOMIC_SYMBOL[Z] + str(A)
if isomeric_state > 0:
name += '_m{}'.format(isomeric_state)
yield_j = ufloat(values[4*j + 2], values[4*j + 3])
yields[name] = yield_j
data.append(yields)
return energies, data
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
# Assign basic nuclide properties
self.nuclide = {
'name': ev.gnd_name,
'atomic_number': ev.target['atomic_number'],
'mass_number': ev.target['mass_number'],
'isomeric_state': ev.target['isomeric_state']
}
# Read independent yields
if (8, 454) in ev.section:
file_obj = StringIO(ev.section[8, 454])
self.energies, self.independent = get_yields(file_obj)
# Read cumulative yields
if (8, 459) in ev.section:
file_obj = StringIO(ev.section[8, 459])
energies, self.cumulative = get_yields(file_obj)
assert np.all(energies == self.energies)
@classmethod
def from_endf(cls, ev_or_filename):
"""Generate fission product yield data from an ENDF evaluation
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Returns
-------
openmc.data.FissionProductYields
Fission product yield data
"""
return cls(ev_or_filename)
class DecayMode(EqualityMixin):
"""Radioactive decay mode.
Parameters
----------
parent : str
Parent decaying nuclide
modes : list of str
Successive decay modes
daughter_state : int
Metastable state of the daughter nuclide
energy : uncertainties.UFloat
Total decay energy in eV available in the decay process.
branching_ratio : uncertainties.UFloat
Fraction of the decay of the parent nuclide which proceeds by this mode.
Attributes
----------
branching_ratio : uncertainties.UFloat
Fraction of the decay of the parent nuclide which proceeds by this mode.
daughter : str
Name of daughter nuclide produced from decay
energy : uncertainties.UFloat
Total decay energy in eV available in the decay process.
modes : list of str
Successive decay modes
parent : str
Parent decaying nuclide
"""
def __init__(self, parent, modes, daughter_state, energy,
branching_ratio):
self._daughter_state = daughter_state
self.parent = parent
self.modes = modes
self.energy = energy
self.branching_ratio = branching_ratio
def __repr__(self):
return ('<DecayMode: ({}), {} -> {}, {}>'.format(
','.join(self.modes), self.parent, self.daughter,
self.branching_ratio))
@property
def branching_ratio(self):
return self._branching_ratio
@property
def daughter(self):
# Determine atomic number and mass number of parent
symbol, A = re.match(r'([A-Zn][a-z]*)(\d+)', self.parent).groups()
A = int(A)
Z = ATOMIC_NUMBER[symbol]
# Process changes
for mode in self.modes:
for name, changes in _DECAY_MODES.values():
if name == mode:
if changes is not None:
delta_A, delta_Z = changes
A += delta_A
Z += delta_Z
if self._daughter_state > 0:
return '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, self._daughter_state)
else:
return '{}{}'.format(ATOMIC_SYMBOL[Z], A)
@property
def energy(self):
return self._energy
@property
def modes(self):
return self._modes
@property
def parent(self):
return self._parent
@branching_ratio.setter
def branching_ratio(self, branching_ratio):
cv.check_type('branching ratio', branching_ratio, UFloat)
cv.check_greater_than('branching ratio',
branching_ratio.nominal_value, 0.0, True)
if branching_ratio.nominal_value == 0.0:
warn('Decay mode {} of parent {} has a zero branching ratio.'
.format(self.modes, self.parent))
cv.check_greater_than('branching ratio uncertainty',
branching_ratio.std_dev, 0.0, True)
self._branching_ratio = branching_ratio
@energy.setter
def energy(self, energy):
cv.check_type('decay energy', energy, UFloat)
cv.check_greater_than('decay energy', energy.nominal_value, 0.0, True)
cv.check_greater_than('decay energy uncertainty',
energy.std_dev, 0.0, True)
self._energy = energy
@modes.setter
def modes(self, modes):
cv.check_type('decay modes', modes, Iterable, string_types)
self._modes = modes
@parent.setter
def parent(self, parent):
cv.check_type('parent nuclide', parent, string_types)
self._parent = parent
class Decay(EqualityMixin):
"""Radioactive decay data.
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Attributes
----------
average_energies : dict
Average decay energies in eV of each type of radiation for decay heat
applications.
decay_constant : uncertainties.UFloat
Decay constant in inverse seconds.
half_life : uncertainties.UFloat
Half-life of the decay in seconds.
modes : list
Decay mode information for each mode of decay.
nuclide : dict
Dictionary describing decaying nuclide with keys 'name',
'excited_state', 'mass', 'stable', 'spin', and 'parity'.
spectra : dict
Resulting radiation spectra for each radiation type.
"""
def __init__(self, ev_or_filename):
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
file_obj = StringIO(ev.section[8, 457])
self.nuclide = {}
self.modes = []
self.spectra = {}
self.average_energies = {}
# Get head record
items = get_head_record(file_obj)
Z, A = divmod(items[0], 1000)
metastable = items[3]
self.nuclide['atomic_number'] = Z
self.nuclide['mass_number'] = A
self.nuclide['isomeric_state'] = metastable
if metastable > 0:
self.nuclide['name'] = '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A,
metastable)
else:
self.nuclide['name'] = '{}{}'.format(ATOMIC_SYMBOL[Z], A)
self.nuclide['mass'] = items[1] # AWR
self.nuclide['excited_state'] = items[2] # State of the original nuclide
self.nuclide['stable'] = (items[4] == 1) # Nucleus stability flag
# Determine if radioactive/stable
if not self.nuclide['stable']:
NSP = items[5] # Number of radiation types
# Half-life and decay energies
items, values = get_list_record(file_obj)
self.half_life = ufloat(items[0], items[1])
NC = items[4]//2
pairs = [x for x in zip(values[::2], values[1::2])]
ex = self.average_energies
ex['light'] = ufloat(*pairs[0])
ex['electromagnetic'] = ufloat(*pairs[1])
ex['heavy'] = ufloat(*pairs[2])
if NC == 17:
ex['beta-'] = ufloat(*pairs[3])
ex['beta+'] = ufloat(*pairs[4])
ex['auger'] = ufloat(*pairs[5])
ex['conversion'] = ufloat(*pairs[6])
ex['gamma'] = ufloat(*pairs[7])
ex['xray'] = ufloat(*pairs[8])
ex['Bremsstrahlung'] = ufloat(*pairs[9])
ex['annihilation'] = ufloat(*pairs[10])
ex['alpha'] = ufloat(*pairs[11])
ex['recoil'] = ufloat(*pairs[12])
ex['SF'] = ufloat(*pairs[13])
ex['neutron'] = ufloat(*pairs[14])
ex['proton'] = ufloat(*pairs[15])
ex['neutrino'] = ufloat(*pairs[16])
items, values = get_list_record(file_obj)
spin = items[0]
if spin == -77.777:
self.nuclide['spin'] = None
else:
self.nuclide['spin'] = spin
self.nuclide['parity'] = items[1] # Parity of the nuclide
# Decay mode information
n_modes = items[5] # Number of decay modes
for i in range(n_modes):
decay_type = get_decay_modes(values[6*i])
isomeric_state = int(values[6*i + 1])
energy = ufloat(*values[6*i + 2:6*i + 4])
branching_ratio = ufloat(*values[6*i + 4:6*(i + 1)])
mode = DecayMode(self.nuclide['name'], decay_type, isomeric_state,
energy, branching_ratio)
self.modes.append(mode)
discrete_type = {0.0: None, 1.0: 'allowed', 2.0: 'first-forbidden',
3.0: 'second-forbidden', 4.0: 'third-forbidden',
5.0: 'fourth-forbidden', 6.0: 'fifth-forbidden'}
# Read spectra
for i in range(NSP):
spectrum = {}
items, values = get_list_record(file_obj)
# Decay radiation type
spectrum['type'] = _RADIATION_TYPES[items[1]]
# Continuous spectrum flag
spectrum['continuous_flag'] = {0: 'discrete', 1: 'continuous',
2: 'both'}[items[2]]
spectrum['discrete_normalization'] = ufloat(*values[0:2])
spectrum['energy_average'] = ufloat(*values[2:4])
spectrum['continuous_normalization'] = ufloat(*values[4:6])
NER = items[5] # Number of tabulated discrete energies
if not spectrum['continuous_flag'] == 'continuous':
# Information about discrete spectrum
spectrum['discrete'] = []
for j in range(NER):
items, values = get_list_record(file_obj)
di = {}
di['energy'] = ufloat(*items[0:2])
di['from_mode'] = get_decay_modes(values[0])
di['type'] = discrete_type[values[1]]
di['intensity'] = ufloat(*values[2:4])
if spectrum['type'] == 'ec/beta+':
di['positron_intensity'] = ufloat(*values[4:6])
elif spectrum['type'] == 'gamma':
di['internal_pair'] = ufloat(*values[4:6])
if len(values) >= 8:
di['total_internal_conversion'] = ufloat(*values[6:8])
if len(values) == 12:
di['k_shell_conversion'] = ufloat(*values[8:10])
di['l_shell_conversion'] = ufloat(*values[10:12])
spectrum['discrete'].append(di)
if not spectrum['continuous_flag'] == 'discrete':
# Read continuous spectrum
ci = {}
params, ci['probability'] = get_tab1_record(file_obj)
ci['type'] = get_decay_modes(params[0])
# Read covariance (Ek, Fk) table
LCOV = params[3]
if LCOV != 0:
items, values = get_list_record(file_obj)
ci['covariance_lb'] = items[3]
ci['covariance'] = zip(values[0::2], values[1::2])
spectrum['continuous'] = ci
# Add spectrum to dictionary
self.spectra[spectrum['type']] = spectrum
else:
items, values = get_list_record(file_obj)
items, values = get_list_record(file_obj)
self.nuclide['spin'] = items[0]
self.nuclide['parity'] = items[1]
@property
def decay_constant(self):
if hasattr(self.half_life, 'n'):
return log(2.)/self.half_life
else:
mu, sigma = self.half_life
return ufloat(log(2.)/mu, log(2.)/mu**2*sigma)
@classmethod
def from_endf(cls, ev_or_filename):
"""Generate radioactive decay data from an ENDF evaluation
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Returns
-------
openmc.data.Decay
Radioactive decay data
"""
return cls(ev_or_filename)

60
openmc/data/endf.py
View file

@ -14,9 +14,11 @@ import os
from math import pi
from collections import OrderedDict, Iterable
from six import string_types
import numpy as np
from numpy.polynomial.polynomial import Polynomial
from .data import ATOMIC_SYMBOL
from .function import Tabulated1D, INTERPOLATION_SCHEME
from openmc.stats.univariate import Uniform, Tabular, Legendre
@ -47,16 +49,6 @@ SUM_RULES = {1: [2, 3],
_ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
def radiation_type(value):
p = {0: 'gamma', 1: 'beta-', 2: 'ec/beta+', 3: 'IT',
4: 'alpha', 5: 'neutron', 6: 'sf', 7: 'proton',
8: 'e-', 9: 'xray', 10: 'unknown'}
if value % 1.0 == 0:
return p[int(value)]
else:
return (p[int(value)], p[int(10*value % 10)])
def float_endf(s):
"""Convert string of floating point number in ENDF to float.
@ -258,14 +250,40 @@ def get_tab2_record(file_obj):
return params, Tabulated2D(breakpoints, interpolation)
def get_evaluations(filename):
"""Return a list of all evaluations within an ENDF file.
Parameters
----------
filename : str
Path to ENDF-6 formatted file
Returns
-------
list
A list of :class:`openmc.data.endf.Evaluation` instances.
"""
evaluations = []
with open(filename, 'r') as fh:
while True:
pos = fh.tell()
line = fh.readline()
if line[66:70] == ' -1':
break
fh.seek(pos)
evaluations.append(Evaluation(fh))
return evaluations
class Evaluation(object):
"""ENDF material evaluation with multiple files/sections
Parameters
----------
filename : str
Path to ENDF file to read
filename_or_obj : str or file-like
Path to ENDF file to read or an open file positioned at the start of an
ENDF material
Attributes
----------
@ -282,8 +300,11 @@ class Evaluation(object):
indicator (MOD).
"""
def __init__(self, filename):
fh = open(filename, 'r')
def __init__(self, filename_or_obj):
if isinstance(filename_or_obj, string_types):
fh = open(filename_or_obj, 'r')
else:
fh = filename_or_obj
self.section = {}
self.info = {}
self.target = {}
@ -313,6 +334,7 @@ class Evaluation(object):
# If end of material reached, exit loop
if MAT == 0:
fh.readline()
break
section_data = ''
@ -396,6 +418,16 @@ class Evaluation(object):
mod = 0
self.reaction_list.append((mf, mt, nc, mod))
@property
def gnd_name(self):
symbol = ATOMIC_SYMBOL[self.target['atomic_number']]
A = self.target['mass_number']
m = self.target['isomeric_state']
if m > 0:
return '{}{}_m{}'.format(symbol, A, m)
else:
return '{}{}'.format(symbol, A)
class Tabulated2D(object):
"""Metadata for a two-dimensional function.

3
openmc/data/laboratory.py
View file

@ -137,3 +137,6 @@ class LaboratoryAngleEnergy(AngleEnergy):
energy_out.append(energy_out_i)
return cls(tab2.breakpoints, tab2.interpolation, energy, mu, energy_out)
def to_hdf5(self, group):
raise NotImplementedError

17
openmc/data/reaction.py
View file

@ -55,13 +55,14 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)',
198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)',
649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
849: '(n,ac)'}
REACTION_NAME.update({i: '(n,n{})'.format(i-50) for i in range(50, 91)})
REACTION_NAME.update({i: '(n,p{})'.format(i-600) for i in range(600, 649)})
REACTION_NAME.update({i: '(n,d{})'.format(i-650) for i in range(650, 699)})
REACTION_NAME.update({i: '(n,t{})'.format(i-700) for i in range(700, 749)})
REACTION_NAME.update({i: '(n,3He{})'.format(i-750) for i in range(750, 799)})
REACTION_NAME.update({i: '(n,a{})'.format(i-800) for i in range(800, 849)})
849: '(n,ac)', 891: '(n,2nc)'}
REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)})
REACTION_NAME.update({i: '(n,p{})'.format(i - 600) for i in range(600, 649)})
REACTION_NAME.update({i: '(n,d{})'.format(i - 650) for i in range(650, 699)})
REACTION_NAME.update({i: '(n,t{})'.format(i - 700) for i in range(700, 749)})
REACTION_NAME.update({i: '(n,3He{})'.format(i - 750) for i in range(750, 799)})
REACTION_NAME.update({i: '(n,a{})'.format(i - 800) for i in range(800, 849)})
REACTION_NAME.update({i: '(n,2n{})'.format(i - 875) for i in range(875, 891)})
def _get_products(ev, mt):
@ -817,7 +818,7 @@ class Reaction(EqualityMixin):
Parameters
----------
group : h5py.Group
HDF5 group to write to
HDF5 group to read from
energy : dict
Dictionary whose keys are temperatures (e.g., '300K') and values are
arrays of energies at which cross sections are tabulated at.

336
openmc/filter.py
View file

@ -1,8 +1,7 @@
from abc import ABCMeta, abstractproperty
from abc import ABCMeta
from collections import Iterable, OrderedDict
import copy
from numbers import Real, Integral
import sys
from xml.etree import ElementTree as ET
from six import add_metaclass
@ -100,9 +99,13 @@ class Filter(object):
@classmethod
def _recursive_subclasses(cls):
"""Return all subclasses and their subclasses, etc."""
subs = cls.__subclasses__()
subsubs = [grand for s in subs for grand in s.__subclasses__()]
return subs + subsubs
all_subclasses = []
for subclass in cls.__subclasses__():
all_subclasses.append(subclass)
all_subclasses.extend(subclass._recursive_subclasses())
return all_subclasses
@classmethod
def from_hdf5(cls, group, **kwargs):
@ -774,7 +777,114 @@ class MeshFilter(Filter):
return df
class EnergyFilter(Filter):
class RealFilter(Filter):
"""Tally modifier that describes phase-space and other characteristics
Parameters
----------
bins : Iterable of Real
A grid of bin values.
Attributes
----------
bins : Iterable of Real
A grid of bin values.
num_bins : Integral
The number of filter bins
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
"""
def __gt__(self, other):
if type(self) is type(other):
# Compare largest/smallest bin edges in filters
# This logic is used when merging tallies with real filters
return self.bins[0] >= other.bins[-1]
else:
return super(RealFilter, self).__gt__(other)
@property
def num_bins(self):
return len(self.bins) - 1
@num_bins.setter
def num_bins(self, num_bins):
cv.check_type('filter num_bins', num_bins, Integral)
cv.check_greater_than('filter num_bins', num_bins, 0, equality=True)
self._num_bins = num_bins
def can_merge(self, other):
if type(self) is not type(other):
return False
if self.bins[0] == other.bins[-1]:
# This low edge coincides with other's high edge
return True
elif self.bins[-1] == other.bins[0]:
# This high edge coincides with other's low edge
return True
else:
return False
def merge(self, other):
if not self.can_merge(other):
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(type(self), type(other))
raise ValueError(msg)
# Merge unique filter bins
merged_bins = np.concatenate((self.bins, other.bins))
merged_bins = np.unique(merged_bins)
# Create a new filter with these bins
return type(self)(sorted(merged_bins))
def is_subset(self, other):
"""Determine if another filter is a subset of this filter.
If all of the bins in the other filter are included as bins in this
filter, then it is a subset of this filter.
Parameters
----------
other : openmc.Filter
The filter to query as a subset of this filter
Returns
-------
bool
Whether or not the other filter is a subset of this filter
"""
if type(self) is not type(other):
return False
elif len(self.bins) != len(other.bins):
return False
else:
return np.allclose(self.bins, other.bins)
def get_bin_index(self, filter_bin):
i = np.where(self.bins == filter_bin[1])[0]
if len(i) == 0:
msg = 'Unable to get the bin index for Filter since "{0}" ' \
'is not one of the bins'.format(filter_bin)
raise ValueError(msg)
else:
return i[0] - 1
def get_bin(self, bin_index):
cv.check_type('bin_index', bin_index, Integral)
cv.check_greater_than('bin_index', bin_index, 0, equality=True)
cv.check_less_than('bin_index', bin_index, self.num_bins)
# Construct 2-tuple of lower, upper bins for real-valued filters
return (self.bins[bin_index], self.bins[bin_index + 1])
class EnergyFilter(RealFilter):
"""Bins tally events based on incident particle energy.
Parameters
@ -794,23 +904,16 @@ class EnergyFilter(Filter):
"""
def __gt__(self, other):
if type(self) is type(other):
# Compare largest/smallest energy bin edges in energy filters
# This logic is used when merging tallies with energy filters
return self.bins[0] >= other.bins[-1]
def get_bin_index(self, filter_bin):
# Use lower energy bound to find index for RealFilters
deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
min_delta = np.min(deltas)
if min_delta < 1E-3:
return deltas.argmin() - 1
else:
return super(EnergyFilter, self).__gt__(other)
@property
def num_bins(self):
return len(self.bins) - 1
@num_bins.setter
def num_bins(self, num_bins):
cv.check_type('filter num_bins', num_bins, Integral)
cv.check_greater_than('filter num_bins', num_bins, 0, equality=True)
self._num_bins = num_bins
msg = 'Unable to get the bin index for Filter since "{0}" ' \
'is not one of the bins'.format(filter_bin)
raise ValueError(msg)
def check_bins(self, bins):
for edge in bins:
@ -829,62 +932,9 @@ class EnergyFilter(Filter):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, self.type)
'increasing'.format(bins, type(self))
raise ValueError(msg)
def can_merge(self, other):
if type(self) is not type(other):
return False
if self.bins[0] == other.bins[-1]:
# This low energy edge coincides with other's high energy edge
return True
elif self.bins[-1] == other.bins[0]:
# This high energy edge coincides with other's low energy edge
return True
else:
return False
def merge(self, other):
if not self.can_merge(other):
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(self.type, other.type)
raise ValueError(msg)
# Merge unique filter bins
merged_bins = np.concatenate((self.bins, other.bins))
merged_bins = np.unique(merged_bins)
# Create a new filter with these bins
return type(self)(sorted(merged_bins))
def is_subset(self, other):
if type(self) is not type(other):
return False
elif len(self.bins) != len(other.bins):
return False
else:
return np.allclose(self.bins, other.bins)
def get_bin_index(self, filter_bin):
# Use lower energy bound to find index for energy Filters
deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
min_delta = np.min(deltas)
if min_delta < 1E-3:
return deltas.argmin() - 1
else:
msg = 'Unable to get the bin index for Filter since "{0}" ' \
'is not one of the bins'.format(filter_bin)
raise ValueError(msg)
def get_bin(self, bin_index):
cv.check_type('bin_index', bin_index, Integral)
cv.check_greater_than('bin_index', bin_index, 0, equality=True)
cv.check_less_than('bin_index', bin_index, self.num_bins)
# Construct 2-tuple of lower, upper energies for energy(out) filters
return (self.bins[bin_index], self.bins[bin_index+1])
def get_pandas_dataframe(self, data_size, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
@ -1228,7 +1278,7 @@ class DistribcellFilter(Filter):
return df
class MuFilter(Filter):
class MuFilter(RealFilter):
"""Bins tally events based on particle scattering angle.
Parameters
@ -1256,8 +1306,82 @@ class MuFilter(Filter):
"""
def check_bins(self, bins):
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, type(self))
raise ValueError(msg)
elif edge < -1.:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is less than -1'.format(edge, type(self))
raise ValueError(msg)
elif edge > 1.:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is greater than 1'.format(edge, type(self))
raise ValueError(msg)
class PolarFilter(Filter):
# Check that bin edges are monotonically increasing
for index in range(1, len(bins)):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
This method constructs a Pandas DataFrame object for the filter with
columns annotated by filter bin information. This is a helper method
for :meth:`Tally.get_pandas_dataframe`.
Parameters
----------
data_size : Integral
The total number of bins in the tally corresponding to this filter
Returns
-------
pandas.DataFrame
A Pandas DataFrame with one column of the lower energy bound and one
column of upper energy bound for each filter bin. The number of
rows in the DataFrame is the same as the total number of bins in the
corresponding tally, with the filter bin appropriately tiled to map
to the corresponding tally bins.
Raises
------
ImportError
When Pandas is not installed
See also
--------
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Extract the lower and upper energy bounds, then repeat and tile
# them as necessary to account for other filters.
lo_bins = np.repeat(self.bins[:-1], self.stride)
hi_bins = np.repeat(self.bins[1:], self.stride)
tile_factor = data_size / len(lo_bins)
lo_bins = np.tile(lo_bins, tile_factor)
hi_bins = np.tile(hi_bins, tile_factor)
# Add the new energy columns to the DataFrame.
df.loc[:, self.short_name.lower() + ' low'] = lo_bins
df.loc[:, self.short_name.lower() + ' high'] = hi_bins
return df
class PolarFilter(RealFilter):
"""Bins tally events based on the incident particle's direction.
Parameters
@ -1285,6 +1409,30 @@ class PolarFilter(Filter):
"""
def check_bins(self, bins):
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, type(self))
raise ValueError(msg)
elif edge < 0.:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is less than 0'.format(edge, type(self))
raise ValueError(msg)
elif edge > np.pi:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is greater than pi'.format(edge, type(self))
raise ValueError(msg)
# Check that bin edges are monotonically increasing
for index in range(1, len(bins)):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
@ -1330,12 +1478,13 @@ class PolarFilter(Filter):
hi_bins = np.tile(hi_bins, tile_factor)
# Add the new angle columns to the DataFrame.
df.loc[:, self.type + ' low'] = lo_bins
df.loc[:, 'polar low'] = lo_bins
df.loc[:, 'polar high'] = hi_bins
return df
class AzimuthalFilter(Filter):
class AzimuthalFilter(RealFilter):
"""Bins tally events based on the incident particle's direction.
Parameters
@ -1363,6 +1512,30 @@ class AzimuthalFilter(Filter):
"""
def check_bins(self, bins):
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, type(self))
raise ValueError(msg)
elif edge < -np.pi:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is less than -pi'.format(edge, type(self))
raise ValueError(msg)
elif edge > np.pi:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is greater than pi'.format(edge, type(self))
raise ValueError(msg)
# Check that bin edges are monotonically increasing
for index in range(1, len(bins)):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, distribcell_paths=True):
"""Builds a Pandas DataFrame for the Filter's bins.
@ -1408,7 +1581,8 @@ class AzimuthalFilter(Filter):
hi_bins = np.tile(hi_bins, tile_factor)
# Add the new angle columns to the DataFrame.
df.loc[:, self.type + ' low'] = lo_bins
df.loc[:, 'azimuthal low'] = lo_bins
df.loc[:, 'azimuthal high'] = hi_bins
return df

4
openmc/macroscopic.py
View file

@ -1,5 +1,3 @@
import sys
from six import string_types
from openmc.checkvalue import check_type
@ -45,7 +43,7 @@ class Macroscopic(object):
return hash((self._name))
def __repr__(self):
string = 'Nuclide - {0}\n'.format(self._name)
string = 'Macroscopic - {0}\n'.format(self._name)
return string
@property

138
openmc/mgxs/library.py
View file

@ -59,15 +59,23 @@ class Library(object):
The spatial domain(s) for which MGXS in the Library are computed
correction : {'P0', None}
Apply the P0 correction to scattering matrices if set to 'P0'
scatter_format : {'legendre', 'histogram'}
Representation of the angular scattering distribution (default is
'legendre')
legendre_order : int
The highest legendre moment in the scattering matrices (default is 0)
The highest Legendre moment in the scattering matrix; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0)
histogram_bins : int
The number of equally-spaced bins for the histogram representation of
the angular scattering distribution; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16)
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure for energy condensation
num_delayed_groups : int
Number of delayed groups
estimator : str or None
The tally estimator used to compute multi-group cross sections. If None,
the default for each MGXS type is used.
The tally estimator used to compute multi-group cross sections.
If None, the default for each MGXS type is used.
tally_trigger : openmc.Trigger
An (optional) tally precision trigger given to each tally used to
compute the cross section
@ -101,7 +109,9 @@ class Library(object):
self._energy_groups = None
self._num_delayed_groups = 0
self._correction = 'P0'
self._scatter_format = 'legendre'
self._legendre_order = 0
self._histogram_bins = 16
self._tally_trigger = None
self._all_mgxs = OrderedDict()
self._sp_filename = None
@ -130,7 +140,9 @@ class Library(object):
clone._domain_type = self.domain_type
clone._domains = copy.deepcopy(self.domains)
clone._correction = self.correction
clone._scatter_format = self.scatter_format
clone._legendre_order = self.legendre_order
clone._histogram_bins = self.histogram_bins
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._num_delayed_groups = self.num_delayed_groups
clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
@ -209,10 +221,18 @@ class Library(object):
def correction(self):
return self._correction
@property
def scatter_format(self):
return self._scatter_format
@property
def legendre_order(self):
return self._legendre_order
@property
def histogram_bins(self):
return self._histogram_bins
@property
def tally_trigger(self):
return self._tally_trigger
@ -267,7 +287,7 @@ class Library(object):
def by_nuclide(self, by_nuclide):
cv.check_type('by_nuclide', by_nuclide, bool)
if by_nuclide == True and self.domain_type == 'mesh':
if by_nuclide and self.domain_type == 'mesh':
raise ValueError('Unable to create MGXS library by nuclide with '
'mesh domain')
@ -277,7 +297,7 @@ class Library(object):
def domain_type(self, domain_type):
cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES)
if self.by_nuclide == True and domain_type == 'mesh':
if self.by_nuclide and domain_type == 'mesh':
raise ValueError('Unable to create MGXS library by nuclide with '
'mesh domain')
@ -337,26 +357,72 @@ class Library(object):
def correction(self, correction):
cv.check_value('correction', correction, ('P0', None))
if correction == 'P0' and self.legendre_order > 0:
warn('The P0 correction will be ignored since the scattering '
'order "{}" is greater than zero'.format(self.legendre_order))
if self.scatter_format == 'legendre':
if correction == 'P0' and self.legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warn(msg)
elif self.scatter_format == 'histogram':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self._correction = correction
@scatter_format.setter
def scatter_format(self, scatter_format):
cv.check_value('scatter_format', scatter_format,
openmc.mgxs.MU_TREATMENTS)
if scatter_format == 'histogram' and self.correction == 'P0':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self.correction = None
self._scatter_format = scatter_format
@legendre_order.setter
def legendre_order(self, legendre_order):
cv.check_type('legendre_order', legendre_order, Integral)
cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
cv.check_greater_than('legendre_order', legendre_order, 0,
equality=True)
cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the scattering ' \
'order {} is greater than zero'.format(self.legendre_order)
warn(msg, RuntimeWarning)
self.correction = None
if self.scatter_format == 'legendre':
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warn(msg, RuntimeWarning)
self.correction = None
elif self.scatter_format == 'histogram':
msg = 'The legendre order will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self._legendre_order = legendre_order
@histogram_bins.setter
def histogram_bins(self, histogram_bins):
cv.check_type('histogram_bins', histogram_bins, Integral)
cv.check_greater_than('histogram_bins', histogram_bins, 0)
if self.scatter_format == 'legendre':
msg = 'The histogram bins will be ignored since the ' \
'scatter format is set to legendre'
warn(msg)
elif self.scatter_format == 'histogram':
if self.correction == 'P0':
msg = 'The P0 correction will be ignored since ' \
'a histogram representation of the scattering '\
'kernel is requested'
warn(msg, RuntimeWarning)
self.correction = None
self._histogram_bins = histogram_bins
@tally_trigger.setter
def tally_trigger(self, tally_trigger):
cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
@ -420,7 +486,7 @@ class Library(object):
mgxs.delayed_groups = None
else:
delayed_groups \
= list(range(1,self.num_delayed_groups+1))
= list(range(1, self.num_delayed_groups + 1))
mgxs.delayed_groups = delayed_groups
# If a tally trigger was specified, add it to the MGXS
@ -430,7 +496,9 @@ class Library(object):
# Specify whether to use a transport ('P0') correction
if isinstance(mgxs, openmc.mgxs.ScatterMatrixXS):
mgxs.correction = self.correction
mgxs.scatter_format = self.scatter_format
mgxs.legendre_order = self.legendre_order
mgxs.histogram_bins = self.histogram_bins
self.all_mgxs[domain.id][mgxs_type] = mgxs
@ -809,13 +877,13 @@ class Library(object):
return pickle.load(open(full_filename, 'rb'))
def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
order=None, subdomain=None):
subdomain=None):
"""Generates an openmc.XSdata object describing a multi-group cross section
dataset for writing to an openmc.MGXSLibrary object.
Note that this method does not build an XSdata
object with nested temperature tables. The temperature of each
XSdata object will be left at the default value of 300K.
XSdata object will be left at the default value of 294K.
Parameters
----------
@ -830,10 +898,6 @@ class Library(object):
Provide the macro or micro cross section in units of cm^-1 or
barns. Defaults to 'macro'. If the Library object is not tallied by
nuclide this will be set to 'macro' regardless.
order : int
Scattering order for this data entry. Default is None,
which will set the XSdata object to use the order of the
Library.
subdomain : iterable of int
This parameter is not used unless using a mesh domain. In that
case, the subdomain is an [i,j,k] index (1-based indexing) of the
@ -863,10 +927,6 @@ class Library(object):
cv.check_type('xsdata_name', xsdata_name, string_types)
cv.check_type('nuclide', nuclide, string_types)
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
cv.check_type('order', order, (type(None), Integral))
if order is not None:
cv.check_greater_than('order', order, 0, equality=True)
cv.check_less_than('order', order, 10, equality=True)
if subdomain is not None:
cv.check_iterable_type('subdomain', subdomain, Integral,
max_depth=3)
@ -888,16 +948,7 @@ class Library(object):
xsdata = openmc.XSdata(name, self.energy_groups)
xsdata.num_delayed_groups = self.num_delayed_groups
if order is None:
# Set the order to the Library's order (the defualt behavior)
xsdata.order = self.legendre_order
else:
# Set the order of the xsdata object to the minimum of
# the provided order or the Library's order.
xsdata.order = min(order, self.legendre_order)
# Right now only 'legendre' data and isotropic weighting is supported
self.scatter_format = 'legendre'
# Right now only isotropic weighting is supported
self.representation = 'isotropic'
if nuclide != 'total':
@ -1037,10 +1088,17 @@ class Library(object):
# accounted for approximately by using an adjusted
# absorption cross section.
if 'total' in self.mgxs_types or 'transport' in self.mgxs_types:
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] \
= np.subtract(xsdata._total[i], np.sum(
xsdata._scatter_matrix[i][0, :, :], axis=1))
if xsdata.scatter_format == 'legendre':
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] = \
np.subtract(xsdata._total[i], np.sum(
xsdata._scatter_matrix[i][0, :, :], axis=1))
elif xsdata.scatter_format == 'histogram':
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] = \
np.subtract(xsdata._total[i], np.sum(np.sum(
xsdata._scatter_matrix[i][:, :, :], axis=0),
axis=1))
return xsdata
@ -1321,7 +1379,7 @@ class Library(object):
'scattering matrix is not provided.')
# Total or transport can be present, but if using
# self.correction=="P0", then we should use transport.
if (((self.correction is "P0") and
if (((self.correction == "P0") and
('nu-transport' not in self.mgxs_types))):
error_flag = True
warn('A "nu-transport" MGXS type is required since a "P0" '

263
openmc/mgxs/mgxs.py
View file

@ -59,6 +59,12 @@ _DOMAINS = (openmc.Cell,
openmc.Material,
openmc.Mesh)
# Supported ScatterMatrixXS and NuScatterMatrixXS angular distribution types
MU_TREATMENTS = ('legendre', 'histogram')
# Maximum Legendre order supported by OpenMC
_MAX_LEGENDRE = 10
@add_metaclass(ABCMeta)
class MGXS(object):
@ -1526,9 +1532,19 @@ class MGXS(object):
else:
df = df.drop('score', axis=1)
# Determine if change-in-angle bins are included in the MGXS to
# properly tile the group boundaries
if 'mu low' in df:
# Find the length of the mu filters indirectly from the number
# of times the mu bins repeats.
num_mu = int(df.shape[0] /
df[df['mu low'] == df['mu low'][0]].shape[0])
else:
num_mu = 1
# Override energy groups bounds with indices
all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
all_groups = np.repeat(all_groups, len(query_nuclides))
all_groups = np.repeat(all_groups, len(query_nuclides) * num_mu)
if 'energy low [eV]' in df and 'energyout low [eV]' in df:
df.rename(columns={'energy low [eV]': 'group in'},
inplace=True)
@ -3207,8 +3223,8 @@ class NuScatterXS(MGXS):
class ScatterMatrixXS(MatrixMGXS):
r"""A scattering matrix multi-group cross section for one or more Legendre
moments.
r"""A scattering matrix multi-group cross section with the cosine of the
change-in-angle represented as one or more Legendre moments or a histogram.
This class can be used for both OpenMC input generation and tally data
post-processing to compute spatially-homogenized and energy-integrated
@ -3226,7 +3242,7 @@ class ScatterMatrixXS(MatrixMGXS):
For a spatial domain :math:`V`, incoming energy group
:math:`[E_{g'},E_{g'-1}]`, and outgoing energy group :math:`[E_g,E_{g-1}]`,
the scattering moments are calculated as:
the Legendre scattering moments are calculated as:
.. math::
@ -3266,9 +3282,18 @@ class ScatterMatrixXS(MatrixMGXS):
Attributes
----------
correction : 'P0' or None
Apply the P0 correction to scattering matrices if set to 'P0'
Apply the P0 correction to scattering matrices if set to 'P0'; this is
used only if :attr:`ScatterMatrixXS.scatter_format` is 'legendre'
scatter_format : {'legendre', or 'histogram'}
Representation of the angular scattering distribution (default is
'legendre')
legendre_order : int
The highest Legendre moment in the scattering matrix (default is 0)
The highest Legendre moment in the scattering matrix; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0)
histogram_bins : int
The number of equally-spaced bins for the histogram representation of
the angular scattering distribution; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16)
name : str, optional
Name of the multi-group cross section
rxn_type : str
@ -3335,7 +3360,9 @@ class ScatterMatrixXS(MatrixMGXS):
groups, by_nuclide, name)
self._rxn_type = 'scatter'
self._correction = 'P0'
self._scatter_format = 'legendre'
self._legendre_order = 0
self._histogram_bins = 16
self._hdf5_key = 'scatter matrix'
self._estimator = 'analog'
self._valid_estimators = ['analog']
@ -3343,26 +3370,39 @@ class ScatterMatrixXS(MatrixMGXS):
def __deepcopy__(self, memo):
clone = super(ScatterMatrixXS, self).__deepcopy__(memo)
clone._correction = self.correction
clone._scatter_format = self.scatter_format
clone._legendre_order = self.legendre_order
clone._histogram_bins = self.histogram_bins
return clone
@property
def correction(self):
return self._correction
@property
def scatter_format(self):
return self._scatter_format
@property
def legendre_order(self):
return self._legendre_order
@property
def histogram_bins(self):
return self._histogram_bins
@property
def scores(self):
scores = ['flux']
if self.correction == 'P0' and self.legendre_order == 0:
scores += ['{}-0'.format(self.rxn_type),
'{}-1'.format(self.rxn_type)]
else:
scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)]
if self.scatter_format == 'legendre':
if self.correction == 'P0' and self.legendre_order == 0:
scores += ['{}-0'.format(self.rxn_type),
'{}-1'.format(self.rxn_type)]
else:
scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)]
elif self.scatter_format == 'histogram':
scores += [self.rxn_type]
return scores
@ -3372,10 +3412,15 @@ class ScatterMatrixXS(MatrixMGXS):
energy = openmc.EnergyFilter(group_edges)
energyout = openmc.EnergyoutFilter(group_edges)
if self.correction == 'P0' and self.legendre_order == 0:
filters = [[energy], [energy, energyout], [energyout]]
else:
filters = [[energy], [energy, energyout]]
if self.scatter_format == 'legendre':
if self.correction == 'P0' and self.legendre_order == 0:
filters = [[energy], [energy, energyout], [energyout]]
else:
filters = [[energy], [energy, energyout]]
elif self.scatter_format == 'histogram':
bins = np.linspace(-1., 1., num=self.histogram_bins + 1,
endpoint=True)
filters = [[energy], [energy, energyout, openmc.MuFilter(bins)]]
return filters
@ -3383,20 +3428,24 @@ class ScatterMatrixXS(MatrixMGXS):
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
if self.scatter_format == 'legendre':
# If using P0 correction subtract scatter-1 from the diagonal
if self.correction == 'P0' and self.legendre_order == 0:
scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
energy_filter = copy.deepcopy(energy_filter)
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
self._rxn_rate_tally = scatter_p0 - scatter_p1
# If using P0 correction subtract scatter-1 from the diagonal
if self.correction == 'P0' and self.legendre_order == 0:
scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
energy_filter = copy.deepcopy(energy_filter)
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
self._rxn_rate_tally = scatter_p0 - scatter_p1
# Extract scattering moment reaction rate Tally
else:
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
self._rxn_rate_tally = self.tallies[tally_key]
# Extract scattering moment reaction rate Tally
else:
tally_key = '{}-P{}'.format(self.rxn_type,
self.legendre_order)
self._rxn_rate_tally = self.tallies[tally_key]
elif self.scatter_format == 'histogram':
# Extract scattering rate distribution tally
self._rxn_rate_tally = self.tallies[self.rxn_type]
self._rxn_rate_tally.sparse = self.sparse
@ -3406,27 +3455,53 @@ class ScatterMatrixXS(MatrixMGXS):
def correction(self, correction):
cv.check_value('correction', correction, ('P0', None))
if correction == 'P0' and self.legendre_order > 0:
msg = 'The P0 correction will be ignored since the scattering ' \
'order {} is greater than zero'.format(self.legendre_order)
if self.scatter_format == 'legendre':
if correction == 'P0' and self.legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warnings.warn(msg)
elif self.scatter_format == 'histogram':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warnings.warn(msg)
self._correction = correction
@scatter_format.setter
def scatter_format(self, scatter_format):
cv.check_value('scatter_format', scatter_format, MU_TREATMENTS)
self._scatter_format = scatter_format
@legendre_order.setter
def legendre_order(self, legendre_order):
cv.check_type('legendre_order', legendre_order, Integral)
cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
cv.check_greater_than('legendre_order', legendre_order, 0,
equality=True)
cv.check_less_than('legendre_order', legendre_order, _MAX_LEGENDRE,
equality=True)
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the scattering ' \
'order {} is greater than zero'.format(self.legendre_order)
warnings.warn(msg, RuntimeWarning)
self.correction = None
if self.scatter_format == 'legendre':
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warnings.warn(msg, RuntimeWarning)
self.correction = None
elif self.scatter_format == 'histogram':
msg = 'The legendre order will be ignored since the ' \
'scatter format is set to histogram'
warnings.warn(msg)
self._legendre_order = legendre_order
@histogram_bins.setter
def histogram_bins(self, histogram_bins):
cv.check_type('histogram_bins', histogram_bins, Integral)
cv.check_greater_than('histogram_bins', histogram_bins, 0)
self._histogram_bins = histogram_bins
def load_from_statepoint(self, statepoint):
"""Extracts tallies in an OpenMC StatePoint with the data needed to
compute multi-group cross sections.
@ -3456,12 +3531,16 @@ class ScatterMatrixXS(MatrixMGXS):
self._rxn_rate_tally = None
self._loaded_sp = False
# Expand scores to match the format in the statepoint
# e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
if self.correction != 'P0' or self.legendre_order != 0:
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
self.tallies[tally_key].scores = \
[self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)]
if self.scatter_format == 'legendre':
# Expand scores to match the format in the statepoint
# e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
if self.correction != 'P0' or self.legendre_order != 0:
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
self.tallies[tally_key].scores = \
[self.rxn_type + '-{}'.format(i)
for i in range(self.legendre_order + 1)]
elif self.scatter_format == 'histogram':
self.tallies[self.rxn_type].scores = [self.rxn_type]
super(ScatterMatrixXS, self).load_from_statepoint(statepoint)
@ -3508,7 +3587,7 @@ class ScatterMatrixXS(MatrixMGXS):
slice_xs._xs_tally = None
# Slice the Legendre order if needed
if legendre_order != 'same':
if legendre_order != 'same' and self.scatter_format == 'legendre':
cv.check_type('legendre_order', legendre_order, Integral)
cv.check_less_than('legendre_order', legendre_order,
self.legendre_order, equality=True)
@ -3517,7 +3596,8 @@ class ScatterMatrixXS(MatrixMGXS):
# Slice the scattering tally
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
expand_scores = \
[self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)]
[self.rxn_type + '-{}'.format(i)
for i in range(self.legendre_order + 1)]
slice_xs.tallies[tally_key] = \
slice_xs.tallies[tally_key].get_slice(scores=expand_scores)
@ -3548,7 +3628,8 @@ class ScatterMatrixXS(MatrixMGXS):
This method constructs a 5D NumPy array for the requested
multi-group cross section data for one or more subdomains
(1st dimension), energy groups in (2nd dimension), energy groups out
(3rd dimension), nuclides (4th dimension), and moments (5th dimension).
(3rd dimension), nuclides (4th dimension), and moments/histograms
(5th dimension).
NOTE: The scattering moments are not multiplied by the :math:`(2l+1)/2`
prefactor in the expansion of the scattering source into Legendre
@ -3637,7 +3718,7 @@ class ScatterMatrixXS(MatrixMGXS):
filter_bins.append((self.energy_groups.get_group_bounds(group),))
# Construct CrossScore for requested scattering moment
if moment != 'all':
if moment != 'all' and self.scatter_format == 'legendre':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
@ -3687,9 +3768,19 @@ class ScatterMatrixXS(MatrixMGXS):
else:
num_out_groups = len(out_groups)
if self.scatter_format == 'histogram':
num_mu_bins = self.histogram_bins
else:
num_mu_bins = 1
# Reshape tally data array with separate axes for domain and energy
num_subdomains = int(xs.shape[0] / (num_in_groups * num_out_groups))
new_shape = (num_subdomains, num_in_groups, num_out_groups)
num_subdomains = int(xs.shape[0] /
(num_mu_bins * num_in_groups * num_out_groups))
if self.scatter_format == 'histogram':
new_shape = (num_subdomains, num_in_groups, num_out_groups,
num_mu_bins)
else:
new_shape = (num_subdomains, num_in_groups, num_out_groups)
new_shape += xs.shape[1:]
xs = np.reshape(xs, new_shape)
@ -3700,11 +3791,22 @@ class ScatterMatrixXS(MatrixMGXS):
# Reverse data if user requested increasing energy groups since
# tally data is stored in order of increasing energies
if order_groups == 'increasing':
xs = xs[:, ::-1, ::-1, :]
xs = xs[:, ::-1, ::-1, ...]
if squeeze:
xs = np.squeeze(xs)
xs = np.atleast_2d(xs)
# We want to squeeze out everything but the in_groups, out_groups,
# and, if needed, num_mu_bins dimension. These must not be squeezed
# so 1-group problems have the correct shape.
if self.scatter_format == 'histogram':
axes = (5, 4, 0)
else:
axes = (4, 3, 0)
# Squeeze will return a ValueError if the axis has a size greater
# than 1, so try each axis in axes one at a time, catching the
# ValueError as needed.
for axis in axes:
if xs.shape[axis] == 1:
xs = np.squeeze(xs, axis=axis)
return xs
@ -3755,28 +3857,41 @@ class ScatterMatrixXS(MatrixMGXS):
df = super(ScatterMatrixXS, self).get_pandas_dataframe(
groups, nuclides, xs_type, distribcell_paths)
# Add a moment column to dataframe
if self.legendre_order > 0:
# Insert a column corresponding to the Legendre moments
moments = ['P{}'.format(i) for i in range(self.legendre_order+1)]
moments = np.tile(moments, int(df.shape[0] / len(moments)))
df['moment'] = moments
if self.scatter_format == 'legendre':
# Add a moment column to dataframe
if self.legendre_order > 0:
# Insert a column corresponding to the Legendre moments
moments = ['P{}'.format(i)
for i in range(self.legendre_order + 1)]
moments = np.tile(moments, int(df.shape[0] / len(moments)))
df['moment'] = moments
# Place the moment column before the mean column
columns = df.columns.tolist()
mean_index = [i for i, s in enumerate(columns) if 'mean' in s][0]
if self.domain_type == 'mesh':
df = df[columns[:mean_index] + [('moment', '')] + columns[mean_index:-1]]
else:
df = df[columns[:mean_index] + ['moment'] + columns[mean_index:-1]]
# Place the moment column before the mean column
columns = df.columns.tolist()
mean_index \
= [i for i, s in enumerate(columns) if 'mean' in s][0]
if self.domain_type == 'mesh':
df = df[columns[:mean_index] + [('moment', '')] +
columns[mean_index:-1]]
else:
df = df[columns[:mean_index] + ['moment'] +
columns[mean_index:-1]]
# Select rows corresponding to requested scattering moment
if moment != 'all':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
df = df[df['moment'] == 'P{}'.format(moment)]
# Select rows corresponding to requested scattering moment
if moment != 'all':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
df = df[df['moment'] == 'P{}'.format(moment)]
elif self.scatter_format == 'histogram':
# Replace the mu low and mu high columns with a single mu bin
del df['mu high']
df.rename(columns={'mu low': 'mu bins'}, inplace=True)
bins = [i + 1 for i in range(self.histogram_bins)]
bins = np.tile(bins, int(df.shape[0] / len(bins)))
df['mu bins'] = bins
return df
@ -3827,7 +3942,7 @@ class ScatterMatrixXS(MatrixMGXS):
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
if self.correction != 'P0':
if self.correction != 'P0' and self.scatter_format == 'legendre':
rxn_type = '{0} (P{1})'.format(self.rxn_type, moment)
else:
rxn_type = self.rxn_type

464
openmc/mgxs_library.py
View file

@ -1,6 +1,6 @@
from collections import Iterable
from numbers import Real, Integral
import sys
import os
from six import string_types
import numpy as np
@ -15,7 +15,7 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \
# Supported incoming particle MGXS angular treatment representations
_REPRESENTATIONS = ['isotropic', 'angle']
_SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
"[G'][DG]", "[G][G'][DG]"]
@ -42,7 +42,7 @@ class XSdata(object):
----------
name : str
Unique identifier for the xsdata object
aromic_weight_ratio : float
atomic_weight_ratio : float
Atomic weight ratio of an isotope. That is, the ratio of the mass
of the isotope to the mass of a single neutron.
temperatures : numpy.ndarray
@ -71,50 +71,50 @@ class XSdata(object):
Number of equal width angular bins that the polar angular domain is
subdivided into. This only applies when :attr:`XSdata.representation`
is "angle".
total : dict of numpy.ndarray
total : list of numpy.ndarray
Group-wise total cross section.
absorption : dict of numpy.ndarray
absorption : list of numpy.ndarray
Group-wise absorption cross section.
scatter_matrix : dict of numpy.ndarray
scatter_matrix : list of numpy.ndarray
Scattering moment matrices presented with the columns representing
incoming group and rows representing the outgoing group. That is,
down-scatter will be above the diagonal of the resultant matrix.
multiplicity_matrix : dict of numpy.ndarray
multiplicity_matrix : list of numpy.ndarray
Ratio of neutrons produced in scattering collisions to the neutrons
which undergo scattering collisions; that is, the multiplicity provides
the code with a scaling factor to account for neutrons produced in
(n,xn) reactions.
fission : dict of numpy.ndarray
fission : list of numpy.ndarray
Group-wise fission cross section.
kappa_fission : dict of numpy.ndarray
kappa_fission : list of numpy.ndarray
Group-wise kappa_fission cross section.
chi : dict of numpy.ndarray
chi : list of numpy.ndarray
Group-wise fission spectra ordered by increasing group index (i.e.,
fast to thermal). This attribute should be used if making the common
approximation that the fission spectra does not depend on incoming
energy. If the user does not wish to make this approximation, then
this should not be provided and this information included in the
:attr:`XSdata.nu_fission` attribute instead.
chi_prompt : dict of numpy.ndarray
chi_prompt : list of numpy.ndarray
Group-wise prompt fission spectra ordered by increasing group index
(i.e., fast to thermal). This attribute should be used if chi from
prompt and delayed neutrons is being set separately.
chi_delayed : dict of numpy.ndarray
chi_delayed : list of numpy.ndarray
Group-wise delayed fission spectra ordered by increasing group index
(i.e., fast to thermal). This attribute should be used if chi from
prompt and delayed neutrons is being set separately.
nu_fission : dict of numpy.ndarray
nu_fission : list of numpy.ndarray
Group-wise fission production cross section vector (i.e., if ``chi`` is
provided), or is the group-wise fission production matrix.
prompt_nu_fission : dict of numpy.ndarray
prompt_nu_fission : list of numpy.ndarray
Group-wise prompt fission production cross section vector.
delayed_nu_fission : dict of numpy.ndarray
delayed_nu_fission : list of numpy.ndarray
Group-wise delayed fission production cross section vector.
beta : dict of numpy.ndarray
beta : list of numpy.ndarray
Delayed-group-wise delayed neutron fraction cross section vector.
decay_rate : dict of numpy.ndarray
decay_rate : list of numpy.ndarray
Delayed-group-wise decay rate vector.
inverse_velocity : dict of numpy.ndarray
inverse_velocity : list of numpy.ndarray
Inverse of velocity, in units of sec/cm.
xs_shapes : dict of iterable of int
Dictionary with keys of _XS_SHAPES and iterable of int values with the
@ -134,7 +134,7 @@ class XSdata(object):
Note that some cross sections can be input in more than one shape so they
are listed multiple times:
[Order][G][G']: scatter_matrix
[G][G'][Order]: scatter_matrix
[G]: total, absorption, fission, kappa_fission, nu_fission,
prompt_nu_fission, delayed_nu_fission, inverse_velocity
@ -304,13 +304,13 @@ class XSdata(object):
self.energy_groups.num_groups,
self.num_delayed_groups)
self._xs_shapes["[Order][G][G']"] \
= (self.num_orders, self.energy_groups.num_groups,
self.energy_groups.num_groups)
self._xs_shapes["[G][G'][Order]"] \
= (self.energy_groups.num_groups,
self.energy_groups.num_groups, self.num_orders)
# If representation is by angle prepend num polar and num azim
if self.representation == 'angle':
for key,shapes in self._xs_shapes.items():
for key, shapes in self._xs_shapes.items():
self._xs_shapes[key] \
= (self.num_polar, self.num_azimuthal) + shapes
@ -338,7 +338,7 @@ class XSdata(object):
def num_delayed_groups(self, num_delayed_groups):
# Check validity of num_delayed_groups
check_type('num_delayed_groups', num_delayed_groups, int)
check_type('num_delayed_groups', num_delayed_groups, Integral)
check_less_than('num_delayed_groups', num_delayed_groups,
openmc.mgxs.MAX_DELAYED_GROUPS, equality=True)
check_greater_than('num_delayed_groups', num_delayed_groups, 0,
@ -724,7 +724,7 @@ class XSdata(object):
"""
# Get the accepted shapes for this xs
shapes = [self.xs_shapes["[Order][G][G']"]]
shapes = [self.xs_shapes["[G][G'][Order]"]]
# Convert to a numpy array so we can easily get the shape for checking
scatter = np.asarray(scatter)
@ -1521,33 +1521,41 @@ class XSdata(object):
check_type('temperature', temperature, Real)
check_value('temperature', temperature, self.temperatures)
if self.scatter_format != 'legendre':
msg = 'Anisotropic scattering representations other than ' \
'Legendre expansions have not yet been implemented in ' \
'openmc.mgxs.'
raise ValueError(msg)
# Set the value of scatter_format based on the same value within
# scatter
self.scatter_format = scatter.scatter_format
# If the user has not defined XSdata.order, then we will set
# the order based on the data within scatter.
# Otherwise, we will check to see that XSdata.order to match
# Otherwise, we will check to see that XSdata.order matches
# the order of scatter
if self.order is None:
self.order = scatter.legendre_order
else:
check_value('legendre_order', scatter.legendre_order,
[self.order])
if self.scatter_format == 'legendre':
if self.order is None:
self.order = scatter.legendre_order
else:
check_value('legendre_order', scatter.legendre_order,
[self.order])
elif self.scatter_format == 'histogram':
if self.order is None:
self.order = scatter.histogram_bins
else:
check_value('histogram_bins', scatter.histogram_bins,
[self.order])
i = np.where(self.temperatures == temperature)[0][0]
if self.representation == 'isotropic':
# Get the scattering orders in the outermost dimension
self._scatter_matrix[i] = np.zeros((self.num_orders,
self.energy_groups.num_groups,
self.energy_groups.num_groups))
for moment in range(self.num_orders):
self._scatter_matrix[i][moment, :, :] = \
if self.scatter_format == 'legendre':
# Get the scattering orders in the outermost dimension
self._scatter_matrix[i] = \
np.zeros(self.xs_shapes["[G][G'][Order]"])
for moment in range(self.num_orders):
self._scatter_matrix[i][:, :, moment] = \
scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
moment=moment, subdomains=subdomain)
else:
self._scatter_matrix[i] = \
scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
moment=moment, subdomains=subdomain)
subdomains=subdomain)
elif self.representation == 'angle':
msg = 'Angular-Dependent MGXS have not yet been implemented'
raise ValueError(msg)
@ -1625,6 +1633,10 @@ class XSdata(object):
scatt = scatter.get_xs(nuclides=nuclide,
xs_type=xs_type, moment=0,
subdomains=subdomain)
if scatter.scatter_format == 'histogram':
scatt = np.sum(scatt, axis=0)
if nuscatter.scatter_format == 'histogram':
nuscatt = np.sum(nuscatt, axis=0)
self._multiplicity_matrix[i] = np.divide(nuscatt, scatt)
elif self.representation == 'angle':
msg = 'Angular-Dependent MGXS have not yet been implemented'
@ -1641,6 +1653,7 @@ class XSdata(object):
HDF5 File (a root Group) to write to
"""
grp = file.create_group(self.name)
if self.atomic_weight_ratio is not None:
grp.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
@ -1656,7 +1669,7 @@ class XSdata(object):
if self.num_polar is not None:
grp.attrs['num_polar'] = self.num_polar
grp.attrs['scatter_shape'] = np.string_("[Order][G][G']")
grp.attrs['scatter_shape'] = np.string_("[G][G'][Order]")
if self.scatter_format is not None:
grp.attrs['scatter_format'] = np.string_(self.scatter_format)
if self.order is not None:
@ -1708,11 +1721,12 @@ class XSdata(object):
(self._delayed_nu_fission[i] is None or \
self._prompt_nu_fission[i] is None):
raise ValueError('nu-fission or prompt-nu-fission and '
'delayed-nu-fission data must be provided '
'when writing the HDF5 library')
'delayed-nu-fission data must be '
'provided when writing the HDF5 library')
if self._nu_fission[i] is not None:
xs_grp.create_dataset("nu-fission", data=self._nu_fission[i])
xs_grp.create_dataset("nu-fission",
data=self._nu_fission[i])
if self._prompt_nu_fission[i] is not None:
xs_grp.create_dataset("prompt-nu-fission",
@ -1726,7 +1740,8 @@ class XSdata(object):
xs_grp.create_dataset("beta", data=self._beta[i])
if self._decay_rate[i] is not None:
xs_grp.create_dataset("decay rate", data=self._decay_rate[i])
xs_grp.create_dataset("decay rate",
data=self._decay_rate[i])
if self._scatter_matrix[i] is None:
raise ValueError('Scatter matrix must be provided when '
@ -1735,97 +1750,88 @@ class XSdata(object):
# Get the sparse scattering data to print to the library
G = self.energy_groups.num_groups
if self.representation == 'isotropic':
g_out_bounds = np.zeros((G, 2), dtype=np.int)
for g_in in range(G):
nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
g_out_bounds[g_in, 0] = nz[0][0]
g_out_bounds[g_in, 1] = nz[0][-1]
# Now create the flattened scatter matrix array
matrix = self._scatter_matrix[i]
flat_scatt = []
for g_in in range(G):
for g_out in range(g_out_bounds[g_in, 0],
g_out_bounds[g_in, 1] + 1):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
# And write it.
scatt_grp = xs_grp.create_group('scatter_data')
scatt_grp.create_dataset("scatter_matrix",
data=np.array(flat_scatt))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
matrix = self._multiplicity_matrix[i][:, :]
flat_mult = []
for g_in in range(G):
for g_out in range(g_out_bounds[g_in, 0],
g_out_bounds[g_in, 1] + 1):
flat_mult.append(matrix[g_in, g_out])
scatt_grp.create_dataset("multiplicity matrix",
data=np.array(flat_mult))
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :] += 1
scatt_grp.create_dataset("g_min", data=g_out_bounds[:, 0])
scatt_grp.create_dataset("g_max", data=g_out_bounds[:, 1])
Np = 1
Na = 1
elif self.representation == 'angle':
Np = self.num_polar
Na = self.num_azimuthal
g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
for p in range(Np):
for a in range(Na):
for g_in in range(G):
matrix = self._scatter_matrix[i][p, a, 0, g_in, :]
nz = np.nonzero(matrix)
g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
for p in range(Np):
for a in range(Na):
for g_in in range(G):
if self.scatter_format == 'legendre':
if self.representation == 'isotropic':
matrix = \
self._scatter_matrix[i][g_in, :, 0]
elif self.representation == 'angle':
matrix = \
self._scatter_matrix[i][p, a, g_in, :, 0]
elif self.scatter_format == 'histogram':
if self.representation == 'isotropic':
matrix = \
np.sum(self._scatter_matrix[i][g_in, :, :],
axis=1)
elif self.representation == 'angle':
matrix = \
np.sum(self._scatter_matrix[i][p, a, g_in, :, :],
axis=1)
nz = np.nonzero(matrix)
# It is possible that there only zeros in matrix
# and therefore nz will be empty, in that case set
# g_out_bounds to 0s
if len(nz[0]) == 0:
g_out_bounds[p, a, g_in, :] = 0
else:
g_out_bounds[p, a, g_in, 0] = nz[0][0]
g_out_bounds[p, a, g_in, 1] = nz[0][-1]
# Now create the flattened scatter matrix array
flat_scatt = []
for p in range(Np):
for a in range(Na):
if self.representation == 'isotropic':
matrix = self._scatter_matrix[i][:, :, :]
elif self.representation == 'angle':
matrix = self._scatter_matrix[i][p, a, :, :, :]
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
for l in range(len(matrix[g_in, g_out, :])):
flat_scatt.append(matrix[g_in, g_out, l])
# And write it.
scatt_grp = xs_grp.create_group('scatter_data')
scatt_grp.create_dataset("scatter_matrix",
data=np.array(flat_scatt))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
flat_scatt = []
flat_mult = []
for p in range(Np):
for a in range(Na):
matrix = self._scatter_matrix[i][p, a, :, :, :]
if self.representation == 'isotropic':
matrix = self._multiplicity_matrix[i][:, :]
elif self.representation == 'angle':
matrix = self._multiplicity_matrix[i][p, a, :, :]
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
flat_mult.append(matrix[g_in, g_out])
# And write it.
scatt_grp = xs_grp.create_group('scatter_data')
scatt_grp.create_dataset("scatter_matrix",
data=np.array(flat_scatt))
scatt_grp.create_dataset("multiplicity_matrix",
data=np.array(flat_mult))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
flat_mult = []
for p in range(Np):
for a in range(Na):
matrix = self._multiplicity_matrix[i][p, a, :, :]
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
flat_mult.append(matrix[g_in, g_out])
# And write it.
scatt_grp.create_dataset("multiplicity_matrix",
data=np.array(flat_mult))
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :, :, :] += 1
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :, :, :] += 1
if self.representation == 'isotropic':
scatt_grp.create_dataset("g_min", data=g_out_bounds[0, 0, :, 0])
scatt_grp.create_dataset("g_max", data=g_out_bounds[0, 0, :, 1])
elif self.representation == 'angle':
scatt_grp.create_dataset("g_min", data=g_out_bounds[:, :, :, 0])
scatt_grp.create_dataset("g_max", data=g_out_bounds[:, :, :, 1])
@ -1834,6 +1840,143 @@ class XSdata(object):
xs_grp.create_dataset("inverse-velocity",
data=self._inverse_velocity[i])
@classmethod
def from_hdf5(cls, group, name, energy_groups, num_delayed_groups):
"""Generate XSdata object from an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to read from
name : str
Name of the mgxs data set.
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure
num_delayed_groups : int
Number of delayed groups
Returns
-------
openmc.XSdata
Multi-group cross section data
"""
# Get a list of all the subgroups which will contain our temperature
# strings
subgroups = group.keys()
temperatures = []
for subgroup in subgroups:
if subgroup != 'kTs':
temperatures.append(subgroup)
# To ensure the actual floating point temperature used when creating
# the new library is consistent with that used when originally creating
# the file, get the floating point temperatures straight from the kTs
# group.
kTs_group = group['kTs']
float_temperatures = []
for temperature in temperatures:
kT = kTs_group[temperature].value
float_temperatures.append(kT / openmc.data.K_BOLTZMANN)
attrs = group.attrs.keys()
if 'representation' in attrs:
representation = group.attrs['representation'].decode()
else:
representation = 'isotropic'
data = cls(name, energy_groups, float_temperatures, representation,
num_delayed_groups)
if 'scatter_format' in attrs:
data.scatter_format = group.attrs['scatter_format'].decode()
# Get the remaining optional attributes
if 'atomic_weight_ratio' in attrs:
data.atomic_weight_ratio = group.attrs['atomic_weight_ratio']
if 'order' in attrs:
data.order = group.attrs['order']
if data.representation == 'angle':
data.num_azimuthal = group.attrs['num_azimuthal']
data.num_polar = group.attrs['num_polar']
# Read the temperature-dependent datasets
for temp, float_temp in zip(temperatures, float_temperatures):
xs_types = ['total', 'absorption', 'fission', 'kappa-fission',
'chi', 'chi-prompt', 'chi-delayed', 'nu-fission',
'prompt-nu-fission', 'delayed-nu-fission', 'beta',
'decay rate', 'inverse-velocity']
temperature_group = group[temp]
for xs_type in xs_types:
set_func = 'set_' + xs_type.replace(' ', '_').replace('-', '_')
if xs_type in temperature_group:
getattr(data, set_func)(temperature_group[xs_type].value,
float_temp)
scatt_group = temperature_group['scatter_data']
# Get scatter matrix and 'un-flatten' it
g_max = scatt_group['g_max']
g_min = scatt_group['g_min']
flat_scatter = scatt_group['scatter_matrix'].value
scatter_matrix = np.zeros(data.xs_shapes["[G][G'][Order]"])
G = data.energy_groups.num_groups
if data.representation == 'isotropic':
Np = 1
Na = 1
elif data.representation == 'angle':
Np = data.num_polar
Na = data.num_azimuthal
flat_index = 0
for p in range(Np):
for a in range(Na):
for g_in in range(G):
if data.representation == 'isotropic':
g_mins = g_min[g_in]
g_maxs = g_max[g_in]
elif data.representation == 'angle':
g_mins = g_min[p, a, g_in]
g_maxs = g_max[p, a, g_in]
for g_out in range(g_mins - 1, g_maxs):
for ang in range(data.num_orders):
if data.representation == 'isotropic':
scatter_matrix[g_in, g_out, ang] = \
flat_scatter[flat_index]
elif data.representation == 'angle':
scatter_matrix[p, a, g_in, g_out, ang] = \
flat_scatter[flat_index]
flat_index += 1
data.set_scatter_matrix(scatter_matrix, float_temp)
# Repeat for multiplicity
if 'multiplicity_matrix' in scatt_group:
flat_mult = scatt_group['multiplicity_matrix'].value
mult_matrix = np.zeros(data.xs_shapes["[G][G']"])
flat_index = 0
for p in range(Np):
for a in range(Na):
for g_in in range(G):
if data.representation == 'isotropic':
g_mins = g_min[g_in]
g_maxs = g_max[g_in]
elif data.representation == 'angle':
g_mins = g_min[p, a, g_in]
g_maxs = g_max[p, a, g_in]
for g_out in range(g_mins - 1, g_maxs):
if data.representation == 'isotropic':
mult_matrix[g_in, g_out] = \
flat_mult[flat_index]
elif data.representation == 'angle':
mult_matrix[p, a, g_in, g_out] = \
flat_mult[flat_index]
flat_index += 1
data.set_multiplicity_matrix(mult_matrix, float_temp)
return data
class MGXSLibrary(object):
"""Multi-Group Cross Sections file used for an OpenMC simulation.
@ -1870,14 +2013,14 @@ class MGXSLibrary(object):
def num_delayed_groups(self):
return self._num_delayed_groups
@property
def temperatures(self):
return self._temperatures
@property
def xsdatas(self):
return self._xsdatas
@property
def names(self):
return [xsdata.name for xsdata in self.xsdatas]
@energy_groups.setter
def energy_groups(self, energy_groups):
check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups)
@ -1885,7 +2028,7 @@ class MGXSLibrary(object):
@num_delayed_groups.setter
def num_delayed_groups(self, num_delayed_groups):
check_type('num_delayed_groups', num_delayed_groups, int)
check_type('num_delayed_groups', num_delayed_groups, Integral)
check_greater_than('num_delayed_groups', num_delayed_groups, 0,
equality=True)
check_less_than('num_delayed_groups', num_delayed_groups,
@ -1914,7 +2057,7 @@ class MGXSLibrary(object):
self._xsdatas.append(xsdata)
def add_xsdatas(self, xsdatas):
"""Add multiple xsdatas to the file.
"""Add multiple XSdatas to the file.
Parameters
----------
@ -1945,6 +2088,27 @@ class MGXSLibrary(object):
self._xsdatas.remove(xsdata)
def get_by_name(self, name):
"""Access the XSdata objects by name
Parameters
----------
name : str
Name of openmc.XSdata object to obtain
Returns
-------
result : openmc.XSdata or None
Provides the matching XSdata object or None, if not found
"""
check_type("name", name, str)
result = None
for xsdata in self.xsdatas:
if name == xsdata.name:
result = xsdata
return result
def export_to_hdf5(self, filename='mgxs.h5'):
"""Create an hdf5 file that can be used for a simulation.
@ -1967,3 +2131,43 @@ class MGXSLibrary(object):
xsdata.to_hdf5(file)
file.close()
@classmethod
def from_hdf5(cls, filename=None):
"""Generate an MGXS Library from an HDF5 group or file
Parameters
----------
filename : str, optional
Name of HDF5 file containing MGXS data. Default is None.
If not provided, the value of the OPENMC_MG_CROSS_SECTIONS
environmental variable will be used
Returns
-------
openmc.MGXSLibrary
Multi-group cross section data object.
"""
# If filename is None, get the cross sections from the
# OPENMC_CROSS_SECTIONS environment variable
if filename is None:
filename = os.environ.get('OPENMC_MG_CROSS_SECTIONS')
# Check to make sure there was an environmental variable.
if filename is None:
raise ValueError("Either path or OPENMC_MG_CROSS_SECTIONS "
"environmental variable must be set")
check_type('filename', filename, str)
file = h5py.File(filename, 'r')
group_structure = file.attrs['group structure']
num_delayed_groups = file.attrs['delayed_groups']
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
data = cls(energy_groups, num_delayed_groups)
for group_name, group in file.items():
data.add_xsdata(openmc.XSdata.from_hdf5(group, group_name,
energy_groups,
num_delayed_groups))
return data

3
openmc/tallies.py
View file

@ -1320,7 +1320,8 @@ class Tally(object):
# Create list of 2-tuples for energy boundary bins
elif isinstance(self_filter, (openmc.EnergyFilter,
openmc.EnergyoutFilter)):
openmc.EnergyoutFilter, openmc.MuFilter,
openmc.PolarFilter, openmc.AzimuthalFilter)):
bins = []
for k in range(self_filter.num_bins):
bins.append((self_filter.bins[k], self_filter.bins[k+1]))

40
scripts/openmc-update-mgxs
View file

@ -23,12 +23,10 @@ optional arguments:
from __future__ import print_function
import os
from shutil import move
import warnings
import xml.etree.ElementTree as ET
import argparse
import h5py
import numpy as np
import openmc.mgxs_library
@ -51,7 +49,7 @@ def parse_args():
if args['output'] == '':
filename = args['input'].name
extension = filenameos.path.splitext()
extension = os.path.splitext(filename)
if extension == '.xml':
filename = filename[:filename.rfind('.')] + '.h5'
args['output'] = filename
@ -84,6 +82,8 @@ if __name__ == '__main__':
temp = tree.find('group_structure').text.strip()
temp = np.array(temp.split())
group_structure = temp.astype(np.float)
# Convert from MeV to eV
group_structure *= 1.e6
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
temp = tree.find('inverse-velocity')
if temp is not None:
@ -103,7 +103,7 @@ if __name__ == '__main__':
temperature = get_data(xsdata_elem, 'kT')
if temperature is not None:
temperature = \
float(temperature) / openmc.data.K_BOLTZMANN
float(temperature) / openmc.data.K_BOLTZMANN * 1.E6
else:
temperature = 294.
temperatures = [temperature]
@ -163,7 +163,7 @@ if __name__ == '__main__':
if temp is not None:
temp = np.array(temp.split(), dtype=float)
total = temp.astype(np.float)
total.shape = xsd[i].vector_shape
total.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_total(total, temperature)
if inverse_velocity is not None:
@ -172,41 +172,55 @@ if __name__ == '__main__':
temp = get_data(xsdata_elem, 'absorption')
temp = np.array(temp.split())
absorption = temp.astype(np.float)
absorption.shape = xsd[i].vector_shape
absorption.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_absorption(absorption, temperature)
temp = get_data(xsdata_elem, 'scatter')
temp = np.array(temp.split())
scatter = temp.astype(np.float)
scatter.shape = xsd[i].pn_matrix_shape
# This is now a flattened-array of something that started with a
# shape of [Order][G][G']; we need to unflatten and then switch the
# ordering
in_shape = (order_dim, energy_groups.num_groups,
energy_groups.num_groups)
if representation == 'angle':
in_shape = (n_pol, n_azi) + in_shape
scatter.shape = in_shape
scatter = np.swapaxes(scatter, 2, 3)
scatter = np.swapaxes(scatter, 3, 4)
else:
scatter.shape = in_shape
scatter = np.swapaxes(scatter, 0, 1)
scatter = np.swapaxes(scatter, 1, 2)
xsd[i].set_scatter_matrix(scatter, temperature)
temp = get_data(xsdata_elem, 'multiplicity')
if temp is not None:
temp = np.array(temp.split())
multiplicity = temp.astype(np.float)
multiplicity.shape = xsd[i].matrix_shape
multiplicity.shape = xsd[i].xs_shapes["[G][G']"]
xsd[i].set_multiplicity_matrix(multiplicity, temperature)
temp = get_data(xsdata_elem, 'fission')
if temp is not None:
temp = np.array(temp.split())
fission = temp.astype(np.float)
fission.shape = xsd[i].vector_shape
fission.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_fission(fission, temperature)
temp = get_data(xsdata_elem, 'kappa_fission')
if temp is not None:
temp = np.array(temp.split())
kappa_fission = temp.astype(np.float)
kappa_fission.shape = xsd[i].vector_shape
kappa_fission.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_kappa_fission(kappa_fission, temperature)
temp = get_data(xsdata_elem, 'chi')
if temp is not None:
temp = np.array(temp.split())
chi = temp.astype(np.float)
chi.shape = xsd[i].vector_shape
chi.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_chi(chi, temperature)
else:
chi = None
@ -216,9 +230,9 @@ if __name__ == '__main__':
temp = np.array(temp.split())
nu_fission = temp.astype(np.float)
if chi is not None:
nu_fission.shape = xsd[i].vector_shape
nu_fission.shape = xsd[i].xs_shapes['[G]']
else:
nu_fission.shape = xsd[i].matrix_shape
nu_fission.shape = xsd[i].xs_shapes["[G][G']"]
xsd[i].set_nu_fission(nu_fission, temperature)
# Build library as we go, but first we have enough to initialize it

1
setup.py
View file

@ -43,6 +43,7 @@ if have_setuptools:
# Optional dependencies
'extras_require': {
'decay': ['uncertainties'],
'pandas': ['pandas>=0.17.0'],
'sparse' : ['scipy'],
'vtk': ['vtk', 'silomesh'],

6
src/mgxs_header.F90
View file

@ -355,7 +355,7 @@ module mgxs_header
if (attribute_exists(xs_id, "scatter_shape")) then
call read_attribute(temp_str, xs_id, "scatter_shape")
temp_str = trim(temp_str)
if (to_lower(temp_str) /= "[order][g][g']") then
if (to_lower(temp_str) /= "[g][g'][order]") then
call fatal_error("Invalid scatter_shape option!")
end if
end if
@ -1995,8 +1995,8 @@ module mgxs_header
order_dim = order + 1
end if
! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing
! to ScattData
! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for
! passing to ScattData
allocate(input_scatt(energy_groups, this % n_azi, this % n_pol))
index = 1

4
src/scattdata_header.F90
View file

@ -274,7 +274,7 @@ contains
! Get this by summing the un-normalized P0 coefficient in matrix
! over all outgoing groups
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
this % scattxs(gin) = sum(matrix(gin) % data(:, :))
end do
allocate(energy(groups))
@ -282,7 +282,7 @@ contains
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = 1, groups
do gout = gmin(gin), gmax(gin)
norm = sum(matrix(gin) % data(:, gout))
energy(gin) % data(gout) = norm
if (norm /= ZERO) then

BIN
tests/1d_mgxs.h5
View file

Binary file not shown.

2
tests/test_mg_tallies/results_true.dat
View file

@ -1 +1 @@
3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989
864328b2c4f3c4bfa9c80c756acbedac6fbdf3d502c03c2f2a3e7461b774729c47a55341c9515eff246e1bd445485f0d0a05a5712663ee499046afdbe0c77aef

2
tests/test_tallies/inputs_true.dat
View file

@ -1 +1 @@
b1f616bc0e342e3886b4b27eca22c4d1e55e445979254d5ba5f9ea01048f53d0e9a01c575745ae130eb34108080d2256ab9e4cb926ea0050891a44ae37ecdd88
ac118c5796593df0f53491920875bbe77921882253e7f24942e365e5812d726fb3d3bab86da02b3aad93f6f8f20a090125ffa51f9fdf9aef84009702f9cc363d

2
tests/test_tallies/results_true.dat
View file

@ -1 +1 @@
ceb1420d48e097185dd0798a495676c8e968c37c6e75fd6232137a4df8f6fae33d8232b95c68f6479815a014d20ccfc7611f96b22fcb7da159b42346e00cdf4e
a8172f7492fcc69d21f53c20523ffd0ad2882a753a3e1d263d135153cc6b3ee8d8e38f8823a4a2eb3288d730579892875d2d19f85f083dd659f91414988115ea

52
tests/test_tallies/test_tallies.py
View file

@ -3,12 +3,14 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
from openmc.filter import *
from openmc import Mesh, Tally, Tallies
from openmc.source import Source
from openmc.stats import Box
class TalliesTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Build default materials/geometry
@ -21,33 +23,27 @@ class TalliesTestHarness(PyAPITestHarness):
self._input_set.settings.source = Source(space=Box(
[-160, -160, -183], [160, 160, 183]))
azimuthal_bins = (-3.1416, -1.8850, -0.6283, 0.6283, 1.8850, 3.1416)
azimuthal_filter1 = AzimuthalFilter(azimuthal_bins)
azimuthal_bins = (-3.14159, -1.8850, -0.6283, 0.6283, 1.8850, 3.14159)
azimuthal_filter = AzimuthalFilter(azimuthal_bins)
azimuthal_tally1 = Tally()
azimuthal_tally1.filters = [azimuthal_filter1]
azimuthal_tally1.filters = [azimuthal_filter]
azimuthal_tally1.scores = ['flux']
azimuthal_tally1.estimator = 'tracklength'
azimuthal_tally2 = Tally()
azimuthal_tally2.filters = [azimuthal_filter1]
azimuthal_tally2.filters = [azimuthal_filter]
azimuthal_tally2.scores = ['flux']
azimuthal_tally2.estimator = 'analog'
azimuthal_filter2 = AzimuthalFilter(5)
azimuthal_tally3 = Tally()
azimuthal_tally3.filters = [azimuthal_filter2]
azimuthal_tally3.scores = ['flux']
azimuthal_tally3.estimator = 'tracklength'
mesh_2x2 = Mesh(mesh_id=1)
mesh_2x2.lower_left = [-182.07, -182.07]
mesh_2x2.upper_right = [182.07, 182.07]
mesh_2x2.dimension = [2, 2]
mesh_filter = MeshFilter(mesh_2x2)
azimuthal_tally4 = Tally()
azimuthal_tally4.filters = [azimuthal_filter2, mesh_filter]
azimuthal_tally4.scores = ['flux']
azimuthal_tally4.estimator = 'tracklength'
azimuthal_tally3 = Tally()
azimuthal_tally3.filters = [azimuthal_filter, mesh_filter]
azimuthal_tally3.scores = ['flux']
azimuthal_tally3.estimator = 'tracklength'
cellborn_tally = Tally()
cellborn_tally.filters = [CellbornFilter((10, 21, 22, 23))]
@ -76,20 +72,17 @@ class TalliesTestHarness(PyAPITestHarness):
material_tally.filters = [MaterialFilter((1, 2, 3, 4))]
material_tally.scores = ['total']
mu_bins = (-1.0, -0.5, 0.0, 0.5, 1.0)
mu_filter = MuFilter(mu_bins)
mu_tally1 = Tally()
mu_tally1.filters = [MuFilter((-1.0, -0.5, 0.0, 0.5, 1.0))]
mu_tally1.filters = [mu_filter]
mu_tally1.scores = ['scatter', 'nu-scatter']
mu_filter = MuFilter(5)
mu_tally2 = Tally()
mu_tally2.filters = [mu_filter]
mu_tally2.filters = [mu_filter, mesh_filter]
mu_tally2.scores = ['scatter', 'nu-scatter']
mu_tally3 = Tally()
mu_tally3.filters = [mu_filter, mesh_filter]
mu_tally3.scores = ['scatter', 'nu-scatter']
polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.1416)
polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.14159)
polar_filter = PolarFilter(polar_bins)
polar_tally1 = Tally()
polar_tally1.filters = [polar_filter]
@ -101,17 +94,11 @@ class TalliesTestHarness(PyAPITestHarness):
polar_tally2.scores = ['flux']
polar_tally2.estimator = 'analog'
polar_filter2 = PolarFilter((5,))
polar_tally3 = Tally()
polar_tally3.filters = [polar_filter2]
polar_tally3.filters = [polar_filter, mesh_filter]
polar_tally3.scores = ['flux']
polar_tally3.estimator = 'tracklength'
polar_tally4 = Tally()
polar_tally4.filters = [polar_filter2, mesh_filter]
polar_tally4.scores = ['flux']
polar_tally4.estimator = 'tracklength'
universe_tally = Tally()
universe_tally.filters = [UniverseFilter((1, 2, 3, 4, 6, 8))]
universe_tally.scores = ['total']
@ -173,10 +160,9 @@ class TalliesTestHarness(PyAPITestHarness):
self._input_set.tallies = Tallies()
self._input_set.tallies += (
[azimuthal_tally1, azimuthal_tally2, azimuthal_tally3,
azimuthal_tally4, cellborn_tally, dg_tally, energy_tally,
energyout_tally, transfer_tally, material_tally, mu_tally1,
mu_tally2, mu_tally3, polar_tally1, polar_tally2, polar_tally3,
polar_tally4, universe_tally])
cellborn_tally, dg_tally, energy_tally, energyout_tally,
transfer_tally, material_tally, mu_tally1, mu_tally2,
polar_tally1, polar_tally2, polar_tally3, universe_tally])
self._input_set.tallies += score_tallies
self._input_set.tallies += flux_tallies
self._input_set.tallies += (scatter_tally1, scatter_tally2)