mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge remote-tracking branch 'upstream/develop' into linlinenergy
This commit is contained in:
Sterling Harper
commit
e8c9fe664a
42 changed files with 2456 additions and 663 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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>`_
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
File diff suppressed because one or more lines are too long
|
|
@ -374,6 +374,8 @@ Core Classes
|
|||
openmc.data.CoherentElastic
|
||||
openmc.data.FissionEnergyRelease
|
||||
openmc.data.DataLibrary
|
||||
openmc.data.Decay
|
||||
openmc.data.FissionProductYields
|
||||
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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])
|
||||
|
|
|
|||
Binary file not shown.
|
|
@ -25,6 +25,7 @@ from openmc.statepoint import *
|
|||
from openmc.summary import *
|
||||
from openmc.particle_restart import *
|
||||
from openmc.mixin import *
|
||||
from openmc.plotter import *
|
||||
|
||||
try:
|
||||
from openmc.opencg_compatible import *
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -368,7 +368,7 @@ class CMFD(object):
|
|||
@cmfd_mesh.setter
|
||||
def cmfd_mesh(self, mesh):
|
||||
check_type('CMFD mesh', mesh, CMFDMesh)
|
||||
self._mesh = mesh
|
||||
self._cmfd_mesh = mesh
|
||||
|
||||
@norm.setter
|
||||
def norm(self, norm):
|
||||
|
|
@ -446,8 +446,8 @@ class CMFD(object):
|
|||
element.text = str(self._ktol)
|
||||
|
||||
def _create_mesh_subelement(self):
|
||||
if self._mesh is not None:
|
||||
xml_element = self._mesh._get_xml_element()
|
||||
if self._cmfd_mesh is not None:
|
||||
xml_element = self._cmfd_mesh._get_xml_element()
|
||||
self._cmfd_file.append(xml_element)
|
||||
|
||||
def _create_norm_subelement(self):
|
||||
|
|
|
|||
|
|
@ -6,6 +6,7 @@ HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
|
|||
|
||||
from .data import *
|
||||
from .neutron import *
|
||||
from .decay import *
|
||||
from .reaction import *
|
||||
from .ace import *
|
||||
from .angle_distribution import *
|
||||
|
|
|
|||
|
|
@ -186,3 +186,9 @@ K_BOLTZMANN = 8.6173324e-5
|
|||
|
||||
# Used for converting units in ACE data
|
||||
EV_PER_MEV = 1.0e6
|
||||
|
||||
# Avogadro's constant from CODATA 2010
|
||||
AVOGADRO = 6.02214129E23
|
||||
|
||||
# Neutron mass from CODATA 2010 in units of amu
|
||||
NEUTRON_MASS = 1.008664916
|
||||
|
|
|
|||
485
openmc/data/decay.py
Normal file
485
openmc/data/decay.py
Normal 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)
|
||||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,10 +1,12 @@
|
|||
import os
|
||||
import xml.etree.ElementTree as ET
|
||||
from six import string_types
|
||||
|
||||
import h5py
|
||||
|
||||
from openmc.mixin import EqualityMixin
|
||||
from openmc.clean_xml import clean_xml_indentation
|
||||
from openmc.checkvalue import check_type
|
||||
|
||||
|
||||
class DataLibrary(EqualityMixin):
|
||||
|
|
@ -95,13 +97,14 @@ class DataLibrary(EqualityMixin):
|
|||
method='xml')
|
||||
|
||||
@classmethod
|
||||
def from_xml(cls, path):
|
||||
def from_xml(cls, path=None):
|
||||
"""Read cross section data library from an XML file.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
path : str
|
||||
Path to XML file to read.
|
||||
path : str, optional
|
||||
Path to XML file to read. If not provided, the
|
||||
`OPENMC_CROSS_SECTIONS` environment variable will be used.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -112,6 +115,18 @@ class DataLibrary(EqualityMixin):
|
|||
|
||||
data = cls()
|
||||
|
||||
# If path is None, get the cross sections from the
|
||||
# OPENMC_CROSS_SECTIONS environment variable
|
||||
if path is None:
|
||||
path = os.environ.get('OPENMC_CROSS_SECTIONS')
|
||||
|
||||
# Check to make sure there was an environmental variable.
|
||||
if path is None:
|
||||
raise ValueError("Either path or OPENMC_CROSS_SECTIONS "
|
||||
"environmental variable must be set")
|
||||
|
||||
check_type('path', path, string_types)
|
||||
|
||||
tree = ET.parse(path)
|
||||
root = tree.getroot()
|
||||
if root.find('directory') is not None:
|
||||
|
|
|
|||
|
|
@ -511,6 +511,9 @@ class IncidentNeutron(EqualityMixin):
|
|||
rxs = [data[mt] for mt in SUM_RULES[mt_sum] if mt in data]
|
||||
if len(rxs) > 0:
|
||||
data.summed_reactions[mt_sum] = rx = Reaction(mt_sum)
|
||||
if rx.mt == 18 and 'total_nu' in group:
|
||||
tgroup = group['total_nu']
|
||||
rx.derived_products.append(Product.from_hdf5(tgroup))
|
||||
for T in data.temperatures:
|
||||
rx.xs[T] = Sum([rx_i.xs[T] for rx_i in rxs])
|
||||
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -1,13 +1,12 @@
|
|||
from collections import OrderedDict
|
||||
import re
|
||||
import sys
|
||||
import os
|
||||
|
||||
from six import string_types
|
||||
from xml.etree import ElementTree as ET
|
||||
|
||||
import openmc
|
||||
from openmc.checkvalue import check_type, check_length
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data import NATURAL_ABUNDANCE, atomic_mass
|
||||
|
||||
|
||||
|
|
@ -80,8 +79,8 @@ class Element(object):
|
|||
|
||||
@name.setter
|
||||
def name(self, name):
|
||||
check_type('element name', name, string_types)
|
||||
check_length('element name', name, 1, 2)
|
||||
cv.check_type('element name', name, string_types)
|
||||
cv.check_length('element name', name, 1, 2)
|
||||
self._name = name
|
||||
|
||||
@scattering.setter
|
||||
|
|
@ -254,6 +253,6 @@ class Element(object):
|
|||
for nuclide, abundance in abundances.items():
|
||||
nuc = openmc.Nuclide(nuclide)
|
||||
nuc.scattering = self.scattering
|
||||
isotopes.append((nuc, percent*abundance, percent_type))
|
||||
isotopes.append((nuc, percent * abundance, percent_type))
|
||||
|
||||
return isotopes
|
||||
|
|
|
|||
336
openmc/filter.py
336
openmc/filter.py
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -3,9 +3,9 @@ from copy import deepcopy
|
|||
from numbers import Real, Integral
|
||||
import warnings
|
||||
from xml.etree import ElementTree as ET
|
||||
import sys
|
||||
|
||||
from six import string_types
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
import openmc.data
|
||||
|
|
@ -595,7 +595,7 @@ class Material(object):
|
|||
nuclides = OrderedDict()
|
||||
|
||||
for nuclide, density, density_type in self._nuclides:
|
||||
nuclides[nuclide.name] = (nuclide, density)
|
||||
nuclides[nuclide.name] = (nuclide, density, density_type)
|
||||
|
||||
for ele, ele_pct, ele_pct_type, enr in self._elements:
|
||||
|
||||
|
|
@ -606,6 +606,80 @@ class Material(object):
|
|||
|
||||
return nuclides
|
||||
|
||||
def get_nuclide_atom_densities(self):
|
||||
"""Returns all nuclides in the material and their atomic densities in
|
||||
units of atom/b-cm
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuclides : dict
|
||||
Dictionary whose keys are nuclide names and values are tuples of
|
||||
(nuclide, density in atom/b-cm)
|
||||
|
||||
"""
|
||||
|
||||
# Expand elements in to nuclides
|
||||
nuclides = self.get_nuclide_densities()
|
||||
|
||||
sum_density = False
|
||||
if self.density_units == 'sum':
|
||||
sum_density = True
|
||||
density = 0.
|
||||
elif self.density_units == 'macro':
|
||||
density = self.density
|
||||
elif self.density_units == 'g/cc' or self.density_units == 'g/cm3':
|
||||
density = -self.density
|
||||
elif self.density_units == 'kg/m3':
|
||||
density = -0.001 * self.density
|
||||
elif self.density_units == 'atom/b-cm':
|
||||
density = self.density
|
||||
elif self.density_units == 'atom/cm3' or self.density_units == 'atom/cc':
|
||||
density = 1.E-24 * self.density
|
||||
|
||||
# For ease of processing split out nuc, nuc_density,
|
||||
# and nuc_density_type in to separate arrays
|
||||
nucs = []
|
||||
nuc_densities = []
|
||||
nuc_density_types = []
|
||||
for nuclide in nuclides.items():
|
||||
nuc, nuc_density, nuc_density_type = nuclide[1]
|
||||
nucs.append(nuc)
|
||||
nuc_densities.append(nuc_density)
|
||||
nuc_density_types.append(nuc_density_type)
|
||||
|
||||
if sum_density:
|
||||
density = np.sum(nuc_densities)
|
||||
percent_in_atom = np.all(nuc_density_types == 'ao')
|
||||
density_in_atom = density > 0.
|
||||
sum_percent = 0.
|
||||
|
||||
awrs = []
|
||||
for n, nuclide in enumerate(nuclides.items()):
|
||||
awr = openmc.data.atomic_mass(nuclide[0])
|
||||
if awr is not None:
|
||||
awrs.append(awr / openmc.data.NEUTRON_MASS)
|
||||
else:
|
||||
raise ValueError(nuclide[0] + " is invalid")
|
||||
|
||||
# Now that we have the awr, lets finish calculating densities
|
||||
sum_percent = np.sum(nuc_densities)
|
||||
nuc_densities = nuc_densities / sum_percent
|
||||
if not density_in_atom:
|
||||
sum_percent = 0.
|
||||
for n, nuc in enumerate(nucs):
|
||||
x = nuc_densities[n]
|
||||
sum_percent += x * awrs[n]
|
||||
sum_percent = 1. / sum_percent
|
||||
density = -density * sum_percent * \
|
||||
openmc.data.AVOGADRO / openmc.data.NEUTRON_MASS * 1.E-24
|
||||
nuc_densities = density * nuc_densities
|
||||
|
||||
nuclides = OrderedDict()
|
||||
for n, nuc in enumerate(nucs):
|
||||
nuclides[nuc] = (nuc, nuc_densities[n])
|
||||
|
||||
return nuclides
|
||||
|
||||
def _get_nuclide_xml(self, nuclide, distrib=False):
|
||||
xml_element = ET.Element("nuclide")
|
||||
xml_element.set("name", nuclide[0].name)
|
||||
|
|
|
|||
|
|
@ -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" '
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
@ -192,7 +198,9 @@ class MGXS(object):
|
|||
clone._rxn_rate_tally = copy.deepcopy(self._rxn_rate_tally, memo)
|
||||
clone._xs_tally = copy.deepcopy(self._xs_tally, memo)
|
||||
clone._sparse = self.sparse
|
||||
clone._loaded_sp = self._loaded_sp
|
||||
clone._derived = self.derived
|
||||
clone._hdf5_key = self._hdf5_key
|
||||
|
||||
clone._tallies = OrderedDict()
|
||||
for tally_type, tally in self.tallies.items():
|
||||
|
|
@ -325,7 +333,7 @@ class MGXS(object):
|
|||
|
||||
@property
|
||||
def num_subdomains(self):
|
||||
if self.domain_type.startswith('avg('):
|
||||
if self.domain_type.startswith('sum('):
|
||||
domain_type = self.domain_type[4:-1]
|
||||
else:
|
||||
domain_type = self.domain_type
|
||||
|
|
@ -789,16 +797,22 @@ class MGXS(object):
|
|||
if not isinstance(subdomains, string_types):
|
||||
cv.check_iterable_type('subdomains', subdomains, Integral,
|
||||
max_depth=3)
|
||||
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
subdomain_bins = []
|
||||
for subdomain in subdomains:
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
filter_bins.append((subdomain,))
|
||||
subdomain_bins.append(subdomain)
|
||||
filter_bins.append(tuple(subdomain_bins))
|
||||
|
||||
# Construct list of energy group bounds tuples for all requested groups
|
||||
if not isinstance(groups, string_types):
|
||||
cv.check_iterable_type('groups', groups, Integral)
|
||||
filters.append(openmc.EnergyFilter)
|
||||
energy_bins = []
|
||||
for group in groups:
|
||||
filters.append(openmc.EnergyFilter)
|
||||
filter_bins.append((self.energy_groups.get_group_bounds(group),))
|
||||
energy_bins.append(
|
||||
(self.energy_groups.get_group_bounds(group),))
|
||||
filter_bins.append(tuple(energy_bins))
|
||||
|
||||
# Construct a collection of the nuclides to retrieve from the xs tally
|
||||
if self.by_nuclide:
|
||||
|
|
@ -958,30 +972,27 @@ class MGXS(object):
|
|||
# Construct a collection of the subdomain filter bins to average across
|
||||
if not isinstance(subdomains, string_types):
|
||||
cv.check_iterable_type('subdomains', subdomains, Integral)
|
||||
subdomains = [(subdomain,) for subdomain in subdomains]
|
||||
subdomains = [tuple(subdomains)]
|
||||
elif self.domain_type == 'distribcell':
|
||||
subdomains = np.arange(self.num_subdomains)
|
||||
subdomains = [i for i in range(self.num_subdomains)]
|
||||
subdomains = [tuple(subdomains)]
|
||||
else:
|
||||
subdomains = None
|
||||
|
||||
# Clone this MGXS to initialize the subdomain-averaged version
|
||||
avg_xs = copy.deepcopy(self)
|
||||
avg_xs._rxn_rate_tally = None
|
||||
avg_xs._xs_tally = None
|
||||
|
||||
if self.derived:
|
||||
avg_xs._rxn_rate_tally = avg_xs.rxn_rate_tally.average(
|
||||
filter_type=_DOMAIN_TO_FILTER[self.domain_type],
|
||||
filter_bins=subdomains)
|
||||
else:
|
||||
avg_xs._rxn_rate_tally = None
|
||||
avg_xs._xs_tally = None
|
||||
# Average each of the tallies across subdomains
|
||||
for tally_type, tally in avg_xs.tallies.items():
|
||||
filt_type = _DOMAIN_TO_FILTER[self.domain_type]
|
||||
tally_avg = tally.summation(filter_type=filt_type,
|
||||
filter_bins=subdomains)
|
||||
avg_xs.tallies[tally_type] = tally_avg
|
||||
|
||||
# Average each of the tallies across subdomains
|
||||
for tally_type, tally in avg_xs.tallies.items():
|
||||
filt_type = _DOMAIN_TO_FILTER[self.domain_type]
|
||||
tally_avg = tally.average(filter_type=filt_type,
|
||||
filter_bins=subdomains)
|
||||
avg_xs.tallies[tally_type] = tally_avg
|
||||
|
||||
avg_xs._domain_type = 'avg({0})'.format(self.domain_type)
|
||||
avg_xs._domain_type = 'sum({0})'.format(self.domain_type)
|
||||
avg_xs.sparse = self.sparse
|
||||
return avg_xs
|
||||
|
||||
|
|
@ -1304,8 +1315,8 @@ class MGXS(object):
|
|||
cv.check_iterable_type('subdomains', subdomains, Integral)
|
||||
elif self.domain_type == 'distribcell':
|
||||
subdomains = np.arange(self.num_subdomains, dtype=np.int)
|
||||
elif self.domain_type == 'avg(distribcell)':
|
||||
domain_filter = self.xs_tally.find_filter('avg(distribcell)')
|
||||
elif self.domain_type == 'sum(distribcell)':
|
||||
domain_filter = self.xs_tally.find_filter('sum(distribcell)')
|
||||
subdomains = domain_filter.bins
|
||||
elif self.domain_type == 'mesh':
|
||||
xyz = [range(1, x+1) for x in self.domain.dimension]
|
||||
|
|
@ -1526,9 +1537,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)
|
||||
|
|
@ -1784,17 +1805,19 @@ class MatrixMGXS(MGXS):
|
|||
if not isinstance(subdomains, string_types):
|
||||
cv.check_iterable_type('subdomains', subdomains, Integral,
|
||||
max_depth=3)
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
subdomain_bins = []
|
||||
for subdomain in subdomains:
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
filter_bins.append((subdomain,))
|
||||
subdomain_bins.append(subdomain)
|
||||
filter_bins.append(tuple(subdomain_bins))
|
||||
|
||||
# Construct list of energy group bounds tuples for all requested groups
|
||||
if not isinstance(in_groups, string_types):
|
||||
cv.check_iterable_type('groups', in_groups, Integral)
|
||||
filters.append(openmc.EnergyFilter)
|
||||
for group in in_groups:
|
||||
filters.append(openmc.EnergyFilter)
|
||||
filter_bins.append((
|
||||
self.energy_groups.get_group_bounds(group),))
|
||||
energy_bins.append((self.energy_groups.get_group_bounds(group),))
|
||||
filter_bins.append(tuple(energy_bins))
|
||||
|
||||
# Construct list of energy group bounds tuples for all requested groups
|
||||
if not isinstance(out_groups, string_types):
|
||||
|
|
@ -3207,8 +3230,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 +3249,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 +3289,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 +3367,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 +3377,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 +3419,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 +3435,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 +3462,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 +3538,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 +3594,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 +3603,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 +3635,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
|
||||
|
|
@ -3618,16 +3706,21 @@ class ScatterMatrixXS(MatrixMGXS):
|
|||
# Construct a collection of the domain filter bins
|
||||
if not isinstance(subdomains, string_types):
|
||||
cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3)
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
subdomain_bins = []
|
||||
for subdomain in subdomains:
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
filter_bins.append((subdomain,))
|
||||
subdomain_bins.append(subdomain)
|
||||
filter_bins.append(tuple(subdomain_bins))
|
||||
|
||||
# Construct list of energy group bounds tuples for all requested groups
|
||||
if not isinstance(in_groups, string_types):
|
||||
cv.check_iterable_type('groups', in_groups, Integral)
|
||||
filters.append(openmc.EnergyFilter)
|
||||
energy_bins = []
|
||||
for group in in_groups:
|
||||
filters.append(openmc.EnergyFilter)
|
||||
filter_bins.append((self.energy_groups.get_group_bounds(group),))
|
||||
energy_bins.append(
|
||||
(self.energy_groups.get_group_bounds(group),))
|
||||
filter_bins.append(tuple(energy_bins))
|
||||
|
||||
# Construct list of energy group bounds tuples for all requested groups
|
||||
if not isinstance(out_groups, string_types):
|
||||
|
|
@ -3637,7 +3730,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 +3780,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 +3803,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 +3869,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 +3954,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
|
||||
|
|
@ -4606,16 +4733,21 @@ class Chi(MGXS):
|
|||
# Construct a collection of the domain filter bins
|
||||
if not isinstance(subdomains, string_types):
|
||||
cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3)
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
subdomain_bins = []
|
||||
for subdomain in subdomains:
|
||||
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
|
||||
filter_bins.append((subdomain,))
|
||||
subdomain_bins.append(subdomain)
|
||||
filter_bins.append(tuple(subdomain_bins))
|
||||
|
||||
# Construct list of energy group bounds tuples for all requested groups
|
||||
if not isinstance(groups, string_types):
|
||||
cv.check_iterable_type('groups', groups, Integral)
|
||||
filters.append(openmc.EnergyoutFilter)
|
||||
energy_bins = []
|
||||
for group in groups:
|
||||
filters.append(openmc.EnergyoutFilter)
|
||||
filter_bins.append((self.energy_groups.get_group_bounds(group),))
|
||||
energy_bins.append(
|
||||
(self.energy_groups.get_group_bounds(group),))
|
||||
filter_bins.append(tuple(energy_bins))
|
||||
|
||||
# If chi was computed for each nuclide in the domain
|
||||
if self.by_nuclide:
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,10 +1,8 @@
|
|||
from numbers import Integral
|
||||
import sys
|
||||
import warnings
|
||||
|
||||
from six import string_types
|
||||
|
||||
from openmc.checkvalue import check_type
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
|
||||
class Nuclide(object):
|
||||
|
|
@ -72,7 +70,7 @@ class Nuclide(object):
|
|||
|
||||
@name.setter
|
||||
def name(self, name):
|
||||
check_type('name', name, string_types)
|
||||
cv.check_type('name', name, string_types)
|
||||
self._name = name
|
||||
|
||||
if '-' in name:
|
||||
|
|
|
|||
542
openmc/plotter.py
Normal file
542
openmc/plotter.py
Normal file
|
|
@ -0,0 +1,542 @@
|
|||
from numbers import Integral, Real
|
||||
from six import string_types
|
||||
from itertools import chain
|
||||
|
||||
import numpy as np
|
||||
|
||||
import openmc.checkvalue as cv
|
||||
import openmc.data
|
||||
|
||||
# Supported keywords for xs plotting
|
||||
PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
|
||||
'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity',
|
||||
'slowing-down power', 'damage']
|
||||
|
||||
# Special MT values
|
||||
UNITY_MT = -1
|
||||
XI_MT = -2
|
||||
|
||||
# MTs to combine to generate associated plot_types
|
||||
_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27]
|
||||
PLOT_TYPES_MT = {'total': openmc.data.SUM_RULES[1],
|
||||
'scatter': [2] + _INELASTIC,
|
||||
'elastic': [2],
|
||||
'inelastic': _INELASTIC,
|
||||
'fission': [18],
|
||||
'absorption': [27], 'capture': [101],
|
||||
'nu-fission': [18],
|
||||
'nu-scatter': [2] + _INELASTIC,
|
||||
'unity': [UNITY_MT],
|
||||
'slowing-down power': [2] + _INELASTIC + [XI_MT],
|
||||
'damage': [444]}
|
||||
# Operations to use when combining MTs the first np.add is used in reference
|
||||
# to zero
|
||||
PLOT_TYPES_OP = {'total': (np.add,),
|
||||
'scatter': (np.add,) * (len(PLOT_TYPES_MT['scatter']) - 1),
|
||||
'elastic': (),
|
||||
'inelastic': (np.add,) * (len(PLOT_TYPES_MT['inelastic']) - 1),
|
||||
'fission': (), 'absorption': (),
|
||||
'capture': (), 'nu-fission': (),
|
||||
'nu-scatter': (np.add,) * (len(PLOT_TYPES_MT['nu-scatter']) - 1),
|
||||
'unity': (),
|
||||
'slowing-down power':
|
||||
(np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,),
|
||||
'damage': ()}
|
||||
|
||||
# Types of plots to plot linearly in y
|
||||
PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter',
|
||||
'nu-fission / absorption', 'fission / absorption'}
|
||||
|
||||
|
||||
def plot_xs(this, types, divisor_types=None, temperature=294., axis=None,
|
||||
sab_name=None, cross_sections=None, enrichment=None, **kwargs):
|
||||
"""Creates a figure of continuous-energy cross sections for this item
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : openmc.Element, openmc.Nuclide, or openmc.Material
|
||||
Object to source data from
|
||||
types : Iterable of values of PLOT_TYPES
|
||||
The type of cross sections to include in the plot.
|
||||
divisor_types : Iterable of values of PLOT_TYPES, optional
|
||||
Cross section types which will divide those produced by types
|
||||
before plotting. A type of 'unity' can be used to effectively not
|
||||
divide some types.
|
||||
temperature : float, optional
|
||||
Temperature in Kelvin to plot. If not specified, a default
|
||||
temperature of 294K will be plotted. Note that the nearest
|
||||
temperature in the library for each nuclide will be used as opposed
|
||||
to using any interpolation.
|
||||
axis : matplotlib.axes, optional
|
||||
A previously generated axis to use for plotting. If not specified,
|
||||
a new axis and figure will be generated.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable; only used
|
||||
for items which are instances of openmc.Element or openmc.Nuclide
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
enrichment : float, optional
|
||||
Enrichment for U235 in weight percent. For example, input 4.95 for
|
||||
4.95 weight percent enriched U. Default is None. This is only used for
|
||||
items which are instances of openmc.Element
|
||||
**kwargs
|
||||
All keyword arguments are passed to
|
||||
:func:`matplotlib.pyplot.figure`.
|
||||
|
||||
Returns
|
||||
-------
|
||||
fig : matplotlib.figure.Figure
|
||||
If axis is None, then a Matplotlib Figure of the generated
|
||||
cross section will be returned. Otherwise, a value of
|
||||
None will be returned as the figure and axes have already been
|
||||
generated.
|
||||
|
||||
"""
|
||||
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
if isinstance(this, openmc.Nuclide):
|
||||
data_type = 'nuclide'
|
||||
elif isinstance(this, openmc.Element):
|
||||
data_type = 'element'
|
||||
elif isinstance(this, openmc.Material):
|
||||
data_type = 'material'
|
||||
else:
|
||||
raise TypeError("Invalid type for plotting")
|
||||
|
||||
E, data = calculate_xs(this, types, temperature, sab_name, cross_sections,
|
||||
enrichment)
|
||||
|
||||
if divisor_types:
|
||||
cv.check_length('divisor types', divisor_types, len(types),
|
||||
len(types))
|
||||
Ediv, data_div = calculate_xs(this, divisor_types, temperature,
|
||||
sab_name, cross_sections, enrichment)
|
||||
|
||||
# Create a new union grid, interpolate data and data_div on to that
|
||||
# grid, and then do the actual division
|
||||
Enum = E[:]
|
||||
E = np.union1d(Enum, Ediv)
|
||||
data_new = np.zeros((len(types), len(E)))
|
||||
|
||||
for line in range(len(types)):
|
||||
data_new[line, :] = \
|
||||
np.divide(np.interp(E, Enum, data[line, :]),
|
||||
np.interp(E, Ediv, data_div[line, :]))
|
||||
if divisor_types[line] != 'unity':
|
||||
types[line] = types[line] + ' / ' + divisor_types[line]
|
||||
data = data_new
|
||||
|
||||
# Generate the plot
|
||||
if axis is None:
|
||||
fig = plt.figure(**kwargs)
|
||||
ax = fig.add_subplot(111)
|
||||
else:
|
||||
fig = None
|
||||
ax = axis
|
||||
# Set to loglog or semilogx depending on if we are plotting a data
|
||||
# type which we expect to vary linearly
|
||||
if set(types).issubset(PLOT_TYPES_LINEAR):
|
||||
plot_func = ax.semilogx
|
||||
else:
|
||||
plot_func = ax.loglog
|
||||
# Plot the data
|
||||
for i in range(len(data)):
|
||||
data[i, :] = np.nan_to_num(data[i, :])
|
||||
if np.sum(data[i, :]) > 0.:
|
||||
plot_func(E, data[i, :], label=types[i])
|
||||
|
||||
ax.set_xlabel('Energy [eV]')
|
||||
if divisor_types:
|
||||
if data_type == 'nuclide':
|
||||
ylabel = 'Nuclidic Microscopic Data'
|
||||
elif data_type == 'element':
|
||||
ylabel = 'Elemental Microscopic Data'
|
||||
elif data_type == 'material':
|
||||
ylabel = 'Macroscopic Data'
|
||||
else:
|
||||
if data_type == 'nuclide':
|
||||
ylabel = 'Microscopic Cross Section [b]'
|
||||
elif data_type == 'element':
|
||||
ylabel = 'Elemental Cross Section [b]'
|
||||
elif data_type == 'material':
|
||||
ylabel = 'Macroscopic Cross Section [1/cm]'
|
||||
ax.set_ylabel(ylabel)
|
||||
ax.legend(loc='best')
|
||||
# Set to the most likely expected range
|
||||
ax.set_xlim((1.E-5, 20.E6))
|
||||
if this.name is not None:
|
||||
ax.set_title('Cross Section for ' + this.name)
|
||||
|
||||
return fig
|
||||
|
||||
|
||||
def calculate_xs(this, types, temperature=294., sab_name=None,
|
||||
cross_sections=None, enrichment=None):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : openmc.Element, openmc.Nuclide, or openmc.Material
|
||||
Object to source data from
|
||||
types : Iterable of values of PLOT_TYPES
|
||||
The type of cross sections to calculate
|
||||
temperature : float, optional
|
||||
Temperature in Kelvin to plot. If not specified, a default
|
||||
temperature of 294K will be plotted. Note that the nearest
|
||||
temperature in the library for each nuclide will be used as opposed
|
||||
to using any interpolation.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
enrichment : float, optional
|
||||
Enrichment for U235 in weight percent. For example, input 4.95 for
|
||||
4.95 weight percent enriched U. Default is None
|
||||
(natural composition).
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy_grid : numpy.ndarray
|
||||
Energies at which cross sections are calculated, in units of eV
|
||||
data : numpy.ndarray
|
||||
Cross sections calculated at the energy grid described by energy_grid
|
||||
|
||||
"""
|
||||
|
||||
# Check types
|
||||
cv.check_type('temperature', temperature, Real)
|
||||
if sab_name:
|
||||
cv.check_type('sab_name', sab_name, string_types)
|
||||
if enrichment:
|
||||
cv.check_type('enrichment', enrichment, Real)
|
||||
|
||||
if isinstance(this, openmc.Nuclide):
|
||||
energy_grid, xs = _calculate_xs_nuclide(this, types, temperature,
|
||||
sab_name, cross_sections)
|
||||
# Convert xs (Iterable of Callable) to a grid of cross section values
|
||||
# calculated on @ the points in energy_grid for consistency with the
|
||||
# element and material functions.
|
||||
data = np.zeros((len(types), len(energy_grid)))
|
||||
for line in range(len(types)):
|
||||
data[line, :] = xs[line](energy_grid)
|
||||
elif isinstance(this, openmc.Element):
|
||||
energy_grid, data = _calculate_xs_elem_mat(this, types, temperature,
|
||||
cross_sections, sab_name,
|
||||
enrichment)
|
||||
elif isinstance(this, openmc.Material):
|
||||
energy_grid, data = _calculate_xs_elem_mat(this, types, temperature,
|
||||
cross_sections)
|
||||
else:
|
||||
raise TypeError("Invalid type")
|
||||
|
||||
return energy_grid, data
|
||||
|
||||
|
||||
def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None,
|
||||
cross_sections=None):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : openmc.Nuclide
|
||||
Nuclide object to source data from
|
||||
types : Iterable of str or Integral
|
||||
The type of cross sections to calculate; values can either be those
|
||||
in openmc.PLOT_TYPES or integers which correspond to reaction
|
||||
channel (MT) numbers.
|
||||
temperature : float, optional
|
||||
Temperature in Kelvin to plot. If not specified, a default
|
||||
temperature of 294K will be plotted. Note that the nearest
|
||||
temperature in the library for each nuclide will be used as opposed
|
||||
to using any interpolation.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy_grid : numpy.ndarray
|
||||
Energies at which cross sections are calculated, in units of eV
|
||||
data : Iterable of Callable
|
||||
Requested cross section functions
|
||||
|
||||
"""
|
||||
|
||||
# Parse the types
|
||||
mts = []
|
||||
ops = []
|
||||
yields = []
|
||||
for line in types:
|
||||
if line in PLOT_TYPES:
|
||||
mts.append(PLOT_TYPES_MT[line])
|
||||
if line.startswith('nu'):
|
||||
yields.append(True)
|
||||
else:
|
||||
yields.append(False)
|
||||
ops.append(PLOT_TYPES_OP[line])
|
||||
else:
|
||||
# Not a built-in type, we have to parse it ourselves
|
||||
cv.check_type('MT in types', line, Integral)
|
||||
cv.check_greater_than('MT in types', line, 0)
|
||||
mts.append((line,))
|
||||
ops.append(())
|
||||
yields.append(False)
|
||||
|
||||
# Load the library
|
||||
library = openmc.data.DataLibrary.from_xml(cross_sections)
|
||||
|
||||
# Convert temperature to format needed for access in the library
|
||||
strT = "{}K".format(int(round(temperature)))
|
||||
T = temperature
|
||||
|
||||
# Now we can create the data sets to be plotted
|
||||
energy_grid = []
|
||||
xs = []
|
||||
lib = library.get_by_material(this.name)
|
||||
if lib is not None:
|
||||
nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path'])
|
||||
# Obtain the nearest temperature
|
||||
if strT in nuc.temperatures:
|
||||
nucT = strT
|
||||
else:
|
||||
data_Ts = nuc.temperatures
|
||||
for t in range(len(data_Ts)):
|
||||
# Take off the "K" and convert to a float
|
||||
data_Ts[t] = float(data_Ts[t][:-1])
|
||||
min_delta = float('inf')
|
||||
closest_t = -1
|
||||
for t in data_Ts:
|
||||
if abs(data_Ts[t] - T) < min_delta:
|
||||
closest_t = t
|
||||
nucT = "{}K".format(int(round(data_Ts[closest_t])))
|
||||
|
||||
# Prep S(a,b) data if needed
|
||||
if sab_name:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
|
||||
# Obtain the nearest temperature
|
||||
if strT in sab.temperatures:
|
||||
sabT = strT
|
||||
else:
|
||||
data_Ts = sab.temperatures
|
||||
for t in range(len(data_Ts)):
|
||||
# Take off the "K" and convert to a float
|
||||
data_Ts[t] = float(data_Ts[t][:-1])
|
||||
min_delta = np.finfo(np.float64).max
|
||||
closest_t = -1
|
||||
for t in data_Ts:
|
||||
if abs(data_Ts[t] - T) < min_delta:
|
||||
closest_t = t
|
||||
sabT = "{}K".format(int(round(data_Ts[closest_t])))
|
||||
|
||||
# Create an energy grid composed the S(a,b) and
|
||||
# the nuclide's grid
|
||||
grid = nuc.energy[nucT]
|
||||
sab_Emax = 0.
|
||||
sab_funcs = []
|
||||
if sab.elastic_xs:
|
||||
elastic = sab.elastic_xs[sabT]
|
||||
if isinstance(elastic, openmc.data.CoherentElastic):
|
||||
grid = np.union1d(grid, elastic.bragg_edges)
|
||||
if elastic.bragg_edges[-1] > sab_Emax:
|
||||
sab_Emax = elastic.bragg_edges[-1]
|
||||
elif isinstance(elastic, openmc.data.Tabulated1D):
|
||||
grid = np.union1d(grid, elastic.x)
|
||||
if elastic.x[-1] > sab_Emax:
|
||||
sab_Emax = elastic.x[-1]
|
||||
sab_funcs.append(elastic)
|
||||
if sab.inelastic_xs:
|
||||
inelastic = sab.inelastic_xs[sabT]
|
||||
grid = np.union1d(grid, inelastic.x)
|
||||
if inelastic.x[-1] > sab_Emax:
|
||||
sab_Emax = inelastic.x[-1]
|
||||
sab_funcs.append(inelastic)
|
||||
energy_grid = grid
|
||||
else:
|
||||
energy_grid = nuc.energy[nucT]
|
||||
|
||||
for i, mt_set in enumerate(mts):
|
||||
# Get the reaction xs data from the nuclide
|
||||
funcs = []
|
||||
op = ops[i]
|
||||
for mt in mt_set:
|
||||
if mt == 2:
|
||||
if sab_name:
|
||||
# Then we need to do a piece-wise function of
|
||||
# The S(a,b) and non-thermal data
|
||||
sab_sum = openmc.data.Sum(sab_funcs)
|
||||
pw_funcs = openmc.data.Regions1D(
|
||||
[sab_sum, nuc[mt].xs[nucT]],
|
||||
[sab_Emax])
|
||||
funcs.append(pw_funcs)
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt in nuc:
|
||||
if yields[i]:
|
||||
# Get the total yield first if available. This will be
|
||||
# used primarily for fission.
|
||||
for prod in chain(nuc[mt].products,
|
||||
nuc[mt].derived_products):
|
||||
if prod.particle == 'neutron' and \
|
||||
prod.emission_mode == 'total':
|
||||
func = openmc.data.Combination(
|
||||
[nuc[mt].xs[nucT], prod.yield_],
|
||||
[np.multiply])
|
||||
funcs.append(func)
|
||||
break
|
||||
else:
|
||||
# Total doesn't exist so we have to create from
|
||||
# prompt and delayed. This is used for scatter
|
||||
# multiplication.
|
||||
func = None
|
||||
for prod in chain(nuc[mt].products,
|
||||
nuc[mt].derived_products):
|
||||
if prod.particle == 'neutron' and \
|
||||
prod.emission_mode != 'total':
|
||||
if func:
|
||||
func = openmc.data.Combination(
|
||||
[prod.yield_, func], [np.add])
|
||||
else:
|
||||
func = prod.yield_
|
||||
if func:
|
||||
funcs.append(openmc.data.Combination(
|
||||
[func, nuc[mt].xs[nucT]], [np.multiply]))
|
||||
else:
|
||||
# If func is still None, then there were no
|
||||
# products. In that case, assume the yield is
|
||||
# one as its not provided for some summed
|
||||
# reactions like MT=4
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt == UNITY_MT:
|
||||
funcs.append(lambda x: 1.)
|
||||
elif mt == XI_MT:
|
||||
awr = nuc.atomic_weight_ratio
|
||||
alpha = ((awr - 1.) / (awr + 1.))**2
|
||||
xi = 1. + alpha * np.log(alpha) / (1. - alpha)
|
||||
funcs.append(lambda x: xi)
|
||||
else:
|
||||
funcs.append(lambda x: 0.)
|
||||
xs.append(openmc.data.Combination(funcs, op))
|
||||
else:
|
||||
raise ValueError(this.name + " not in library")
|
||||
|
||||
return energy_grid, xs
|
||||
|
||||
|
||||
def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None,
|
||||
sab_name=None, enrichment=None):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : {openmc.Material, openmc.Element}
|
||||
Object to source data from
|
||||
types : Iterable of values of PLOT_TYPES
|
||||
The type of cross sections to calculate
|
||||
temperature : float, optional
|
||||
Temperature in Kelvin to plot. If not specified, a default
|
||||
temperature of 294K will be plotted. Note that the nearest
|
||||
temperature in the library for each nuclide will be used as opposed
|
||||
to using any interpolation.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable.
|
||||
enrichment : float, optional
|
||||
Enrichment for U235 in weight percent. For example, input 4.95 for
|
||||
4.95 weight percent enriched U. Default is None
|
||||
(natural composition).
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy_grid : numpy.ndarray
|
||||
Energies at which cross sections are calculated, in units of eV
|
||||
data : numpy.ndarray
|
||||
Cross sections calculated at the energy grid described by energy_grid
|
||||
|
||||
"""
|
||||
|
||||
if isinstance(this, openmc.Material):
|
||||
if this.temperature is not None:
|
||||
T = this.temperature
|
||||
else:
|
||||
T = temperature
|
||||
else:
|
||||
T = temperature
|
||||
|
||||
# Load the library
|
||||
library = openmc.data.DataLibrary.from_xml(cross_sections)
|
||||
|
||||
if isinstance(this, openmc.Material):
|
||||
# Expand elements in to nuclides with atomic densities
|
||||
nuclides = this.get_nuclide_atom_densities()
|
||||
# For ease of processing split out the nuclide and its fraction
|
||||
nuc_fractions = {nuclide[1][0].name: nuclide[1][1]
|
||||
for nuclide in nuclides.items()}
|
||||
# Create a dict of [nuclide name] = nuclide object to carry forward
|
||||
# with a common nuclides format between openmc.Material and
|
||||
# openmc.Element objects
|
||||
nuclides = {nuclide[1][0].name: nuclide[1][0]
|
||||
for nuclide in nuclides.items()}
|
||||
else:
|
||||
# Expand elements in to nuclides with atomic densities
|
||||
nuclides = this.expand(1., 'ao', enrichment=enrichment,
|
||||
cross_sections=cross_sections)
|
||||
# For ease of processing split out the nuclide and its fraction
|
||||
nuc_fractions = {nuclide[0].name: nuclide[1] for nuclide in nuclides}
|
||||
# Create a dict of [nuclide name] = nuclide object to carry forward
|
||||
# with a common nuclides format between openmc.Material and
|
||||
# openmc.Element objects
|
||||
nuclides = {nuclide[0].name: nuclide[0] for nuclide in nuclides}
|
||||
|
||||
# Identify the nuclides which have S(a,b) data
|
||||
sabs = {}
|
||||
for nuclide in nuclides.items():
|
||||
sabs[nuclide[0]] = None
|
||||
if isinstance(this, openmc.Material):
|
||||
for sab_name in this._sab:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(
|
||||
library.get_by_material(sab_name)['path'])
|
||||
for nuc in sab.nuclides:
|
||||
sabs[nuc] = library.get_by_material(sab_name)['path']
|
||||
else:
|
||||
if sab_name:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
|
||||
for nuc in sab.nuclides:
|
||||
sabs[nuc] = library.get_by_material(sab_name)['path']
|
||||
|
||||
# Now we can create the data sets to be plotted
|
||||
xs = {}
|
||||
E = []
|
||||
for nuclide in nuclides.items():
|
||||
name = nuclide[0]
|
||||
nuc = nuclide[1]
|
||||
sab_tab = sabs[name]
|
||||
temp_E, temp_xs = calculate_xs(nuc, types, T, sab_tab, cross_sections)
|
||||
E.append(temp_E)
|
||||
# Since the energy grids are different, store the cross sections as
|
||||
# a tabulated function so they can be calculated on any grid needed.
|
||||
xs[name] = [openmc.data.Tabulated1D(temp_E, temp_xs[line])
|
||||
for line in range(len(types))]
|
||||
|
||||
# Condense the data for every nuclide
|
||||
# First create a union energy grid
|
||||
energy_grid = E[0]
|
||||
for grid in E[1:]:
|
||||
energy_grid = np.union1d(energy_grid, grid)
|
||||
|
||||
# Now we can combine all the nuclidic data
|
||||
data = np.zeros((len(types), len(energy_grid)))
|
||||
for line in range(len(types)):
|
||||
if types[line] == 'unity':
|
||||
data[line, :] = 1.
|
||||
else:
|
||||
for nuclide in nuclides.items():
|
||||
name = nuclide[0]
|
||||
data[line, :] += (nuc_fractions[name] *
|
||||
xs[name][line](energy_grid))
|
||||
|
||||
return energy_grid, data
|
||||
|
|
@ -387,7 +387,7 @@ class Summary(object):
|
|||
|
||||
# Set the distribcell offsets for the lattice
|
||||
if offsets is not None:
|
||||
lattice.offsets = offsets[:, ::-1, :]
|
||||
lattice.offsets = offsets
|
||||
|
||||
# Add the Lattice to the global dictionary of all Lattices
|
||||
self.lattices[index] = lattice
|
||||
|
|
|
|||
|
|
@ -1317,14 +1317,15 @@ 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]))
|
||||
|
||||
# Create list of cell instance IDs for distribcell Filters
|
||||
elif isinstance(self_filter, openmc.DistribcellFilter):
|
||||
bins = np.arange(self_filter.num_bins)
|
||||
bins = [b for b in range(self_filter.num_bins)]
|
||||
|
||||
# EnergyFunctionFilters don't have bins so just add a None
|
||||
elif isinstance(self_filter, openmc.EnergyFunctionFilter):
|
||||
|
|
@ -2119,14 +2120,14 @@ class Tally(object):
|
|||
# Construct lists of tuples for the bins in each of the two filters
|
||||
filters = [type(filter1), type(filter2)]
|
||||
if isinstance(filter1, openmc.DistribcellFilter):
|
||||
filter1_bins = np.arange(filter1.num_bins)
|
||||
filter1_bins = [b for b in range(filter1.num_bins)]
|
||||
elif isinstance(filter1, openmc.EnergyFunctionFilter):
|
||||
filter1_bins = [None]
|
||||
else:
|
||||
filter1_bins = [filter1.get_bin(i) for i in range(filter1.num_bins)]
|
||||
|
||||
if isinstance(filter2, openmc.DistribcellFilter):
|
||||
filter2_bins = np.arange(filter2.num_bins)
|
||||
filter2_bins = [b for b in range(filter2.num_bins)]
|
||||
elif isinstance(filter2, openmc.EnergyFunctionFilter):
|
||||
filter2_bins = [None]
|
||||
else:
|
||||
|
|
@ -2975,8 +2976,16 @@ class Tally(object):
|
|||
# Sum across the bins in the user-specified filter
|
||||
for i, self_filter in enumerate(self.filters):
|
||||
if isinstance(self_filter, filter_type):
|
||||
shape = mean.shape
|
||||
mean = np.take(mean, indices=bin_indices, axis=i)
|
||||
std_dev = np.take(std_dev, indices=bin_indices, axis=i)
|
||||
|
||||
# NumPy take introduces a new dimension in output array
|
||||
# for some special cases that must be removed
|
||||
if len(mean.shape) > len(shape):
|
||||
mean = np.squeeze(mean, axis=i)
|
||||
std_dev = np.squeeze(std_dev, axis=i)
|
||||
|
||||
mean = np.sum(mean, axis=i, keepdims=True)
|
||||
std_dev = np.sum(std_dev**2, axis=i, keepdims=True)
|
||||
std_dev = np.sqrt(std_dev)
|
||||
|
|
@ -3124,10 +3133,18 @@ class Tally(object):
|
|||
# Average across the bins in the user-specified filter
|
||||
for i, self_filter in enumerate(self.filters):
|
||||
if isinstance(self_filter, filter_type):
|
||||
shape = mean.shape
|
||||
mean = np.take(mean, indices=bin_indices, axis=i)
|
||||
std_dev = np.take(std_dev, indices=bin_indices, axis=i)
|
||||
mean = np.mean(mean, axis=i, keepdims=True)
|
||||
std_dev = np.mean(std_dev**2, axis=i, keepdims=True)
|
||||
|
||||
# NumPy take introduces a new dimension in output array
|
||||
# for some special cases that must be removed
|
||||
if len(mean.shape) > len(shape):
|
||||
mean = np.squeeze(mean, axis=i)
|
||||
std_dev = np.squeeze(std_dev, axis=i)
|
||||
|
||||
mean = np.nanmean(mean, axis=i, keepdims=True)
|
||||
std_dev = np.nanmean(std_dev**2, axis=i, keepdims=True)
|
||||
std_dev /= len(bin_indices)
|
||||
std_dev = np.sqrt(std_dev)
|
||||
|
||||
|
|
@ -3151,8 +3168,8 @@ class Tally(object):
|
|||
axis_index = self.num_filters
|
||||
mean = np.take(mean, indices=nuclide_bins, axis=axis_index)
|
||||
std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index)
|
||||
mean = np.mean(mean, axis=axis_index, keepdims=True)
|
||||
std_dev = np.mean(std_dev**2, axis=axis_index, keepdims=True)
|
||||
mean = np.nanmean(mean, axis=axis_index, keepdims=True)
|
||||
std_dev = np.nanmean(std_dev**2, axis=axis_index, keepdims=True)
|
||||
std_dev /= len(nuclide_bins)
|
||||
std_dev = np.sqrt(std_dev)
|
||||
|
||||
|
|
@ -3170,8 +3187,8 @@ class Tally(object):
|
|||
axis_index = self.num_filters + 1
|
||||
mean = np.take(mean, indices=score_bins, axis=axis_index)
|
||||
std_dev = np.take(std_dev, indices=score_bins, axis=axis_index)
|
||||
mean = np.sum(mean, axis=axis_index, keepdims=True)
|
||||
std_dev = np.sum(std_dev**2, axis=axis_index, keepdims=True)
|
||||
mean = np.nanmean(mean, axis=axis_index, keepdims=True)
|
||||
std_dev = np.nanmean(std_dev**2, axis=axis_index, keepdims=True)
|
||||
std_dev /= len(score_bins)
|
||||
std_dev = np.sqrt(std_dev)
|
||||
|
||||
|
|
|
|||
|
|
@ -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
1
setup.py
|
|
@ -43,6 +43,7 @@ if have_setuptools:
|
|||
|
||||
# Optional dependencies
|
||||
'extras_require': {
|
||||
'decay': ['uncertainties'],
|
||||
'pandas': ['pandas>=0.17.0'],
|
||||
'sparse' : ['scipy'],
|
||||
'vtk': ['vtk', 'silomesh'],
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
BIN
tests/1d_mgxs.h5
Binary file not shown.
|
|
@ -1 +1 @@
|
|||
e18c2318bab6c42a263e5079fd796b1bee609e4274884fbc7bfd8b33e59aeb5e5a167da8e093f339f766815e007bd606c002939e3f730af523cbc9cf75c53faa
|
||||
b70886031e22db9e3f0332eac703a7356504750c1e90d7083ffd16b8884d00661d0e20c6d8bead3c93369b2e7c105ca3280c7858ca6a147fa6669a5d3d530461
|
||||
|
|
@ -1 +1 @@
|
|||
3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989
|
||||
864328b2c4f3c4bfa9c80c756acbedac6fbdf3d502c03c2f2a3e7461b774729c47a55341c9515eff246e1bd445485f0d0a05a5712663ee499046afdbe0c77aef
|
||||
|
|
@ -1,79 +1,79 @@
|
|||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.457353 0.010474
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.405649 0.015784
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.405641 0.015787
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.066556 0.00251
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.028979 0.002712
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.037577 0.001487
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.092377 0.003628
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 7.276707e+06 287579.26286
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.390797 0.008717
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.387332 0.014241
|
||||
avg(distribcell) group in group out nuclide moment mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P0 0.387009 0.014230
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P1 0.047179 0.004923
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P2 0.015713 0.003654
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P3 0.005378 0.003137
|
||||
avg(distribcell) group in group out nuclide moment mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P0 0.387332 0.014241
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P1 0.047187 0.004933
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P2 0.015727 0.003654
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P3 0.005387 0.003141
|
||||
avg(distribcell) group in group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 1.000834 0.037242
|
||||
avg(distribcell) group in group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.094516 0.0059
|
||||
avg(distribcell) group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 1.0 0.080455
|
||||
avg(distribcell) group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 1.0 0.080541
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 5.139437e-07 2.133314e-08
|
||||
avg(distribcell) group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.091725 0.003604
|
||||
avg(distribcell) group in group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.093985 0.005872
|
||||
avg(distribcell) delayedgroup group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000021 8.253907e-07
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.000112 4.284000e-06
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 0.000109 4.105197e-06
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.000252 9.271420e-06
|
||||
4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.000112 3.888625e-06
|
||||
5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 0.000047 1.625563e-06
|
||||
avg(distribcell) delayedgroup group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.0 0.000000
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 1.0 1.414214
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 1.0 1.414214
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.0 0.000000
|
||||
4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.0 0.000000
|
||||
5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 1.0 1.414214
|
||||
avg(distribcell) delayedgroup group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000227 0.000012
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.001209 0.000061
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 0.001177 0.000059
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.002727 0.000135
|
||||
4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.001210 0.000058
|
||||
5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 0.000504 0.000024
|
||||
avg(distribcell) delayedgroup group in nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000000 0.000000
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.032739 0.046300
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 0.120780 0.170809
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.000000 0.000000
|
||||
4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.000000 0.000000
|
||||
5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 2.853000 4.034751
|
||||
avg(distribcell) delayedgroup group in group out nuclide mean std. dev.
|
||||
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 1 total 0.000000 0.000000
|
||||
1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 1 total 0.000175 0.000175
|
||||
2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 1 total 0.000178 0.000178
|
||||
3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 1 total 0.000000 0.000000
|
||||
4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 1 total 0.000000 0.000000
|
||||
5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 1 total 0.000178 0.000178
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.457353 0.010474
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.405649 0.015784
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.405641 0.015787
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.066556 0.00251
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.028979 0.002712
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.037577 0.001487
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.092377 0.003628
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 7.276707e+06 287579.26286
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.390797 0.008717
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.387332 0.014241
|
||||
sum(distribcell) group in group out nuclide moment mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P0 0.387009 0.014230
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P1 0.047179 0.004923
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P2 0.015713 0.003654
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P3 0.005378 0.003137
|
||||
sum(distribcell) group in group out nuclide moment mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P0 0.387332 0.014241
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P1 0.047187 0.004933
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P2 0.015727 0.003654
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P3 0.005387 0.003141
|
||||
sum(distribcell) group in group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 1.000834 0.037242
|
||||
sum(distribcell) group in group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.094516 0.0059
|
||||
sum(distribcell) group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 1.0 0.080455
|
||||
sum(distribcell) group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 1.0 0.080541
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 5.139437e-07 2.133314e-08
|
||||
sum(distribcell) group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.091725 0.003604
|
||||
sum(distribcell) group in group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.093985 0.005872
|
||||
sum(distribcell) delayedgroup group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.000021 8.253907e-07
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 0.000112 4.284000e-06
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 0.000109 4.105197e-06
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.000252 9.271420e-06
|
||||
4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.000112 3.888625e-06
|
||||
5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 0.000047 1.625563e-06
|
||||
sum(distribcell) delayedgroup group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.0 0.000000
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 1.0 1.414214
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 1.0 1.414214
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.0 0.000000
|
||||
4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.0 0.000000
|
||||
5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 1.0 1.414214
|
||||
sum(distribcell) delayedgroup group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.000227 0.000012
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 0.001209 0.000061
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 0.001177 0.000059
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.002727 0.000135
|
||||
4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.001210 0.000058
|
||||
5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 0.000504 0.000024
|
||||
sum(distribcell) delayedgroup group in nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.000000 0.000000
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 0.032739 0.046300
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 0.120780 0.170809
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.000000 0.000000
|
||||
4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.000000 0.000000
|
||||
5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 2.853000 4.034751
|
||||
sum(distribcell) delayedgroup group in group out nuclide mean std. dev.
|
||||
0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 1 total 0.000000 0.000000
|
||||
1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 1 total 0.000175 0.000175
|
||||
2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 1 total 0.000178 0.000178
|
||||
3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 1 total 0.000000 0.000000
|
||||
4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 1 total 0.000000 0.000000
|
||||
5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 1 total 0.000178 0.000178
|
||||
|
|
|
|||
|
|
@ -1 +1 @@
|
|||
b1f616bc0e342e3886b4b27eca22c4d1e55e445979254d5ba5f9ea01048f53d0e9a01c575745ae130eb34108080d2256ab9e4cb926ea0050891a44ae37ecdd88
|
||||
ac118c5796593df0f53491920875bbe77921882253e7f24942e365e5812d726fb3d3bab86da02b3aad93f6f8f20a090125ffa51f9fdf9aef84009702f9cc363d
|
||||
|
|
@ -1 +1 @@
|
|||
ceb1420d48e097185dd0798a495676c8e968c37c6e75fd6232137a4df8f6fae33d8232b95c68f6479815a014d20ccfc7611f96b22fcb7da159b42346e00cdf4e
|
||||
a8172f7492fcc69d21f53c20523ffd0ad2882a753a3e1d263d135153cc6b3ee8d8e38f8823a4a2eb3288d730579892875d2d19f85f083dd659f91414988115ea
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue