crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

fixing merged conflicts in cmakelist and simulation

Browse source 
Merge branch 'develop' of https://github.com/mit-crpg/openmc into new-update

Conflicts:
	CMakeLists.txt
	src/simulation.F90
This commit is contained in:
Jose Salcedo Perez 2018-03-29 01:10:11 +00:00
commit 5230703779
710 changed files with 13821 additions and 9742 deletions
Download patch file Download diff file Expand all files Collapse all files

9
.gitignore vendored
View file

@ -42,11 +42,8 @@ results_error.dat
inputs_error.dat
results_test.dat
# Test build files
tests/build/
tests/coverage/
tests/memcheck/
tests/ctestscript.run
# Test
.pytest_cache/
# HDF5 files
*.h5
@ -101,3 +98,5 @@ examples/jupyter/plots
.cache/
.tox/
.python-version
.coverage
htmlcov

21
.travis.yml
View file

@ -2,8 +2,9 @@ sudo: required
dist: trusty
language: python
python:
- "2.7"
- "3.4"
- "3.5"
- "3.6"
addons:
apt:
packages:
@ -18,29 +19,27 @@ env:
global:
- FC=gfortran
- MPI_DIR=/usr
- PHDF5_DIR=/usr
- HDF5_DIR=/usr
- HDF5_ROOT=/usr
- OMP_NUM_THREADS=2
- OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
- OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
- OPENMC_MULTIPOLE_LIBRARY=$HOME/multipole_lib
- PATH=$PATH:$HOME/NJOY2016/build
- DISPLAY=:99.0
matrix:
- OPENMC_CONFIG="^hdf5-debug$"
- OPENMC_CONFIG="^omp-hdf5-debug$"
- OPENMC_CONFIG="^mpi-hdf5-debug$"
- OPENMC_CONFIG="^phdf5-debug$"
- OMP=n MPI=n PHDF5=n
- OMP=y MPI=n PHDF5=n
- OMP=n MPI=y PHDF5=n
- OMP=n MPI=y PHDF5=y
before_install:
- sudo add-apt-repository ppa:nschloe/hdf5-backports -y
- sudo apt-get update -q
- sudo apt-get install libhdf5-serial-dev libhdf5-mpich-dev -y
install:
- ./tools/ci/travis-install.sh
before_script:
- ./tools/ci/travis-before-script.sh
script:
- ./tools/ci/travis-script.sh
after_success:
- coveralls

103
CMakeLists.txt
View file

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 2.8.12 FATAL_ERROR)
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
project(openmc Fortran C CXX)
# Setup output directories
@ -21,11 +21,6 @@ if (${UNIX})
add_definitions(-DUNIX)
endif()
# Set MACOSX_RPATH
if(POLICY CMP0042)
cmake_policy(SET CMP0042 NEW)
endif()
#===============================================================================
# Command line options
#===============================================================================
@ -48,14 +43,17 @@ add_definitions(-DMAX_COORD=${maxcoord})
set(MPI_ENABLED FALSE)
if($ENV{FC} MATCHES "(mpi[^/]*|ftn)$")
message("-- Detected MPI wrapper: $ENV{FC}")
add_definitions(-DMPI)
add_definitions(-DOPENMC_MPI)
set(MPI_ENABLED TRUE)
# Get directory containing MPI wrapper
get_filename_component(MPI_DIR $ENV{FC} DIRECTORY)
endif()
# Check for Fortran 2008 MPI interface
if(MPI_ENABLED AND mpif08)
message("-- Using Fortran 2008 MPI bindings")
add_definitions(-DMPIF08)
add_definitions(-DOPENMC_MPIF08)
endif()
#===============================================================================
@ -116,10 +114,17 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
endif()
# GCC compiler options
list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2 -fstack-arrays -g)
<<<<<<< HEAD
list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2 -g)
if(debug)
list(REMOVE_ITEM f90flags -O2)
list(APPEND f90flags -g -enable-checking -fbacktrace -Wall -Wno-unused-dummy-argument -pedantic
=======
list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2 -fstack-arrays)
if(debug)
list(REMOVE_ITEM f90flags -O2 -fstack-arrays)
list(APPEND f90flags -g -Wall -Wno-unused-dummy-argument -pedantic
>>>>>>> 47fbf8282ea94c138f75219bd10fdb31501d3fb7
-fbounds-check -ffpe-trap=invalid,overflow,underflow)
list(APPEND ldflags -g -v -da -Q)
endif()
@ -347,7 +352,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/cmfd_prod_operator.F90
src/cmfd_solver.F90
src/constants.F90
src/cross_section.F90
src/dict_header.F90
src/distribution_multivariate.F90
src/distribution_univariate.F90
@ -371,7 +375,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/message_passing.F90
src/mgxs_data.F90
src/mgxs_header.F90
src/multipole.F90
src/multipole_header.F90
src/nuclide_header.F90
src/output.F90
@ -426,6 +429,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/tally_filter_energyfunc.F90
src/tallies/tally_filter_material.F90
src/tallies/tally_filter_mesh.F90
src/tallies/tally_filter_meshsurface.F90
src/tallies/tally_filter_mu.F90
src/tallies/tally_filter_polar.F90
src/tallies/tally_filter_surface.F90
@ -435,8 +439,15 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/trigger_header.F90
)
set(LIBOPENMC_CXX_SRC
src/error.h
src/hdf5_interface.h
src/random_lcg.cpp
src/random_lcg.h
src/random_lcg.cpp)
src/surface.cpp
src/surface.h
src/xml_interface.h
src/pugixml/pugixml.cpp
src/pugixml/pugixml.hpp)
add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc)
add_executable(${program} src/main.F90)
@ -492,72 +503,8 @@ add_custom_command(TARGET libopenmc POST_BUILD
install(TARGETS ${program} libopenmc
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib)
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
install(DIRECTORY src/relaxng DESTINATION share/openmc)
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)
find_package(PythonInterp)
if(PYTHONINTERP_FOUND)
if(debian)
install(CODE "execute_process(
COMMAND ${PYTHON_EXECUTABLE} setup.py install
--root=debian/openmc --install-layout=deb
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})")
else()
install(CODE "set(ENV{PYTHONPATH} \"${CMAKE_INSTALL_PREFIX}/lib/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages\")")
install(CODE "execute_process(
COMMAND ${PYTHON_EXECUTABLE} setup.py install
--prefix=${CMAKE_INSTALL_PREFIX}
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})")
endif()
endif()
#===============================================================================
# Regression tests
#===============================================================================
# This allows for dashboard configuration
include(CTest)
# Get a list of all the tests to run
file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_*.py)
# Loop through all the tests
foreach(test ${TESTS})
# Remove unit tests
if(test MATCHES ".*unit_tests.*")
continue()
endif()
# Get test information
get_filename_component(TEST_NAME ${test} NAME)
get_filename_component(TEST_PATH ${test} PATH)
if (DEFINED ENV{MEM_CHECK})
# Generate input files if needed
if (NOT EXISTS "${TEST_PATH}/geometry.xml")
execute_process(COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --build-inputs
WORKING_DIRECTORY ${TEST_PATH})
endif()
# Add serial test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc>)
else()
# Check serial/parallel
if (${MPI_ENABLED})
# Preform a parallel test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $<TARGET_FILE:openmc>
--mpi_exec $ENV{MPI_DIR}/bin/mpiexec)
else()
# Perform a serial test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $<TARGET_FILE:openmc>)
endif()
endif()
endforeach(test)

21
docs/source/conf.py
View file

@ -18,20 +18,21 @@ on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
# On Read the Docs, we need to mock a few third-party modules so we don't get
# ImportErrors when building documentation
try:
from unittest.mock import MagicMock
except ImportError:
from mock import Mock as MagicMock
from unittest.mock import MagicMock
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.interpolate',
'scipy.integrate', 'scipy.optimize', 'scipy.special',
'scipy.stats', 'h5py', 'pandas', 'uncertainties', 'matplotlib',
'matplotlib.pyplot','openmoc', 'openmc.data.reconstruct']
MOCK_MODULES = [
'numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.sparse.linalg',
'scipy.interpolate', 'scipy.integrate', 'scipy.optimize', 'scipy.special',
'scipy.stats', 'scipy.spatial', 'h5py', 'pandas', 'uncertainties',
'matplotlib', 'matplotlib.pyplot', 'openmoc',
'openmc.data.reconstruct'
]
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
import numpy as np
np.ndarray = MagicMock
np.polynomial.Polynomial = MagicMock
@ -253,6 +254,6 @@ napoleon_use_ivar = True
intersphinx_mapping = {
'python': ('https://docs.python.org/3', None),
'numpy': ('https://docs.scipy.org/doc/numpy/', None),
'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None),
'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
'matplotlib': ('https://matplotlib.org/', None)
}

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

@ -10,9 +10,10 @@ as debugging.
.. toctree::
:numbered:
:maxdepth: 3
:maxdepth: 2
styleguide
workflow
tests
user-input
docbuild

73
docs/source/devguide/tests.rst Normal file
View file

@ -0,0 +1,73 @@
.. _devguide_tests:
==========
Test Suite
==========
Running Tests
-------------
The OpenMC test suite consists of two parts, a regression test suite and a unit
test suite. The regression test suite is based on regression or integrated
testing where different types of input files are configured and the full OpenMC
code is executed. Results from simulations are compared with expected
results. The unit tests are primarily intended to test individual
functions/classes in the OpenMC Python API.
The test suite relies on the third-party `pytest <https://pytest.org>`_
package. To run either or both the regression and unit test suites, it is
assumed that you have OpenMC fully installed, i.e., the :ref:`scripts_openmc`
executable is available on your :envvar:`PATH` and the :mod:`openmc` Python
module is importable. In development where it would be onerous to continually
install OpenMC every time a small change is made, it is recommended to install
OpenMC in development/editable mode. With setuptools, this is accomplished by
running::
python setup.py develop
or using pip (recommended)::
pip install -e .[test]
It is also assumed that you have cross section data available that is pointed to
by the :envvar:`OPENMC_CROSS_SECTIONS` and :envvar:`OPENMC_MULTIPOLE_LIBRARY`
environment variables. Furthermore, to run unit tests for the :mod:`openmc.data`
module, it is necessary to have ENDF/B-VII.1 data available and pointed to by
the :envvar:`OPENMC_ENDF_DATA` environment variable. All data sources can be
obtained using the ``tools/ci/travis-before-script.sh`` script.
To execute the test suite, go to the ``tests/`` directory and run::
pytest
If you want to collect information about source line coverage in the Python API,
you must have the `pytest-cov <https://pypi.python.org/pypi/pytest-cov>`_ plugin
installed and run::
pytest --cov=../openmc --cov-report=html
Adding Tests to the Regression Suite
------------------------------------
To add a new test to the regression test suite, create a sub-directory in the
``tests/regression_tests/`` directory. To configure a test you need to add the
following files to your new test directory:
* OpenMC input XML files, if they are not generated through the Python API
* **test.py** - Python test driver script; please refer to other tests to
see how to construct. Any output files that are generated during testing
must be removed at the end of this script.
* **inputs_true.dat** - ASCII file that contains Python API-generated XML
files concatenated together. When the test is run, inputs that are
generated are compared to this file.
* **results_true.dat** - ASCII file that contains the expected results from
the test. The file *results_test.dat* is compared to this file during the
execution of the python test driver script. When the above files have been
created, generate a *results_test.dat* file and copy it to this name and
commit. It should be noted that this file should be generated with basic
compiler options during openmc configuration and build (e.g., no MPI, no
debug/optimization).
In addition to this description, please see the various types of tests that are
already included in the test suite to see how to create them. If all is
implemented correctly, the new test will automatically be discovered by pytest.

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

@ -89,139 +89,6 @@ features and bug fixes. The general steps for contributing are as follows:
6. After the pull request has been thoroughly vetted, it is merged back into the
*develop* branch of mit-crpg/openmc.
.. _test suite:
OpenMC Test Suite
-----------------
The purpose of this test suite is to ensure that OpenMC compiles using various
combinations of compiler flags and options, and that all user input options can
be used successfully without breaking the code. The test suite is comprised of
regression tests where different types of input files are configured and the
full OpenMC code is executed. Results from simulations are compared with
expected results. The test suite is comprised of many build configurations
(e.g. debug, mpi, hdf5) and the actual tests which reside in sub-directories
in the tests directory. We recommend to developers to test their branches
before submitting a formal pull request using gfortran and Intel compilers
if available.
The test suite is designed to integrate with cmake using ctest_. It is
configured to run with cross sections from NNDC_ augmented with 0 K elastic
scattering data for select nuclides as well as multipole data. To download the
proper data, run the following commands:
.. code-block:: sh
wget -O nndc_hdf5.tar.xz $(cat <openmc_root>/.travis.yml | grep anl.box | awk '{print $2}')
tar xJvf nndc_hdf5.tar.xz
export OPENMC_CROSS_SECTIONS=$(pwd)/nndc_hdf5/cross_sections.xml
git clone --branch=master git://github.com/smharper/windowed_multipole_library.git wmp_lib
tar xzvf wmp_lib/multipole_lib.tar.gz
export OPENMC_MULTIPOLE_LIBRARY=$(pwd)/multipole_lib
The test suite can be run on an already existing build using:
.. code-block:: sh
cd build
make test
or
.. code-block:: sh
cd build
ctest
There are numerous ctest_ command line options that can be set to have
more control over which tests are executed.
Before running the test suite python script, the following environmental
variables should be set if the default paths are incorrect:
* **FC** - The command for a Fortran compiler (e.g. gfotran, ifort).
* Default - *gfortran*
* **CC** - The command for a C compiler (e.g. gcc, icc).
* Default - *gcc*
* **CXX** - The command for a C++ compiler (e.g. g++, icpc).
* Default - *g++*
* **MPI_DIR** - The path to the MPI directory.
* Default - */opt/mpich/3.2-gnu*
* **HDF5_DIR** - The path to the HDF5 directory.
* Default - */opt/hdf5/1.8.16-gnu*
* **PHDF5_DIR** - The path to the parallel HDF5 directory.
* Default - */opt/phdf5/1.8.16-gnu*
To run the full test suite, the following command can be executed in the
tests directory:
.. code-block:: sh
python run_tests.py
A subset of build configurations and/or tests can be run. To see how to use
the script run:
.. code-block:: sh
python run_tests.py --help
As an example, say we want to run all tests with debug flags only on tests
that have cone and plot in their name. Also, we would like to run this on
4 processors. We can run:
.. code-block:: sh
python run_tests.py -j 4 -C debug -R "cone|plot"
Note that standard regular expression syntax is used for selecting build
configurations and tests. To print out a list of build configurations, we
can run:
.. code-block:: sh
python run_tests.py -p
Adding tests to test suite
++++++++++++++++++++++++++
To add a new test to the test suite, create a sub-directory in the tests
directory that conforms to the regular expression *test_*. To configure
a test you need to add the following files to your new test directory,
*test_name* for example:
* OpenMC input XML files
* **test_name.py** - Python test driver script, please refer to other
tests to see how to construct. Any output files that are generated
during testing must be removed at the end of this script.
* **inputs_true.dat** - ASCII file that contains Python API-generated XML
files concatenated together. When the test is run, inputs that are
generated are compared to this file.
* **results_true.dat** - ASCII file that contains the expected results
from the test. The file *results_test.dat* is compared to this file
during the execution of the python test driver script. When the
above files have been created, generate a *results_test.dat* file and
copy it to this name and commit. It should be noted that this file
should be generated with basic compiler options during openmc
configuration and build (e.g., no MPI/HDF5, no debug/optimization).
In addition to this description, please see the various types of tests that
are already included in the test suite to see how to create them. If all is
implemented correctly, the new test directory will automatically be added
to the CTest framework.
Private Development
-------------------
@ -236,6 +103,27 @@ changes you've made in your private repository back to mit-crpg/openmc
repository, simply follow the steps above with an extra step of pulling a branch
from your private repository into a public fork.
.. _devguide_editable:
Working in "Development" Mode
-----------------------------
If you are making changes to the Python API during development, it is highly
suggested to install the Python API in development/editable mode using
pip_. From the root directory of the OpenMC repository, run:
.. code-block:: sh
pip install -e .[test]
This installs the OpenMC Python package in `"editable" mode
<https://pip.pypa.io/en/stable/reference/pip_install/#editable-installs>`_ so
that 1) it can be imported from a Python interpreter and 2) any changes made are
immediately reflected in the installed version (that is, you don't need to keep
reinstalling it). While the same effect can be achieved using the
:envvar:`PYTHONPATH` environment variable, this is generally discouraged as it
can interfere with virtual environments.
.. _git: http://git-scm.com/
.. _GitHub: https://github.com/
.. _git flow: http://nvie.com/git-model
@ -247,3 +135,4 @@ from your private repository into a public fork.
.. _Bitbucket: https://bitbucket.org
.. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html
.. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html
.. _pip: https://pip.pypa.io/en/stable/

25
docs/source/io_formats/data_wmp.rst
View file

@ -28,12 +28,6 @@ Windowed Multipole Library Format
":math:`r`" and ":math:`i`" identifiers, similar to how `h5py`_ does it.
- **end_E** (*double*)
Highest energy the windowed multipole part of the library is valid for.
- **energy_points** (*double[]*)
Energy grid for the pointwise library in the reaction group.
- **fissionable** (*int*)
1 if this nuclide has fission data. 0 if it does not.
- **fit_order** (*int*)
The order of the curve fit.
- **formalism** (*int*)
The formalism of the underlying data. Uses the `ENDF-6`_ format
formalism numbers.
@ -51,18 +45,6 @@ Windowed Multipole Library Format
- **l_value** (*int[]*)
The index for a corresponding pole. Equivalent to the :math:`l` quantum
number of the resonance the pole comes from :math:`+1`.
- **length** (*int*)
Total count of poles in `data`.
- **max_w** (*int*)
Maximum number of poles in a window.
- **MT_count** (*int*)
Number of pointwise tables in the library.
- **MT_list** (*int[]*)
A list of available MT identifiers. See `ENDF-6`_ for meaning.
- **n_grid** (*int*)
Total length of the pointwise data.
- **num_l** (*int*)
Number of possible :math:`l` quantum states for this nuclide.
- **pseudo_K0RS** (*double[]*)
:math:`l` dependent value of
@ -90,13 +72,6 @@ Windowed Multipole Library Format
The pole to start from for each window.
- **w_end** (*int[]*)
The pole to end at for each window.
- **windows** (*int*)
Number of windows.
**/nuclide/reactions/MT<i>**
- **MT_sigma** (*double[]*) -- Cross section value for this reaction.
- **Q_value** (*double*) -- Energy released in this reaction, in eV.
- **threshold** (*int*) -- The first non-zero entry in ``MT_sigma``.
.. _h5py: http://docs.h5py.org/en/latest/
.. _ENDF-6: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf

102
docs/source/io_formats/depletion_chain.rst Normal file
View file

@ -0,0 +1,102 @@
.. _io_depletion_chain:
============================
Depletion Chain -- chain.xml
============================
A depletion chain file has a ``<depletion_chain>`` root element with one or more
``<nuclide>`` child elements. The decay, reaction, and fission product data for
each nuclide appears as child elements of ``<nuclide>``.
---------------------
``<nuclide>`` Element
---------------------
The ``<nuclide>`` element contains information on the decay modes, reactions,
and fission product yields for a given nuclide in the depletion chain. This
element may have the following attributes:
:name:
Name of the nuclide
:half_life:
Half-life of the nuclide in [s]
:decay_modes:
Number of decay modes present
:decay_energy:
Decay energy released in [eV]
:reactions:
Number of reactions present
For each decay mode, a :ref:`io_chain_decay` appears as a child of
``<nuclide>``. For each reaction present, a :ref:`io_chain_reaction` appears as
a child of ``<nuclide>``. If the nuclide is fissionable, a :ref:`io_chain_nfy`
appears as well.
.. _io_chain_decay:
-------------------
``<decay>`` Element
-------------------
The ``<decay>`` element represents a single decay mode and has the following
attributes:
:type:
The type of the decay, e.g. 'ec/beta+'
:target:
The daughter nuclide produced from the decay
:branching_ratio:
The branching ratio for this decay mode
.. _io_chain_reaction:
----------------------
``<reaction>`` Element
----------------------
The ``<reaction>`` element represents a single transmutation reaction. This
element has the following attributes:
:type:
The type of the reaction, e.g., '(n,gamma)'
:Q:
The Q value of the reaction in [eV]
:target:
The nuclide produced in the reaction (absent if the type is 'fission')
:branching_ratio:
The branching ratio for the reaction
.. _io_chain_nfy:
------------------------------------
``<neutron_fission_yields>`` Element
------------------------------------
The ``<neutron_fission_yields>`` element provides yields of fission products for
fissionable nuclides. It has the follow sub-elements:
:energies:
Energies in [eV] at which yields for products are tabulated
:fission_yields:
Fission product yields for a single energy point. This element itself has a
number of attributes/sub-elements:
:energy:
Energy in [eV] at which yields are tabulated
:products:
Names of fission products
:data:
Independent yields for each fission product

42
docs/source/io_formats/depletion_results.rst Normal file
View file

@ -0,0 +1,42 @@
.. _io_depletion_results:
=============================
Depletion Results File Format
=============================
The current version of the depletion results file format is 1.0.
**/**
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the
statepoint file format.
:Datasets: - **eigenvalues** (*double[][]*) -- k-eigenvalues at each
time/stage. This array has shape (number of timesteps, number of
stages).
- **number** (*double[][][][]*) -- Total number of atoms. This array
has shape (number of timesteps, number of stages, number of
materials, number of nuclides).
- **reaction rates** (*double[][][][][]*) -- Reaction rates used to
build depletion matrices. This array has shape (number of
timesteps, number of stages, number of materials, number of
nuclides, number of reactions).
- **time** (*double[][2]*) -- Time in [s] at beginning/end of each
step.
**/materials/<id>/**
:Attributes: - **index** (*int*) -- Index used in results for this material
- **volume** (*double*) -- Volume of this material in [cm^3]
**/nuclides/<name>/**
:Attributes: - **atom number index** (*int*) -- Index in array of total atoms
for this nuclide
- **reaction rate index** (*int*) -- Index in array of reaction
rates for this nuclide
**/reactions/<name>/**
:Attributes: - **index** (*int*) -- Index user in results for this reaction

8
docs/source/io_formats/index.rst
View file

@ -12,7 +12,7 @@ Input Files
.. toctree::
:numbered:
:maxdepth: 2
:maxdepth: 1
geometry
materials
@ -27,9 +27,10 @@ Data Files
.. toctree::
:numbered:
:maxdepth: 2
:maxdepth: 1
cross_sections
depletion_chain
nuclear_data
mgxs_library
data_wmp
@ -41,11 +42,12 @@ Output Files
.. toctree::
:numbered:
:maxdepth: 2
:maxdepth: 1
statepoint
source
summary
depletion_results
particle_restart
track
voxel

15
docs/source/publications.rst
View file

@ -57,6 +57,11 @@ Benchmarking
Coupling and Multi-physics
--------------------------
- Jun Chen, Liangzhi Cao, Chuanqi Zhao, and Zhouyu Liu, "`Development of
Subchannel Code SUBSC for high-fidelity multi-physics coupling application
<https://doi.org/10.1016/j.egypro.2017.08.121>`_", Energy Procedia, **127**,
264-274 (2017).
- Tianliang Hu, Liangzhu Cao, Hongchun Wu, Xianan Du, and Mingtao He, "`Coupled
neutrons and thermal-hydraulics simulation of molten salt reactors based on
OpenMC/TANSY <https://doi.org/10.1016/j.anucene.2017.05.002>`_,"
@ -98,6 +103,11 @@ Coupling and Multi-physics
Geometry and Visualization
--------------------------
- Jin-Yang Li, Long Gu, Hu-Shan Xu, Nadezha Korepanova, Rui Yu, Yan-Lei Zhu, and
Chang-Ping Qin, "`CAD modeling study on FLUKA and OpenMC for accelerator
driven system simulation <https://doi.org/10.1016/j.anucene.2017.12.050>`_",
*Ann. Nucl. Energy*, **114**, 329-341 (2018).
- Logan Abel, William Boyd, Benoit Forget, and Kord Smith, "Interactive
Visualization of Multi-Group Cross Sections on High-Fidelity Spatial Meshes,"
*Trans. Am. Nucl. Soc.*, **114**, 391-394 (2016).
@ -114,6 +124,11 @@ Geometry and Visualization
Miscellaneous
-------------
- Bruno Merk, Dzianis Litskevich, R. Gregg, and A. R. Mount, "`Demand driven
salt clean-up in a molten salt fast reactor -- Defining a priority list
<https://doi.org/10.1371/journal.pone.0192020>`_", *PLOS One*, **13**,
e0192020 (2018).
- Adam G. Nelson, Samuel Shaner, William Boyd, and Paul K. Romano,
"Incorporation of a Multigroup Transport Capability in the OpenMC Monte Carlo
Particle Transport Code," *Trans. Am. Nucl. Soc.*, **117**, 679-681 (2017).

1
docs/source/pythonapi/base.rst
View file

@ -109,6 +109,7 @@ Constructing Tallies
openmc.CellbornFilter
openmc.SurfaceFilter
openmc.MeshFilter
openmc.MeshSurfaceFilter
openmc.EnergyFilter
openmc.EnergyoutFilter
openmc.MuFilter

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

@ -32,9 +32,12 @@ Core Functions
:template: myfunction.rst
openmc.data.atomic_mass
openmc.data.gnd_name
openmc.data.linearize
openmc.data.thin
openmc.data.water_density
openmc.data.write_compact_458_library
openmc.data.zam
Angle-Energy Distributions
--------------------------

85
docs/source/pythonapi/deplete.rst Normal file
View file

@ -0,0 +1,85 @@
.. _pythonapi_deplete:
----------------------------------
:mod:`openmc.deplete` -- Depletion
----------------------------------
.. module:: openmc.deplete
Two functions are provided that implement different time-integration algorithms
for depletion calculations.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
integrator.predictor
integrator.cecm
Each of these functions expects a "transport operator" to be passed. An operator
specific to OpenMC is available using the following class:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
Operator
When running in parallel using `mpi4py <http://mpi4py.scipy.org>`_, the MPI
intercommunicator used can be changed by modifying the following module
variable. If it is not explicitly modified, it defaults to
``mpi4py.MPI.COMM_WORLD``.
.. data:: comm
MPI intercommunicator used to call OpenMC library
:type: mpi4py.MPI.Comm
Internal Classes and Functions
------------------------------
During a depletion calculation, the depletion chain, reaction rates, and number
densities are managed through a series of internal classes that are not normally
visible to a user. However, should you find yourself wondering about these
classes (e.g., if you want to know what decay modes or reactions are present in
a depletion chain), they are documented here. The following classes store data
for a depletion chain:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
Chain
DecayTuple
Nuclide
ReactionTuple
The following classes are used during a depletion simulation and store auxiliary
data, such as number densities and reaction rates for each material.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
AtomNumber
OperatorResult
ReactionRates
Results
ResultsList
TransportOperator
Each of the integrator functions also relies on a number of "helper" functions
as follows:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
integrator.CRAM16
integrator.CRAM48

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

@ -15,14 +15,15 @@ there are many substantial benefits to using the Python API, including:
- The ability to define dimensions using variables.
- Availability of standard-library modules for working with files.
- An entire ecosystem of third-party packages for scientific computing.
- Ability to create materials based on natural elements or uranium enrichment
- Automated multi-group cross section generation (:mod:`openmc.mgxs`)
- A fully-featured nuclear data interface (:mod:`openmc.data`)
- Depletion capability (:mod:`openmc.deplete`)
- Convenience functions (e.g., a function returning a hexagonal region)
- Ability to plot individual universes as geometry is being created
- A :math:`k_\text{eff}` search function (:func:`openmc.search_for_keff`)
- Random sphere packing for generating TRISO particle locations
(:func:`openmc.model.pack_trisos`)
- A fully-featured nuclear data interface (:mod:`openmc.data`)
- Ability to create materials based on natural elements or uranium enrichment
For those new to Python, there are many good tutorials available online. We
recommend going through the modules from `Codecademy
@ -45,6 +46,7 @@ Modules
base
model
examples
deplete
mgxs
stats
data

6
docs/source/pythonapi/model.rst
View file

@ -5,16 +5,12 @@
Convenience Functions
---------------------
Several helper functions are available here. Ther first two create rectangular
and hexagonal prisms defined by the intersection of four and six surface
half-spaces, respectively. The last function takes a sequence of surfaces and
returns the regions that separate them.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.model.borated_water
openmc.model.get_hexagonal_prism
openmc.model.get_rectangular_prism
openmc.model.subdivide

24
docs/source/quickinstall.rst
View file

@ -57,7 +57,6 @@ are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
.. _HDF5: http://www.hdfgroup.org/HDF5/
---------------------------------------
Installing from Source on Ubuntu 15.04+
@ -83,9 +82,12 @@ building and installing OpenMC from source.
Installing from Source on Linux or Mac OS X
-------------------------------------------
All OpenMC source code is hosted on GitHub_. If you have git_, the gfortran_
compiler, CMake_, and HDF5_ installed, you can download and install OpenMC be
entering the following commands in a terminal:
All OpenMC source code is hosted on `GitHub
<https://github.com/mit-crpg/openmc>`_. If you have `git
<https://git-scm.com>`_, the `gcc <https://gcc.gnu.org/>`_ compiler suite,
`CMake <http://www.cmake.org>`_, and `HDF5 <https://www.hdfgroup.org/HDF5/>`_
installed, you can download and install OpenMC be entering the following
commands in a terminal:
.. code-block:: sh
@ -104,11 +106,15 @@ should specify an installation directory where you have write access, e.g.
cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local ..
The :mod:`openmc` Python package must be installed separately. The easiest way
to install it is using `pip <https://pip.pypa.io/en/stable/>`_, which is
included by default in Python 2.7 and Python 3.4+. From the root directory of
the OpenMC distribution/repository, run:
.. code-block:: sh
pip install .
If you want to build a parallel version of OpenMC (using OpenMP or MPI),
directions can be found in the :ref:`detailed installation instructions
<usersguide_build>`.
.. _GitHub: https://github.com/mit-crpg/openmc
.. _git: http://git-scm.com
.. _gfortran: http://gcc.gnu.org/wiki/GFortran
.. _CMake: http://www.cmake.org

4
docs/source/usersguide/beginners.rst
View file

@ -151,8 +151,8 @@ and `Volume II`_. You may also find it helpful to review the following terms:
.. _git: http://git-scm.com/
.. _git tutorials: http://git-scm.com/documentation
.. _Reactor Concepts Manual: http://www.tayloredge.com/periodic/trivia/ReactorConcepts.pdf
.. _Volume I: http://energy.gov/sites/prod/files/2013/06/f2/h1019v1.pdf
.. _Volume II: http://energy.gov/sites/prod/files/2013/06/f2/h1019v2.pdf
.. _Volume I: https://www.standards.doe.gov/standards-documents/1000/1019-bhdbk-1993-v1
.. _Volume II: https://www.standards.doe.gov/standards-documents/1000/1019-bhdbk-1993-v2
.. _OpenMC source code: https://github.com/mit-crpg/openmc
.. _GitHub: https://github.com/
.. _bug reports: https://github.com/mit-crpg/openmc/issues

106
docs/source/usersguide/install.rst
View file

@ -6,17 +6,18 @@ Installation and Configuration
.. currentmodule:: openmc
.. _install_conda:
----------------------------------------
Installing on Linux/Mac with conda-forge
----------------------------------------
`Conda <http://conda.pydata.org/docs/>`_ is an open source package management
system and environment management system for installing multiple versions of
software packages and their dependencies and switching easily between
them. `conda-forge <https://conda-forge.github.io/>`_ is a community-led conda
channel of installable packages. For instructions on installing conda, please
consult their `documentation
<http://conda.pydata.org/docs/install/quick.html>`_.
Conda_ is an open source package management system and environment management
system for installing multiple versions of software packages and their
dependencies and switching easily between them. `conda-forge
<https://conda-forge.github.io/>`_ is a community-led conda channel of
installable packages. For instructions on installing conda, please consult their
`documentation <http://conda.pydata.org/docs/install/quick.html>`_.
Once you have `conda` installed on your system, add the `conda-forge` channel to
your configuration with:
@ -38,6 +39,8 @@ It is possible to list all of the versions of OpenMC available on your platform
conda search openmc --channel conda-forge
.. _install_ppa:
-----------------------------
Installing on Ubuntu with PPA
-----------------------------
@ -68,9 +71,11 @@ are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
--------------------
Building from Source
--------------------
.. _install_source:
----------------------
Installing from Source
----------------------
.. _prerequisites:
@ -191,8 +196,8 @@ switch to the source of the latest stable release, run the following commands::
git checkout master
.. _GitHub: https://github.com/mit-crpg/openmc
.. _git: http://git-scm.com
.. _ssh: http://en.wikipedia.org/wiki/Secure_Shell
.. _git: https://git-scm.com
.. _ssh: https://en.wikipedia.org/wiki/Secure_Shell
.. _usersguide_build:
@ -258,14 +263,15 @@ should be used:
Compiling with MPI
++++++++++++++++++
To compile with MPI, set the :envvar:`FC` and :envvar:`CC` environment variables
to the path to the MPI Fortran and C wrappers, respectively. For example, in a
bash shell:
To compile with MPI, set the :envvar:`FC`, :envvar:`CC`, and :envvar:`CXX`
environment variables to the path to the MPI Fortran, C, and C++ wrappers,
respectively. For example, in a bash shell:
.. code-block:: sh
export FC=mpif90
export FC=mpifort
export CC=mpicc
export CXX=mpicxx
cmake /path/to/openmc
Note that in many shells, environment variables can be set for a single command,
@ -273,7 +279,7 @@ i.e.
.. code-block:: sh
FC=mpif90 CC=mpicc cmake /path/to/openmc
FC=mpifort CC=mpicc CXX=mpicxx cmake /path/to/openmc
Selecting HDF5 Installation
+++++++++++++++++++++++++++
@ -349,7 +355,7 @@ follows:
.. code-block:: sh
mkdir build && cd build
FC=ifort CC=icc FFLAGS=-mmic cmake -Dopenmp=on ..
FC=ifort CC=icc CXX=icpc FFLAGS=-mmic cmake -Dopenmp=on ..
make
Note that unless an HDF5 build for the Intel Xeon Phi (Knights Corner) is
@ -364,26 +370,57 @@ Testing Build
To run the test suite, you will first need to download a pre-generated cross
section library along with windowed multipole data. Please refer to our
:ref:`test suite` documentation for further details.
:ref:`devguide_tests` documentation for further details.
--------------------
Python Prerequisites
--------------------
---------------------
Installing Python API
---------------------
OpenMC's :ref:`Python API <pythonapi>` works with either Python 2.7 or Python
3.2+. In addition to Python itself, the API relies on a number of third-party
packages. All prerequisites can be installed using `conda
<http://conda.pydata.org/docs/>`_ (recommended), `pip
<https://pip.pypa.io/en/stable/>`_, or through the package manager in most Linux
If you installed OpenMC using :ref:`Conda <install_conda>` or :ref:`PPA
<install_ppa>`, no further steps are necessary in order to use OpenMC's
:ref:`Python API <pythonapi>`. However, if you are :ref:`installing from source
<install_source>`, the Python API is not installed by default when ``make
install`` is run because in many situations it doesn't make sense to install a
Python package in the same location as the ``openmc`` executable (for example,
if you are installing the package into a `virtual environment
<https://docs.python.org/3/tutorial/venv.html>`_). The easiest way to install
the :mod:`openmc` Python package is to use pip_, which is included by default in
Python 3.4+. From the root directory of the OpenMC distribution/repository, run:
.. code-block:: sh
pip install .
pip will first check that all :ref:`required third-party packages
<usersguide_python_prereqs>` have been installed, and if they are not present,
they will be installed by downloading the appropriate packages from the Python
Package Index (`PyPI <https://pypi.org/>`_). However, do note that since pip
runs the ``setup.py`` script which requires NumPy, you will have to first
install NumPy:
.. code-block:: sh
pip install numpy
Installing in "Development" Mode
--------------------------------
If you are primarily doing development with OpenMC, it is strongly recommended
to install the Python package in :ref:`"editable" mode <devguide_editable>`.
.. _usersguide_python_prereqs:
Prerequisites
-------------
The Python API works with Python 3.4+. In addition to Python itself, the API
relies on a number of third-party packages. All prerequisites can be installed
using Conda_ (recommended), pip_, or through the package manager in most Linux
distributions.
.. admonition:: Required
:class: error
`six <https://pythonhosted.org/six/>`_
The Python API works with both Python 2.7+ and 3.2+. To do so, the six
compatibility library is used.
`NumPy <http://www.numpy.org/>`_
NumPy is used extensively within the Python API for its powerful
N-dimensional array.
@ -415,6 +452,11 @@ distributions.
.. admonition:: Optional
:class: note
`mpi4py <http://mpi4py.scipy.org/>`_
mpi4py provides Python bindings to MPI for running distributed-memory
parallel runs. This package is needed if you plan on running depletion
simulations in parallel using MPI.
`Cython <http://cython.org/>`_
Cython is used for resonance reconstruction for ENDF data converted to
:class:`openmc.data.IncidentNeutron`.
@ -457,3 +499,5 @@ schemas.xml file in your own OpenMC source directory.
.. _RELAX NG: http://relaxng.org/
.. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html
.. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html
.. _Conda: https://conda.io/docs/
.. _pip: https://pip.pypa.io/en/stable/

14
docs/source/usersguide/materials.rst
View file

@ -43,14 +43,22 @@ of an element, you specify the element itself. For example,
Internally, OpenMC stores data on the atomic masses and natural abundances of
all known isotopes and then uses this data to determine what isotopes should be
added to the material. When the material is later exported to XML for use by the
:ref:`scripts_openmc` executable, you'll see that any natural elements are
:ref:`scripts_openmc` executable, you'll see that any natural elements were
expanded to the naturally-occurring isotopes.
The :meth:`Material.add_element` method can also be used to add uranium at a
specified enrichment through the `enrichment` argument. For example, the
following would add 3.2% enriched uranium to a material::
mat.add_element('U', 1.0, enrichment=3.2)
In addition to U235 and U238, concentrations of U234 and U236 will be present
and are determined through a correlation based on measured data.
Often, cross section libraries don't actually have all naturally-occurring
isotopes for a given element. For example, in ENDF/B-VII.1, cross section
evaluations are given for O16 and O17 but not for O18. If OpenMC is aware of
what cross sections you will be using (either through the
:attr:`Materials.cross_sections` attribute or the
what cross sections you will be using (through the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable), it will attempt to only
put isotopes in your model for which you have cross section data. In the case of
oxygen in ENDF/B-VII.1, the abundance of O18 would end up being lumped with O16.

4
docs/source/usersguide/scripts.rst
View file

@ -168,8 +168,8 @@ ENDF/B-VII.1. It has the following optional arguments:
This script downloads `ENDF/B-VII.1 ACE data
<http://www.nndc.bnl.gov/endf/b7.1/acefiles.html>`_ from NNDC and converts it to
an HDF5 library for use with OpenMC. This data is used for OpenMC's regression
test suite. This script has the following optional arguments:
an HDF5 library for use with OpenMC. This script has the following optional
arguments:
-b, --batch Suppress standard in

2
examples/xml/pincell/tallies.xml
View file

@ -8,7 +8,7 @@
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
<bins>2</bins>
</filter>
<filter id="2" type="energy">

1
include/openmc.h
View file

@ -52,6 +52,7 @@ extern "C" {
int openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n);
int openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins);
int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_next_batch();
int openmc_nuclide_name(int index, char** name);
void openmc_plot_geometry();

49
openmc/_utils.py Normal file
View file

@ -0,0 +1,49 @@
import os.path
from pathlib import Path
from urllib.parse import urlparse
from urllib.request import urlopen
_BLOCK_SIZE = 16384
def download(url):
"""Download file from a URL
Parameters
----------
url : str
URL from which to download
Returns
-------
basename : str
Name of file written locally
"""
req = urlopen(url)
# Get file size from header
file_size = req.length
# Check if file already downloaded
basename = Path(urlparse(url).path).name
if os.path.exists(basename):
if os.path.getsize(basename) == file_size:
print('Skipping {}, already downloaded'.format(basename))
return basename
# Copy file to disk in chunks
print('Downloading {}... '.format(basename), end='')
downloaded = 0
with open(basename, 'wb') as fh:
while True:
chunk = req.read(_BLOCK_SIZE)
if not chunk:
break
fh.write(chunk)
downloaded += len(chunk)
status = '{:10} [{:3.2f}%]'.format(
downloaded, downloaded * 100. / file_size)
print(status + '\b'*len(status), end='')
print('')
return basename

22
openmc/arithmetic.py
View file

@ -2,7 +2,6 @@ import sys
import copy
from collections import Iterable
from six import string_types
import numpy as np
import pandas as pd
@ -86,18 +85,18 @@ class CrossScore(object):
@left_score.setter
def left_score(self, left_score):
cv.check_type('left_score', left_score,
string_types + (CrossScore, AggregateScore))
(str, CrossScore, AggregateScore))
self._left_score = left_score
@right_score.setter
def right_score(self, right_score):
cv.check_type('right_score', right_score,
string_types + (CrossScore, AggregateScore))
(str, CrossScore, AggregateScore))
self._right_score = right_score
@binary_op.setter
def binary_op(self, binary_op):
cv.check_type('binary_op', binary_op, string_types)
cv.check_type('binary_op', binary_op, str)
cv.check_value('binary_op', binary_op, _TALLY_ARITHMETIC_OPS)
self._binary_op = binary_op
@ -202,7 +201,7 @@ class CrossNuclide(object):
@binary_op.setter
def binary_op(self, binary_op):
cv.check_type('binary_op', binary_op, string_types)
cv.check_type('binary_op', binary_op, str)
cv.check_value('binary_op', binary_op, _TALLY_ARITHMETIC_OPS)
self._binary_op = binary_op
@ -335,7 +334,7 @@ class CrossFilter(object):
@binary_op.setter
def binary_op(self, binary_op):
cv.check_type('binary_op', binary_op, string_types)
cv.check_type('binary_op', binary_op, str)
cv.check_value('binary_op', binary_op, _TALLY_ARITHMETIC_OPS)
self._binary_op = binary_op
@ -482,12 +481,12 @@ class AggregateScore(object):
@scores.setter
def scores(self, scores):
cv.check_iterable_type('scores', scores, string_types)
cv.check_iterable_type('scores', scores, str)
self._scores = scores
@aggregate_op.setter
def aggregate_op(self, aggregate_op):
cv.check_type('aggregate_op', aggregate_op, string_types +(CrossScore,))
cv.check_type('aggregate_op', aggregate_op, (str, CrossScore))
cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
self._aggregate_op = aggregate_op
@ -561,13 +560,12 @@ class AggregateNuclide(object):
@nuclides.setter
def nuclides(self, nuclides):
cv.check_iterable_type('nuclides', nuclides,
string_types + (openmc.Nuclide, CrossNuclide))
cv.check_iterable_type('nuclides', nuclides, (str, CrossNuclide))
self._nuclides = nuclides
@aggregate_op.setter
def aggregate_op(self, aggregate_op):
cv.check_type('aggregate_op', aggregate_op, string_types)
cv.check_type('aggregate_op', aggregate_op, str)
cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
self._aggregate_op = aggregate_op
@ -690,7 +688,7 @@ class AggregateFilter(object):
@aggregate_op.setter
def aggregate_op(self, aggregate_op):
cv.check_type('aggregate_op', aggregate_op, string_types)
cv.check_type('aggregate_op', aggregate_op, str)
cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
self._aggregate_op = aggregate_op

9
openmc/capi/__init__.py
View file

@ -15,7 +15,6 @@ objects in the :mod:`openmc.capi` subpackage, for example:
from ctypes import CDLL
import os
import sys
from warnings import warn
import pkg_resources
@ -36,10 +35,7 @@ else:
# available. Instead, we create a mock object so that when the modules
# within the openmc.capi package try to configure arguments and return
# values for symbols, no errors occur
try:
from unittest.mock import Mock
except ImportError:
from mock import Mock
from unittest.mock import Mock
_dll = Mock()
from .error import *
@ -50,6 +46,3 @@ from .cell import *
from .filter import *
from .tally import *
from .settings import settings
warn("The Python bindings to OpenMC's C API are still unstable "
"and may change substantially in future releases.", FutureWarning)

9
openmc/capi/cell.py
View file

@ -1,4 +1,4 @@
from collections import Mapping, Iterable
from collections.abc import Mapping, Iterable
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
@ -65,10 +65,7 @@ class Cell(_FortranObjectWithID):
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A cell with ID={} has already '
@ -81,7 +78,7 @@ class Cell(_FortranObjectWithID):
index = mapping[uid]._index
if index not in cls.__instances:
instance = super(Cell, cls).__new__(cls)
instance = super().__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid

23
openmc/capi/filter.py
View file

@ -1,4 +1,4 @@
from collections import Mapping
from collections.abc import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, \
create_string_buffer
from weakref import WeakValueDictionary
@ -55,6 +55,9 @@ _dll.openmc_material_filter_set_bins.errcheck = _error_handler
_dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32]
_dll.openmc_mesh_filter_set_mesh.restype = c_int
_dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
class Filter(_FortranObjectWithID):
@ -66,10 +69,7 @@ class Filter(_FortranObjectWithID):
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A filter with ID={} has already '
@ -87,7 +87,7 @@ class Filter(_FortranObjectWithID):
index = mapping[uid]._index
if index not in cls.__instances:
instance = super(Filter, cls).__new__(cls)
instance = super().__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
@ -110,7 +110,7 @@ class EnergyFilter(Filter):
filter_type = 'energy'
def __init__(self, bins=None, uid=None, new=True, index=None):
super(EnergyFilter, self).__init__(uid, new, index)
super().__init__(uid, new, index)
if bins is not None:
self.bins = bins
@ -131,7 +131,7 @@ class EnergyFilter(Filter):
self._index, len(energies), energies_p)
class EnergyoutFilter(Filter):
class EnergyoutFilter(EnergyFilter):
filter_type = 'energyout'
@ -167,7 +167,7 @@ class MaterialFilter(Filter):
filter_type = 'material'
def __init__(self, bins=None, uid=None, new=True, index=None):
super(MaterialFilter, self).__init__(uid, new, index)
super().__init__(uid, new, index)
if bins is not None:
self.bins = bins
@ -191,6 +191,10 @@ class MeshFilter(Filter):
filter_type = 'mesh'
class MeshSurfaceFilter(Filter):
filter_type = 'meshsurface'
class MuFilter(Filter):
filter_type = 'mu'
@ -219,6 +223,7 @@ _FILTER_TYPE_MAP = {
'energyfunction': EnergyFunctionFilter,
'material': MaterialFilter,
'mesh': MeshFilter,
'meshsurface': MeshSurfaceFilter,
'mu': MuFilter,
'polar': PolarFilter,
'surface': SurfaceFilter,

7
openmc/capi/material.py
View file

@ -1,4 +1,4 @@
from collections import Mapping
from collections.abc import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
@ -78,10 +78,7 @@ class Material(_FortranObjectWithID):
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A material with ID={} has already '

4
openmc/capi/nuclide.py
View file

@ -1,4 +1,4 @@
from collections import Mapping
from collections.abc import Mapping
from ctypes import c_int, c_char_p, POINTER
from weakref import WeakValueDictionary
@ -58,7 +58,7 @@ class Nuclide(_FortranObject):
def __new__(cls, *args):
if args not in cls.__instances:
instance = super(Nuclide, cls).__new__(cls)
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]

5
openmc/capi/settings.py
View file

@ -11,8 +11,7 @@ _RUN_MODES = {1: 'fixed source',
5: 'volume'}
_dll.openmc_set_seed.argtypes = [c_int64]
_dll.openmc_set_seed.restype = c_int
_dll.openmc_set_seed.errcheck = _error_handler
_dll.openmc_get_seed.restype = c_int64
class _Settings(object):
@ -43,7 +42,7 @@ class _Settings(object):
@property
def seed(self):
return c_int64.in_dll(_dll, 'seed').value
return _dll.openmc_get_seed()
@seed.setter
def seed(self, seed):

9
openmc/capi/tally.py
View file

@ -1,4 +1,4 @@
from collections import Mapping
from collections.abc import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
@ -155,10 +155,7 @@ class Tally(_FortranObjectWithID):
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A tally with ID={} has already '
@ -172,7 +169,7 @@ class Tally(_FortranObjectWithID):
index = mapping[uid]._index
if index not in cls.__instances:
instance = super(Tally, cls).__new__(cls)
instance = super().__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid

63
openmc/cell.py
View file

@ -1,4 +1,5 @@
from collections import OrderedDict, Iterable
from collections import OrderedDict
from collections.abc import Iterable
from copy import deepcopy
from math import cos, sin, pi
from numbers import Real, Integral
@ -6,7 +7,6 @@ from xml.etree import ElementTree as ET
import sys
import warnings
from six import string_types
import numpy as np
import openmc
@ -203,7 +203,7 @@ class Cell(IDManagerMixin):
@name.setter
def name(self, name):
if name is not None:
cv.check_type('cell name', name, string_types)
cv.check_type('cell name', name, str)
self._name = name
else:
self._name = ''
@ -211,14 +211,7 @@ class Cell(IDManagerMixin):
@fill.setter
def fill(self, fill):
if fill is not None:
if isinstance(fill, string_types):
if fill.strip().lower() != 'void':
msg = 'Unable to set Cell ID="{0}" to use a non-Material ' \
'or Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
fill = None
elif isinstance(fill, Iterable):
if isinstance(fill, Iterable):
for i, f in enumerate(fill):
if f is not None:
cv.check_type('cell.fill[i]', f, openmc.Material)
@ -291,50 +284,6 @@ class Cell(IDManagerMixin):
cv.check_type('cell volume', volume, Real)
self._volume = volume
def add_surface(self, surface, halfspace):
"""Add a half-space to the list of half-spaces whose intersection defines the
cell.
.. deprecated:: 0.7.1
Use the :attr:`Cell.region` property to directly specify a Region
expression.
Parameters
----------
surface : openmc.Surface
Quadric surface dividing space
halfspace : {-1, 1}
Indicate whether the negative or positive half-space is to be used
"""
warnings.warn("Cell.add_surface(...) has been deprecated and may be "
"removed in a future version. The region for a Cell "
"should be defined using the region property directly.",
DeprecationWarning)
if not isinstance(surface, openmc.Surface):
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
if halfspace not in [-1, +1]:
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" with halfspace ' \
'"{2}" since it is not +/-1'.format(surface, self._id, halfspace)
raise ValueError(msg)
# If no region has been assigned, simply use the half-space. Otherwise,
# take the intersection of the current region and the half-space
# specified
region = +surface if halfspace == 1 else -surface
if self.region is None:
self.region = region
else:
if isinstance(self.region, Intersection):
self.region &= region
else:
self.region = Intersection(self.region, region)
def add_volume_information(self, volume_calc):
"""Add volume information to a cell.
@ -346,7 +295,7 @@ class Cell(IDManagerMixin):
"""
if volume_calc.domain_type == 'cell':
if self.id in volume_calc.volumes:
self._volume = volume_calc.volumes[self.id][0]
self._volume = volume_calc.volumes[self.id].n
self._atoms = volume_calc.atoms[self.id]
else:
raise ValueError('No volume information found for this cell.')
@ -386,7 +335,7 @@ class Cell(IDManagerMixin):
volume = self.volume
for name, atoms in self._atoms.items():
nuclide = openmc.Nuclide(name)
density = 1.0e-24 * atoms[0]/volume # density in atoms/b-cm
density = 1.0e-24 * atoms.n/volume # density in atoms/b-cm
nuclides[name] = (nuclide, density)
else:
raise RuntimeError(

12
openmc/checkvalue.py
View file

@ -1,5 +1,5 @@
import copy
from collections import Iterable
from collections.abc import Iterable
import numpy as np
@ -246,9 +246,9 @@ def check_filetype_version(obj, expected_type, expected_version):
----------
obj : h5py.File
HDF5 file to check
expected_type
expected_type : str
Expected file type, e.g. 'statepoint'
expected_version
expected_version : int
Expected major version number.
"""
@ -288,7 +288,7 @@ class CheckedList(list):
"""
def __init__(self, expected_type, name, items=[]):
super(CheckedList, self).__init__()
super().__init__()
self.expected_type = expected_type
self.name = name
for item in items:
@ -319,7 +319,7 @@ class CheckedList(list):
"""
check_type(self.name, item, self.expected_type)
super(CheckedList, self).append(item)
super().append(item)
def insert(self, index, item):
"""Insert item before index
@ -333,4 +333,4 @@ class CheckedList(list):
"""
check_type(self.name, item, self.expected_type)
super(CheckedList, self).insert(index, item)
super().insert(index, item)

6
openmc/cmfd.py
View file

@ -10,13 +10,11 @@ References
"""
from collections import Iterable
from collections.abc import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
from six import string_types
from openmc.clean_xml import clean_xml_indentation
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than, check_less_than)
@ -338,7 +336,7 @@ class CMFD(object):
@display.setter
def display(self, display):
check_type('CMFD display', display, string_types)
check_type('CMFD display', display, str)
check_value('CMFD display', display,
['balance', 'dominance', 'entropy', 'source'])
self._display = display

8
openmc/data/ace.py
View file

@ -15,16 +15,14 @@ generates ACE-format cross sections.
"""
from __future__ import division, unicode_literals
from os import SEEK_CUR
import struct
import sys
from six import string_types
import numpy as np
from openmc.mixin import EqualityMixin
from openmc.data.endf import ENDF_FLOAT_RE
from openmc.data.endf import _ENDF_FLOAT_RE
def ascii_to_binary(ascii_file, binary_file):
"""Convert an ACE file in ASCII format (type 1) to binary format (type 2).
@ -153,7 +151,7 @@ class Library(EqualityMixin):
"""
def __init__(self, filename, table_names=None, verbose=False):
if isinstance(table_names, string_types):
if isinstance(table_names, str):
table_names = [table_names]
if table_names is not None:
table_names = set(table_names)
@ -351,7 +349,7 @@ class Library(EqualityMixin):
# after it). If it's too short, then we apply the ENDF float regular
# expression. We don't do this by default because it's expensive!
if xss.size != nxs[1] + 1:
datastr = ENDF_FLOAT_RE.sub(r'\1e\2', datastr)
datastr = _ENDF_FLOAT_RE.sub(r'\1e\2', datastr)
xss = np.fromstring(datastr, sep=' ')
assert xss.size == nxs[1] + 1

4
openmc/data/angle_distribution.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from io import StringIO
from numbers import Real
from warnings import warn
@ -34,7 +34,7 @@ class AngleDistribution(EqualityMixin):
"""
def __init__(self, energy, mu):
super(AngleDistribution, self).__init__()
super().__init__()
self.energy = energy
self.mu = mu

5
openmc/data/angle_energy.py
View file

@ -1,14 +1,11 @@
from abc import ABCMeta, abstractmethod
from io import StringIO
from six import add_metaclass
import openmc.data
from openmc.mixin import EqualityMixin
@add_metaclass(ABCMeta)
class AngleEnergy(EqualityMixin):
class AngleEnergy(EqualityMixin, metaclass=ABCMeta):
"""Distribution in angle and energy of a secondary particle."""
@abstractmethod
def to_hdf5(self, group):

4
openmc/data/correlated.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from numbers import Real, Integral
from warnings import warn
@ -45,7 +45,7 @@ class CorrelatedAngleEnergy(AngleEnergy):
"""
def __init__(self, breakpoints, interpolation, energy, energy_out, mu):
super(CorrelatedAngleEnergy, self).__init__()
super().__init__()
self.breakpoints = breakpoints
self.interpolation = interpolation
self.energy = energy

180
openmc/data/data.py
View file

@ -1,6 +1,8 @@
import itertools
from math import sqrt
import os
import re
from warnings import warn
# Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
@ -133,23 +135,24 @@ ATOMIC_NUMBER = {value: key for key, value in ATOMIC_SYMBOL.items()}
_ATOMIC_MASS = {}
_GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
def atomic_mass(isotope):
"""Return atomic mass of isotope in atomic mass units.
Atomic mass data comes from the Atomic Mass Evaluation 2012, published in
Chinese Physics C 36 (2012), 1287--1602.
Atomic mass data comes from the `Atomic Mass Evaluation 2012
<https://www-nds.iaea.org/amdc/ame2012/AME2012-1.pdf>`_.
Parameters
----------
isotope : str
Name of isotope, e.g. 'Pu239'
Name of isotope, e.g., 'Pu239'
Returns
-------
float or None
Atomic mass of isotope in atomic mass units. If the isotope listed does
not have a known atomic mass, None is returned.
float
Atomic mass of isotope in [amu]
"""
if not _ATOMIC_MASS:
@ -180,7 +183,7 @@ def atomic_mass(isotope):
if '_' in isotope:
isotope = isotope[:isotope.find('_')]
return _ATOMIC_MASS.get(isotope.lower())
return _ATOMIC_MASS[isotope.lower()]
def atomic_weight(element):
@ -196,16 +199,168 @@ def atomic_weight(element):
Returns
-------
float or None
Atomic weight of element in atomic mass units. If the element listed does
not exist, None is returned.
float
Atomic weight of element in [amu]
"""
weight = 0.
for nuclide, abundance in NATURAL_ABUNDANCE.items():
if re.match(r'{}\d+'.format(element), nuclide):
weight += atomic_mass(nuclide) * abundance
return None if weight == 0. else weight
if weight > 0.:
return weight
else:
raise ValueError("No naturally-occurring isotopes for element '{}'."
.format(element))
def water_density(temperature, pressure=0.1013):
"""Return the density of liquid water at a given temperature and pressure.
The density is calculated from a polynomial fit using equations and values
from the 2012 version of the IAPWS-IF97 formulation. Only the equations
for region 1 are implemented here. Region 1 is limited to liquid water
below 100 [MPa] with a temperature above 273.15 [K], below 623.15 [K], and
below saturation.
Reference: International Association for the Properties of Water and Steam,
"Revised Release on the IAPWS Industrial Formulation 1997 for the
Thermodynamic Properties of Water and Steam", IAPWS R7-97(2012).
Parameters
----------
temperature : float
Water temperature in units of [K]
pressure : float
Water pressure in units of [MPa]
Returns
-------
float
Water density in units of [g/cm^3]
"""
# Make sure the temperature and pressure are inside the min/max region 1
# bounds. (Relax the 273.15 bound to 273 in case a user wants 0 deg C data
# but they only use 3 digits for their conversion to K.)
if pressure > 100.0:
warn("Results are not valid for pressures above 100 MPa.")
if pressure < 0.0:
warn("Results are not valid for pressures below zero.")
if temperature < 273:
warn("Results are not valid for temperatures below 273.15 K.")
if temperature > 623.15:
warn("Results are not valid for temperatures above 623.15 K.")
# IAPWS region 4 parameters
n4 = [0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
-0.48232657361591e4, 0.40511340542057e6, -0.23855557567849,
0.65017534844798e3]
# Compute the saturation temperature at the given pressure.
beta = pressure**(0.25)
E = beta**2 + n4[2] * beta + n4[5]
F = n4[0] * beta**2 + n4[3] * beta + n4[6]
G = n4[1] * beta**2 + n4[4] * beta + n4[7]
D = 2.0 * G / (-F - sqrt(F**2 - 4 * E * G))
T_sat = 0.5 * (n4[9] + D
- sqrt((n4[9] + D)**2 - 4.0 * (n4[8] + n4[9] * D)))
# Make sure we aren't above saturation. (Relax this bound by .2 degrees
# for deg C to K conversions.)
if temperature > T_sat + 0.2:
warn("Results are not valid for temperatures above saturation "
"(above the boiling point).")
# IAPWS region 1 parameters
R_GAS_CONSTANT = 0.461526 # kJ / kg / K
ref_p = 16.53 # MPa
ref_T = 1386 # K
n1f = [0.14632971213167, -0.84548187169114, -0.37563603672040e1,
0.33855169168385e1, -0.95791963387872, 0.15772038513228,
-0.16616417199501e-1, 0.81214629983568e-3, 0.28319080123804e-3,
-0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1,
-0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3,
-0.30001780793026e-3, 0.47661393906987e-4, -0.44141845330846e-5,
-0.72694996297594e-15, -0.31679644845054e-4, -0.28270797985312e-5,
-0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6,
-0.14341729937924e-12, -0.40516996860117e-6, -0.12734301741641e-8,
-0.17424871230634e-9, -0.68762131295531e-18, 0.14478307828521e-19,
0.26335781662795e-22, -0.11947622640071e-22, 0.18228094581404e-23,
-0.93537087292458e-25]
I1f = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4,
4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
J1f = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4,
0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
# Nondimensionalize the pressure and temperature.
pi = pressure / ref_p
tau = ref_T / temperature
# Compute the derivative of gamma (dimensionless Gibbs free energy) with
# respect to pi.
gamma1_pi = 0.0
for n, I, J in zip(n1f, I1f, J1f):
gamma1_pi -= n * I * (7.1 - pi)**(I - 1) * (tau - 1.222)**J
# Compute the leading coefficient. This sets the units at
# 1 [MPa] * [kg K / kJ] * [1 / K]
# = 1e6 [N / m^2] * 1e-3 [kg K / N / m] * [1 / K]
# = 1e3 [kg / m^3]
# = 1 [g / cm^3]
coeff = pressure / R_GAS_CONSTANT / temperature
# Compute and return the density.
return coeff / pi / gamma1_pi
def gnd_name(Z, A, m=0):
"""Return nuclide name using GND convention
Parameters
----------
Z : int
Atomic number
A : int
Mass number
m : int, optional
Metastable state
Returns
-------
str
Nuclide name in GND convention, e.g., 'Am242_m1'
"""
if m > 0:
return '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, m)
else:
return '{}{}'.format(ATOMIC_SYMBOL[Z], A)
def zam(name):
"""Return tuple of (atomic number, mass number, metastable state)
Parameters
----------
name : str
Name of nuclide using GND convention, e.g., 'Am242_m1'
Returns
-------
3-tuple of int
Atomic number, mass number, and metastable state
"""
try:
symbol, A, state = _GND_NAME_RE.match(name).groups()
except AttributeError:
raise ValueError("'{}' does not appear to be a nuclide name in GND "
"format.".format(name))
metastable = int(state[2:]) if state else 0
return (ATOMIC_NUMBER[symbol], int(A), metastable)
# Values here are from the Committee on Data for Science and Technology
@ -214,8 +369,9 @@ def atomic_weight(element):
# The value of the Boltzman constant in units of eV / K
K_BOLTZMANN = 8.6173303e-5
# Used for converting units in ACE data
# Unit conversions
EV_PER_MEV = 1.0e6
JOULE_PER_EV = 1.6021766208e-19
# Avogadro's constant
AVOGADRO = 6.022140857e23

14
openmc/data/decay.py
View file

@ -1,16 +1,13 @@
from collections import Iterable, namedtuple
from collections import namedtuple
from collections.abc import Iterable
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'])
from uncertainties import ufloat, unumpy, UFloat
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
@ -278,12 +275,12 @@ class DecayMode(EqualityMixin):
@modes.setter
def modes(self, modes):
cv.check_type('decay modes', modes, Iterable, string_types)
cv.check_type('decay modes', modes, Iterable, str)
self._modes = modes
@parent.setter
def parent(self, parent):
cv.check_type('parent nuclide', parent, string_types)
cv.check_type('parent nuclide', parent, str)
self._parent = parent
@ -457,6 +454,7 @@ class Decay(EqualityMixin):
items, values = get_list_record(file_obj)
self.nuclide['spin'] = items[0]
self.nuclide['parity'] = items[1]
self.half_life = ufloat(float('inf'), float('inf'))
@property
def decay_constant(self):

28
openmc/data/endf.py
View file

@ -6,19 +6,18 @@ Data File ENDF-6". The latest version from June 2009 can be found at
http://www-nds.iaea.org/ndspub/documents/endf/endf102/endf102.pdf
"""
from __future__ import print_function, division, unicode_literals
import io
import re
import os
from math import pi
from collections import OrderedDict, Iterable
from pathlib import PurePath
from collections import OrderedDict
from collections.abc import Iterable
from six import string_types
import numpy as np
from numpy.polynomial.polynomial import Polynomial
from .data import ATOMIC_SYMBOL
from .data import ATOMIC_SYMBOL, gnd_name
from .function import Tabulated1D, INTERPOLATION_SCHEME
from openmc.stats.univariate import Uniform, Tabular, Legendre
@ -46,7 +45,7 @@ SUM_RULES = {1: [2, 3],
106: list(range(750, 800)),
107: list(range(800, 850))}
ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
_ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
def float_endf(s):
@ -69,7 +68,7 @@ def float_endf(s):
The number
"""
return float(ENDF_FLOAT_RE.sub(r'\1e\2', s))
return float(_ENDF_FLOAT_RE.sub(r'\1e\2', s))
def get_text_record(file_obj):
@ -250,6 +249,7 @@ 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.
@ -301,8 +301,8 @@ class Evaluation(object):
"""
def __init__(self, filename_or_obj):
if isinstance(filename_or_obj, string_types):
fh = open(filename_or_obj, 'r')
if isinstance(filename_or_obj, (str, PurePath)):
fh = open(str(filename_or_obj), 'r')
else:
fh = filename_or_obj
self.section = {}
@ -425,13 +425,9 @@ class Evaluation(object):
@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)
return gnd_name(self.target['atomic_number'],
self.target['mass_number'],
self.target['isomeric_state'])
class Tabulated2D(object):

24
openmc/data/energy_distribution.py
View file

@ -1,9 +1,8 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
from collections.abc import Iterable
from numbers import Integral, Real
from warnings import warn
from six import add_metaclass
import numpy as np
from .function import Tabulated1D, INTERPOLATION_SCHEME
@ -14,8 +13,7 @@ from .data import EV_PER_MEV
from .endf import get_tab1_record, get_tab2_record
@add_metaclass(ABCMeta)
class EnergyDistribution(EqualityMixin):
class EnergyDistribution(EqualityMixin, metaclass=ABCMeta):
"""Abstract superclass for all energy distributions."""
def __init__(self):
pass
@ -116,7 +114,7 @@ class ArbitraryTabulated(EnergyDistribution):
"""
def __init__(self, energy, pdf):
super(ArbitraryTabulated, self).__init__()
super().__init__()
self.energy = energy
self.pdf = pdf
@ -184,7 +182,7 @@ class GeneralEvaporation(EnergyDistribution):
"""
def __init__(self, theta, g, u):
super(GeneralEvaporation, self).__init__()
super().__init__()
self.theta = theta
self.g = g
self.u = u
@ -247,7 +245,7 @@ class MaxwellEnergy(EnergyDistribution):
"""
def __init__(self, theta, u):
super(MaxwellEnergy, self).__init__()
super().__init__()
self.theta = theta
self.u = u
@ -380,7 +378,7 @@ class Evaporation(EnergyDistribution):
"""
def __init__(self, theta, u):
super(Evaporation, self).__init__()
super().__init__()
self.theta = theta
self.u = u
@ -516,7 +514,7 @@ class WattEnergy(EnergyDistribution):
"""
def __init__(self, a, b, u):
super(WattEnergy, self).__init__()
super().__init__()
self.a = a
self.b = b
self.u = u
@ -684,7 +682,7 @@ class MadlandNix(EnergyDistribution):
"""
def __init__(self, efl, efh, tm):
super(MadlandNix, self).__init__()
super().__init__()
self.efl = efl
self.efh = efh
self.tm = tm
@ -807,7 +805,7 @@ class DiscretePhoton(EnergyDistribution):
"""
def __init__(self, primary_flag, energy, atomic_weight_ratio):
super(DiscretePhoton, self).__init__()
super().__init__()
self.primary_flag = primary_flag
self.energy = energy
self.atomic_weight_ratio = atomic_weight_ratio
@ -916,7 +914,7 @@ class LevelInelastic(EnergyDistribution):
"""
def __init__(self, threshold, mass_ratio):
super(LevelInelastic, self).__init__()
super().__init__()
self.threshold = threshold
self.mass_ratio = mass_ratio
@ -1021,7 +1019,7 @@ class ContinuousTabular(EnergyDistribution):
"""
def __init__(self, breakpoints, interpolation, energy, energy_out):
super(ContinuousTabular, self).__init__()
super().__init__()
self.breakpoints = breakpoints
self.interpolation = interpolation
self.energy = energy

2
openmc/data/fission_energy.py
View file

@ -1,4 +1,4 @@
from collections import Callable
from collections.abc import Callable
from copy import deepcopy
from io import StringIO
import sys

6
openmc/data/function.py
View file

@ -1,8 +1,7 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable, Callable
from collections.abc import Iterable, Callable
from numbers import Real, Integral
from six import add_metaclass
import numpy as np
import openmc.data
@ -14,8 +13,7 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log',
4: 'log-linear', 5: 'log-log'}
@add_metaclass(ABCMeta)
class Function1D(EqualityMixin):
class Function1D(EqualityMixin, metaclass=ABCMeta):
"""A function of one independent variable with HDF5 support."""
@abstractmethod
def __call__(self): pass

4
openmc/data/kalbach_mann.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from numbers import Real, Integral
from warnings import warn
@ -53,7 +53,7 @@ class KalbachMann(AngleEnergy):
def __init__(self, breakpoints, interpolation, energy, energy_out,
precompound, slope):
super(KalbachMann, self).__init__()
super().__init__()
self.breakpoints = breakpoints
self.interpolation = interpolation
self.energy = energy

4
openmc/data/laboratory.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from numbers import Real, Integral
import numpy as np
@ -44,7 +44,7 @@ class LaboratoryAngleEnergy(AngleEnergy):
"""
def __init__(self, breakpoints, interpolation, energy, mu, energy_out):
super(LaboratoryAngleEnergy, self).__init__()
super().__init__()
self.breakpoints = breakpoints
self.interpolation = interpolation
self.energy = energy

3
openmc/data/library.py
View file

@ -1,6 +1,5 @@
import os
import xml.etree.ElementTree as ET
from six import string_types
import h5py
@ -125,7 +124,7 @@ class DataLibrary(EqualityMixin):
raise ValueError("Either path or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
check_type('path', path, string_types)
check_type('path', path, str)
tree = ET.parse(path)
root = tree.getroot()

171
openmc/data/multipole.py
View file

@ -3,7 +3,6 @@ from math import exp, erf, pi, sqrt
import h5py
import numpy as np
from six import string_types
from . import WMP_VERSION
from .data import K_BOLTZMANN
@ -25,7 +24,7 @@ _RM_RF = 3 # Residue fission
# Multi-level Breit Wigner indices
_MLBW_RT = 1 # Residue total
_MLBW_RX = 2 # Residue compettitive
_MLBW_RX = 2 # Residue competitive
_MLBW_RA = 3 # Residue absorption
_MLBW_RF = 4 # Residue fission
@ -142,6 +141,12 @@ def _broaden_wmp_polynomials(E, dopp, n):
class WindowedMultipole(EqualityMixin):
"""Resonant cross sections represented in the windowed multipole format.
Parameters
----------
formalism : {'MLBW', 'RM'}
The R-matrix formalism used to reconstruct resonances. Either 'MLBW'
for multi-level Breit Wigner or 'RM' for Reich-Moore.
Attributes
----------
num_l : Integral
@ -196,11 +201,9 @@ class WindowedMultipole(EqualityMixin):
a/E + b/sqrt(E) + c + d sqrt(E) + ...
"""
def __init__(self):
self.num_l = None
self.fit_order = None
self.fissionable = None
self.formalism = None
def __init__(self, formalism):
self._num_l = None
self.formalism = formalism
self.spacing = None
self.sqrtAWR = None
self.start_E = None
@ -219,11 +222,15 @@ class WindowedMultipole(EqualityMixin):
@property
def fit_order(self):
return self._fit_order
return self.curvefit.shape[1] - 1
@property
def fissionable(self):
return self._fissionable
if self.formalism == 'RM':
return self.data.shape[1] == 4
else:
# Assume self.formalism == 'MLBW'
return self.data.shape[1] == 5
@property
def formalism(self):
@ -273,35 +280,10 @@ class WindowedMultipole(EqualityMixin):
def curvefit(self):
return self._curvefit
@num_l.setter
def num_l(self, num_l):
if num_l is not None:
cv.check_type('num_l', num_l, Integral)
cv.check_greater_than('num_l', num_l, 1, equality=True)
cv.check_less_than('num_l', num_l, 4, equality=True)
# There is an if block in _evaluate that assumes num_l <= 4.
self._num_l = num_l
@fit_order.setter
def fit_order(self, fit_order):
if fit_order is not None:
cv.check_type('fit_order', fit_order, Integral)
cv.check_greater_than('fit_order', fit_order, 2, equality=True)
# _broaden_wmp_polynomials assumes the curve fit has at least 3
# terms.
self._fit_order = fit_order
@fissionable.setter
def fissionable(self, fissionable):
if fissionable is not None:
cv.check_type('fissionable', fissionable, bool)
self._fissionable = fissionable
@formalism.setter
def formalism(self, formalism):
if formalism is not None:
cv.check_type('formalism', formalism, string_types)
cv.check_value('formalism', formalism, ('MLBW', 'RM'))
cv.check_type('formalism', formalism, str)
cv.check_value('formalism', formalism, ('MLBW', 'RM'))
self._formalism = formalism
@spacing.setter
@ -338,9 +320,20 @@ class WindowedMultipole(EqualityMixin):
cv.check_type('data', data, np.ndarray)
if len(data.shape) != 2:
raise ValueError('Multipole data arrays must be 2D')
if data.shape[1] not in (3, 4, 5): # 3 or 4 for RM, 4 or 5 for MLBW
raise ValueError('The second dimension of multipole data arrays'
' must have a length of 3, 4 or 5')
if self.formalism == 'RM':
if data.shape[1] not in (3, 4):
raise ValueError('For the Reich-Moore formalism, '
'data.shape[1] must be 3 or 4. One value for the pole.'
' One each for the total and absorption residues. '
'Possibly one more for a fission residue.')
else:
# Assume self.formalism == 'MLBW'
if data.shape[1] not in (4, 5):
raise ValueError('For the Multi-level Breit-Wigner '
'formalism, data.shape[1] must be 4 or 5. One value '
'for the pole. One each for the total, competitive, '
'and absorption residues. Possibly one more for a '
'fission residue.')
if not np.issubdtype(data.dtype, complex):
raise TypeError('Multipole data arrays must be complex dtype')
self._data = data
@ -364,6 +357,12 @@ class WindowedMultipole(EqualityMixin):
if not np.issubdtype(l_value.dtype, int):
raise TypeError('Multipole l_value arrays must be integer'
' dtype')
self._num_l = len(np.unique(l_value))
else:
self._num_l = None
self._l_value = l_value
@w_start.setter
@ -429,6 +428,7 @@ class WindowedMultipole(EqualityMixin):
format.
"""
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
@ -443,20 +443,12 @@ class WindowedMultipole(EqualityMixin):
'Python API expects version ' + WMP_VERSION)
group = h5file['nuclide']
out = cls()
# Read scalar values. Note that group['max_w'] is ignored.
length = group['length'].value
windows = group['windows'].value
out.num_l = group['num_l'].value
out.fit_order = group['fit_order'].value
out.fissionable = bool(group['fissionable'].value)
# Read scalars.
if group['formalism'].value == _FORM_MLBW:
out.formalism = 'MLBW'
out = cls('MLBW')
elif group['formalism'].value == _FORM_RM:
out.formalism = 'RM'
out = cls('RM')
else:
raise ValueError('Unrecognized/Unsupported R-matrix formalism')
@ -467,39 +459,36 @@ class WindowedMultipole(EqualityMixin):
# Read arrays.
err = "WMP '{}' array shape is not consistent with the '{}' value"
err = "WMP '{}' array shape is not consistent with the '{}' array shape"
out.data = group['data'].value
if out.data.shape[0] != length:
raise ValueError(err.format('data', 'length'))
out.l_value = group['l_value'].value
if out.l_value.shape[0] != out.data.shape[0]:
raise ValueError(err.format('l_value', 'data'))
out.pseudo_k0RS = group['pseudo_K0RS'].value
if out.pseudo_k0RS.shape[0] != out.num_l:
raise ValueError(err.format('pseudo_k0RS', 'num_l'))
out.l_value = group['l_value'].value
if out.l_value.shape[0] != length:
raise ValueError(err.format('l_value', 'length'))
raise ValueError(err.format('pseudo_k0RS', 'l_value'))
out.w_start = group['w_start'].value
if out.w_start.shape[0] != windows:
raise ValueError(err.format('w_start', 'windows'))
out.w_end = group['w_end'].value
if out.w_end.shape[0] != windows:
raise ValueError(err.format('w_end', 'windows'))
if out.w_end.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('w_end', 'w_start'))
out.broaden_poly = group['broaden_poly'].value.astype(np.bool)
if out.broaden_poly.shape[0] != windows:
raise ValueError(err.format('broaden_poly', 'windows'))
if out.broaden_poly.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('broaden_poly', 'w_start'))
out.curvefit = group['curvefit'].value
if out.curvefit.shape[0] != windows:
raise ValueError(err.format('curvefit', 'windows'))
if out.curvefit.shape[1] != out.fit_order + 1:
raise ValueError(err.format('curvefit', 'fit_order'))
if out.curvefit.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('curvefit', 'w_start'))
# Note that all the file 3 data (group['reactions/MT...']) are ignored.
# _broaden_wmp_polynomials assumes the curve fit has at least 3 terms.
if out.fit_order < 2:
raise ValueError("Windowed multipole is only supported for "
"curvefits with 3 or more terms.")
return out
@ -662,3 +651,47 @@ class WindowedMultipole(EqualityMixin):
fun = np.vectorize(lambda x: self._evaluate(x, T))
return fun(E)
def export_to_hdf5(self, path, libver='earliest'):
"""Export windowed multipole data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
libver : {'earliest', 'latest'}
Compatibility mode for the HDF5 file. 'latest' will produce files
that are less backwards compatible but have performance benefits.
"""
# Open file and write version.
with h5py.File(path, 'w', libver=libver) as f:
f.create_dataset('version', (1, ), dtype='S10')
f['version'][:] = WMP_VERSION.encode('ASCII')
# Make a nuclide group.
g = f.create_group('nuclide')
# Write scalars.
if self.formalism == 'MLBW':
g.create_dataset('formalism',
data=np.array(_FORM_MLBW, dtype=np.int32))
else:
# Assume RM.
g.create_dataset('formalism',
data=np.array(_FORM_RM, dtype=np.int32))
g.create_dataset('spacing', data=np.array(self.spacing))
g.create_dataset('sqrtAWR', data=np.array(self.sqrtAWR))
g.create_dataset('start_E', data=np.array(self.start_E))
g.create_dataset('end_E', data=np.array(self.end_E))
# Write arrays.
g.create_dataset('data', data=self.data)
g.create_dataset('l_value', data=self.l_value)
g.create_dataset('pseudo_K0RS', data=self.pseudo_k0RS)
g.create_dataset('w_start', data=self.w_start)
g.create_dataset('w_end', data=self.w_end)
g.create_dataset('broaden_poly',
data=self.broaden_poly.astype(np.int8))
g.create_dataset('curvefit', data=self.curvefit)

24
openmc/data/neutron.py
View file

@ -1,6 +1,6 @@
from __future__ import division, unicode_literals
import sys
from collections import OrderedDict, Iterable, Mapping, MutableMapping
from collections import OrderedDict
from collections.abc import Iterable, Mapping, MutableMapping
from io import StringIO
from itertools import chain
from math import log10
@ -10,13 +10,12 @@ import shutil
import tempfile
from warnings import warn
from six import string_types
import numpy as np
import h5py
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Library, Table, get_table
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnd_name
from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum, ResonancesWithBackground
@ -94,9 +93,7 @@ def _get_metadata(zaid, metastable_scheme='nndc'):
# Determine name
element = ATOMIC_SYMBOL[Z]
name = '{}{}'.format(element, mass_number)
if metastable > 0:
name += '_m{}'.format(metastable)
name = gnd_name(Z, mass_number, metastable)
return (name, element, Z, mass_number, metastable)
@ -245,7 +242,7 @@ class IncidentNeutron(EqualityMixin):
@name.setter
def name(self, name):
cv.check_type('name', name, string_types)
cv.check_type('name', name, str)
self._name = name
@property
@ -301,7 +298,7 @@ class IncidentNeutron(EqualityMixin):
def urr(self, urr):
cv.check_type('probability table dictionary', urr, MutableMapping)
for key, value in urr:
cv.check_type('probability table temperature', key, string_types)
cv.check_type('probability table temperature', key, str)
cv.check_type('probability tables', value, ProbabilityTables)
self._urr = urr
@ -842,10 +839,7 @@ class IncidentNeutron(EqualityMixin):
Incident neutron continuous-energy data
"""
# Create temporary directory -- it would be preferable to use
# TemporaryDirectory(), but it is only available in Python 3.2
tmpdir = tempfile.mkdtemp()
try:
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
@ -873,8 +867,4 @@ class IncidentNeutron(EqualityMixin):
data.energy['0K'] = xs.x
data[2].xs['0K'] = xs
finally:
# Get rid of temporary files
shutil.rmtree(tmpdir)
return data

8
openmc/data/njoy.py
View file

@ -1,4 +1,3 @@
from __future__ import print_function
import argparse
from collections import namedtuple
from io import StringIO
@ -150,10 +149,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
with open(input_filename, 'w') as f:
f.write(commands)
# Create temporary directory -- it would be preferable to use
# TemporaryDirectory(), but it is only available in Python 3.2
tmpdir = tempfile.mkdtemp()
try:
with tempfile.TemporaryDirectory() as tmpdir:
# Copy evaluations to appropriates 'tapes'
for tape_num, filename in tapein.items():
tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
@ -187,8 +183,6 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
if os.path.isfile(tmpfilename):
shutil.move(tmpfilename, filename)
finally:
shutil.rmtree(tmpdir)
def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):

5
openmc/data/product.py
View file

@ -1,9 +1,8 @@
from collections import Iterable
from collections.abc import Iterable
from io import StringIO
from numbers import Real
import sys
from six import string_types
import numpy as np
import openmc.checkvalue as cv
@ -113,7 +112,7 @@ class Product(EqualityMixin):
@particle.setter
def particle(self, particle):
cv.check_type('product particle type', particle, string_types)
cv.check_type('product particle type', particle, str)
self._particle = particle
@yield_.setter

6
openmc/data/reaction.py
View file

@ -1,11 +1,9 @@
from __future__ import division, unicode_literals
from collections import Iterable, Callable, MutableMapping
from collections.abc import Iterable, Callable, MutableMapping
from copy import deepcopy
from numbers import Real, Integral
from warnings import warn
from io import StringIO
from six import string_types
import numpy as np
import openmc.checkvalue as cv
@ -863,7 +861,7 @@ class Reaction(EqualityMixin):
def xs(self, xs):
cv.check_type('reaction cross section dictionary', xs, MutableMapping)
for key, value in xs.items():
cv.check_type('reaction cross section temperature', key, string_types)
cv.check_type('reaction cross section temperature', key, str)
cv.check_type('reaction cross section', value, Callable)
self._xs = xs

21
openmc/data/resonance.py
View file

@ -1,4 +1,5 @@
from collections import defaultdict, MutableSequence, Iterable
from collections import defaultdict
from collections.abc import MutableSequence, Iterable
import io
import numpy as np
@ -288,8 +289,8 @@ class MultiLevelBreitWigner(ResonanceRange):
"""
def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
super(MultiLevelBreitWigner, self).__init__(
target_spin, energy_min, energy_max, channel, scattering)
super().__init__(target_spin, energy_min, energy_max, channel,
scattering)
self.parameters = None
self.q_value = {}
self.atomic_weight_ratio = None
@ -490,8 +491,8 @@ class SingleLevelBreitWigner(MultiLevelBreitWigner):
"""
def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
super(SingleLevelBreitWigner, self).__init__(
target_spin, energy_min, energy_max, channel, scattering)
super().__init__(target_spin, energy_min, energy_max, channel,
scattering)
# Set resonance reconstruction function
if _reconstruct:
@ -549,8 +550,8 @@ class ReichMoore(ResonanceRange):
"""
def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
super(ReichMoore, self).__init__(
target_spin, energy_min, energy_max, channel, scattering)
super().__init__(target_spin, energy_min, energy_max, channel,
scattering)
self.parameters = None
self.angle_distribution = False
self.num_l_convergence = 0
@ -724,8 +725,7 @@ class RMatrixLimited(ResonanceRange):
"""
def __init__(self, energy_min, energy_max, particle_pairs, spin_groups):
super(RMatrixLimited, self).__init__(0.0, energy_min, energy_max,
None, None)
super().__init__(0.0, energy_min, energy_max, None, None)
self.reduced_width = False
self.formalism = 3
self.particle_pairs = particle_pairs
@ -931,8 +931,7 @@ class Unresolved(ResonanceRange):
"""
def __init__(self, target_spin, energy_min, energy_max, scatter):
super(Unresolved, self).__init__(
target_spin, energy_min, energy_max, None, scatter)
super().__init__(target_spin, energy_min, energy_max, None, scatter)
self.energies = None
self.parameters = None
self.add_to_background = False

10
openmc/data/thermal.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from difflib import get_close_matches
from numbers import Real
import itertools
@ -623,10 +623,7 @@ class ThermalScattering(EqualityMixin):
Thermal scattering data
"""
# Create temporary directory -- it would be preferable to use
# TemporaryDirectory(), but it is only available in Python 3.2
tmpdir = tempfile.mkdtemp()
try:
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
@ -638,8 +635,5 @@ class ThermalScattering(EqualityMixin):
data = cls.from_ace(lib.tables[0])
for table in lib.tables[1:]:
data.add_temperature_from_ace(table)
finally:
# Get rid of temporary files
shutil.rmtree(tmpdir)
return data

2
openmc/data/urr.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from numbers import Integral, Real
import numpy as np

24
openmc/deplete/__init__.py Normal file
View file

@ -0,0 +1,24 @@
"""
openmc.deplete
==============
A depletion front-end tool.
"""
from .dummy_comm import DummyCommunicator
try:
from mpi4py import MPI
comm = MPI.COMM_WORLD
have_mpi = True
except ImportError:
comm = DummyCommunicator()
have_mpi = False
from .nuclide import *
from .chain import *
from .operator import *
from .reaction_rates import *
from .abc import *
from .results import *
from .results_list import *
from .integrator import *

142
openmc/deplete/abc.py Normal file
View file

@ -0,0 +1,142 @@
"""function module.
This module contains the Operator class, which is then passed to an integrator
to run a full depletion simulation.
"""
from collections import namedtuple
import os
from pathlib import Path
from abc import ABCMeta, abstractmethod
from .chain import Chain
OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
OperatorResult.__doc__ = """\
Result of applying transport operator
Parameters
----------
k : float
Resulting eigenvalue
rates : openmc.deplete.ReactionRates
Resulting reaction rates
"""
try:
OperatorResult.k.__doc__ = None
OperatorResult.rates.__doc__ = None
except AttributeError:
# Can't set __doc__ on properties on Python 3.4
pass
class TransportOperator(metaclass=ABCMeta):
"""Abstract class defining a transport operator
Each depletion integrator is written to work with a generic transport
operator that takes a vector of material compositions and returns an
eigenvalue and reaction rates. This abstract class sets the requirements for
such a transport operator. Users should instantiate
:class:`openmc.deplete.Operator` rather than this class.
Parameters
----------
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
Attributes
----------
dilute_initial : float
Initial atom density to add for nuclides that are zero in initial
condition to ensure they exist in the decay chain. Only done for
nuclides with reaction rates. Defaults to 1.0e3.
"""
def __init__(self, chain_file=None):
self.dilute_initial = 1.0e3
self.output_dir = '.'
# Read depletion chain
if chain_file is None:
chain_file = os.environ.get("OPENMC_DEPLETE_CHAIN", None)
if chain_file is None:
raise IOError("No chain specified, either manually or in "
"environment variable OPENMC_DEPLETE_CHAIN.")
self.chain = Chain.from_xml(chain_file)
@abstractmethod
def __call__(self, vec, print_out=True):
"""Runs a simulation.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
print_out : bool, optional
Whether or not to print out time.
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
pass
def __enter__(self):
# Save current directory and move to specific output directory
self._orig_dir = os.getcwd()
if not self.output_dir.exists():
self.output_dir.mkdir() # exist_ok parameter is 3.5+
# In Python 3.6+, chdir accepts a Path directly
os.chdir(str(self.output_dir))
return self.initial_condition()
def __exit__(self, exc_type, exc_value, traceback):
self.finalize()
os.chdir(self._orig_dir)
@property
def output_dir(self):
return self._output_dir
@output_dir.setter
def output_dir(self, output_dir):
self._output_dir = Path(output_dir)
@abstractmethod
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
pass
@abstractmethod
def get_results_info(self):
"""Returns volume list, cell lists, and nuc lists.
Returns
-------
volume : list of float
Volumes corresponding to materials in burn_list
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all cell IDs to be burned. Used for sorting the simulation.
full_burn_list : list of int
All burnable materials in the geometry.
"""
pass
def finalize(self):
pass

214
openmc/deplete/atom_number.py Normal file
View file

@ -0,0 +1,214 @@
"""AtomNumber module.
An ndarray to store atom densities with string, integer, or slice indexing.
"""
from collections import OrderedDict
import numpy as np
class AtomNumber(object):
"""Stores local material compositions (atoms of each nuclide).
Parameters
----------
local_mats : list of str
Material IDs
nuclides : list of str
Nuclides to be tracked
volume : dict
Volume of each material in [cm^3]
n_nuc_burn : int
Number of nuclides to be burned.
Attributes
----------
index_mat : dict
A dictionary mapping material ID as string to index.
index_nuc : dict
A dictionary mapping nuclide name to index.
volume : numpy.ndarray
Volume of each material in [cm^3]. If a volume is not found, it defaults
to 1 so that reading density still works correctly.
number : numpy.ndarray
Array storing total atoms for each material/nuclide
materials : list of str
Material IDs as strings
nuclides : list of str
All nuclide names
burnable_nuclides : list of str
Burnable nuclides names. Used for sorting the simulation.
n_nuc_burn : int
Number of burnable nuclides.
n_nuc : int
Number of nuclides.
"""
def __init__(self, local_mats, nuclides, volume, n_nuc_burn):
self.index_mat = OrderedDict((mat, i) for i, mat in enumerate(local_mats))
self.index_nuc = OrderedDict((nuc, i) for i, nuc in enumerate(nuclides))
self.volume = np.ones(len(local_mats))
for mat, val in volume.items():
if mat in self.index_mat:
ind = self.index_mat[mat]
self.volume[ind] = val
self.n_nuc_burn = n_nuc_burn
self.number = np.zeros((len(local_mats), len(nuclides)))
def __getitem__(self, pos):
"""Retrieves total atom number from AtomNumber.
Parameters
----------
pos : tuple
A two-length tuple containing a material index and a nuc index.
These indexes can be strings (which get converted to integers via
the dictionaries), integers used directly, or slices.
Returns
-------
numpy.ndarray
The value indexed from self.number.
"""
mat, nuc = pos
if isinstance(mat, str):
mat = self.index_mat[mat]
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
return self.number[mat, nuc]
def __setitem__(self, pos, val):
"""Sets total atom number into AtomNumber.
Parameters
----------
pos : tuple
A two-length tuple containing a material index and a nuc index.
These indexes can be strings (which get converted to integers via
the dictionaries), integers used directly, or slices.
val : float
The value [atom] to set the array to.
"""
mat, nuc = pos
if isinstance(mat, str):
mat = self.index_mat[mat]
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
self.number[mat, nuc] = val
@property
def materials(self):
return self.index_mat.keys()
@property
def nuclides(self):
return self.index_nuc.keys()
@property
def n_nuc(self):
return len(self.index_nuc)
@property
def burnable_nuclides(self):
return [nuc for nuc, ind in self.index_nuc.items()
if ind < self.n_nuc_burn]
def get_atom_density(self, mat, nuc):
"""Accesses atom density instead of total number.
Parameters
----------
mat : str, int or slice
Material index.
nuc : str, int or slice
Nuclide index.
Returns
-------
numpy.ndarray
Density in [atom/cm^3]
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
return self[mat, nuc] / self.volume[mat]
def set_atom_density(self, mat, nuc, val):
"""Sets atom density instead of total number.
Parameters
----------
mat : str, int or slice
Material index.
nuc : str, int or slice
Nuclide index.
val : numpy.ndarray
Array of densities to set in [atom/cm^3]
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
self[mat, nuc] = val * self.volume[mat]
def get_mat_slice(self, mat):
"""Gets atom quantity indexed by mats for all burned nuclides
Parameters
----------
mat : str, int or slice
Material index.
Returns
-------
numpy.ndarray
The slice requested in [atom].
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
return self[mat, :self.n_nuc_burn]
def set_mat_slice(self, mat, val):
"""Sets atom quantity indexed by mats for all burned nuclides
Parameters
----------
mat : str, int or slice
Material index.
val : numpy.ndarray
The slice to set in [atom]
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
self[mat, :self.n_nuc_burn] = val
def set_density(self, total_density):
"""Sets density.
Sets the density in the exact same order as total_density_list outputs,
allowing for internal consistency
Parameters
----------
total_density : list of numpy.ndarray
Total atoms.
"""
for i, density_slice in enumerate(total_density):
self.set_mat_slice(i, density_slice)

440
openmc/deplete/chain.py Normal file
View file

@ -0,0 +1,440 @@
"""chain module.
This module contains information about a depletion chain. A depletion chain is
loaded from an .xml file and all the nuclides are linked together.
"""
from collections import OrderedDict, defaultdict
from io import StringIO
from itertools import chain
import math
import re
import os
# Try to use lxml if it is available. It preserves the order of attributes and
# provides a pretty-printer by default. If not available, use OpenMC function to
# pretty print.
try:
import lxml.etree as ET
_have_lxml = True
except ImportError:
import xml.etree.ElementTree as ET
_have_lxml = False
import scipy.sparse as sp
import openmc.data
from openmc.clean_xml import clean_xml_indentation
from .nuclide import Nuclide, DecayTuple, ReactionTuple
# tuple of (reaction name, possible MT values, (dA, dZ)) where dA is the change
# in the mass number and dZ is the change in the atomic number
_REACTIONS = [
('(n,2n)', set(chain([16], range(875, 892))), (-1, 0)),
('(n,3n)', {17}, (-2, 0)),
('(n,4n)', {37}, (-3, 0)),
('(n,gamma)', {102}, (1, 0)),
('(n,p)', set(chain([103], range(600, 650))), (0, -1)),
('(n,a)', set(chain([107], range(800, 850))), (-3, -2))
]
def replace_missing(product, decay_data):
"""Replace missing product with suitable decay daughter.
Parameters
----------
product : str
Name of product in GND format, e.g. 'Y86_m1'.
decay_data : dict
Dictionary of decay data
Returns
-------
product : str
Replacement for missing product in GND format.
"""
# Determine atomic number, mass number, and metastable state
Z, A, state = openmc.data.zam(product)
symbol = openmc.data.ATOMIC_SYMBOL[Z]
# Replace neutron with proton
if Z == 0 and A == 1:
return 'H1'
# First check if ground state is available
if state:
product = '{}{}'.format(symbol, A)
# Find isotope with longest half-life
half_life = 0.0
for nuclide, data in decay_data.items():
m = re.match(r'{}(\d+)(?:_m\d+)?'.format(symbol), nuclide)
if m:
# If we find a stable nuclide, stop search
if data.nuclide['stable']:
mass_longest_lived = int(m.group(1))
break
if data.half_life.nominal_value > half_life:
mass_longest_lived = int(m.group(1))
half_life = data.half_life.nominal_value
# If mass number of longest-lived isotope is less than that of missing
# product, assume it undergoes beta-. Otherwise assume beta+.
beta_minus = (mass_longest_lived < A)
# Iterate until we find an existing nuclide
while product not in decay_data:
if Z > 98:
Z -= 2
A -= 4
else:
if beta_minus:
Z += 1
else:
Z -= 1
product = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
return product
class Chain(object):
"""Full representation of a depletion chain.
A depletion chain can be created by using the :meth:`from_endf` method which
requires a list of ENDF incident neutron, decay, and neutron fission product
yield sublibrary files. The depletion chain used during a depletion
simulation is indicated by either an argument to
:class:`openmc.deplete.Operator` or through the
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable.
Attributes
----------
nuclides : list of openmc.deplete.Nuclide
Nuclides present in the chain.
reactions : list of str
Reactions that are tracked in the depletion chain
nuclide_dict : OrderedDict of str to int
Maps a nuclide name to an index in nuclides.
"""
def __init__(self):
self.nuclides = []
self.reactions = []
self.nuclide_dict = OrderedDict()
def __contains__(self, nuclide):
return nuclide in self.nuclide_dict
def __getitem__(self, name):
"""Get a Nuclide by name."""
return self.nuclides[self.nuclide_dict[name]]
def __len__(self):
"""Number of nuclides in chain."""
return len(self.nuclides)
@classmethod
def from_endf(cls, decay_files, fpy_files, neutron_files):
"""Create a depletion chain from ENDF files.
Parameters
----------
decay_files : list of str
List of ENDF decay sub-library files
fpy_files : list of str
List of ENDF neutron-induced fission product yield sub-library files
neutron_files : list of str
List of ENDF neutron reaction sub-library files
"""
chain = cls()
# Create dictionary mapping target to filename
print('Processing neutron sub-library files...')
reactions = {}
for f in neutron_files:
evaluation = openmc.data.endf.Evaluation(f)
name = evaluation.gnd_name
reactions[name] = {}
for mf, mt, nc, mod in evaluation.reaction_list:
if mf == 3:
file_obj = StringIO(evaluation.section[3, mt])
openmc.data.endf.get_head_record(file_obj)
q_value = openmc.data.endf.get_cont_record(file_obj)[1]
reactions[name][mt] = q_value
# Determine what decay and FPY nuclides are available
print('Processing decay sub-library files...')
decay_data = {}
for f in decay_files:
data = openmc.data.Decay(f)
# Skip decay data for neutron itself
if data.nuclide['atomic_number'] == 0:
continue
decay_data[data.nuclide['name']] = data
print('Processing fission product yield sub-library files...')
fpy_data = {}
for f in fpy_files:
data = openmc.data.FissionProductYields(f)
fpy_data[data.nuclide['name']] = data
print('Creating depletion_chain...')
missing_daughter = []
missing_rx_product = []
missing_fpy = []
missing_fp = []
for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)):
data = decay_data[parent]
nuclide = Nuclide()
nuclide.name = parent
chain.nuclides.append(nuclide)
chain.nuclide_dict[parent] = idx
if not data.nuclide['stable'] and data.half_life.nominal_value != 0.0:
nuclide.half_life = data.half_life.nominal_value
nuclide.decay_energy = sum(E.nominal_value for E in
data.average_energies.values())
sum_br = 0.0
for i, mode in enumerate(data.modes):
type_ = ','.join(mode.modes)
if mode.daughter in decay_data:
target = mode.daughter
else:
print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter))
target = replace_missing(mode.daughter, decay_data)
# Write branching ratio, taking care to ensure sum is unity
br = mode.branching_ratio.nominal_value
sum_br += br
if i == len(data.modes) - 1 and sum_br != 1.0:
br = 1.0 - sum(m.branching_ratio.nominal_value
for m in data.modes[:-1])
# Append decay mode
nuclide.decay_modes.append(DecayTuple(type_, target, br))
if parent in reactions:
reactions_available = set(reactions[parent].keys())
for name, mts, changes in _REACTIONS:
if mts & reactions_available:
delta_A, delta_Z = changes
A = data.nuclide['mass_number'] + delta_A
Z = data.nuclide['atomic_number'] + delta_Z
daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
if name not in chain.reactions:
chain.reactions.append(name)
if daughter not in decay_data:
missing_rx_product.append((parent, name, daughter))
# Store Q value
for mt in sorted(mts):
if mt in reactions[parent]:
q_value = reactions[parent][mt]
break
else:
q_value = 0.0
nuclide.reactions.append(ReactionTuple(
name, daughter, q_value, 1.0))
if any(mt in reactions_available for mt in [18, 19, 20, 21, 38]):
if parent in fpy_data:
q_value = reactions[parent][18]
nuclide.reactions.append(
ReactionTuple('fission', 0, q_value, 1.0))
if 'fission' not in chain.reactions:
chain.reactions.append('fission')
else:
missing_fpy.append(parent)
if parent in fpy_data:
fpy = fpy_data[parent]
if fpy.energies is not None:
nuclide.yield_energies = fpy.energies
else:
nuclide.yield_energies = [0.0]
for E, table in zip(nuclide.yield_energies, fpy.independent):
yield_replace = 0.0
yields = defaultdict(float)
for product, y in table.items():
# Handle fission products that have no decay data available
if product not in decay_data:
daughter = replace_missing(product, decay_data)
product = daughter
yield_replace += y.nominal_value
yields[product] += y.nominal_value
if yield_replace > 0.0:
missing_fp.append((parent, E, yield_replace))
nuclide.yield_data[E] = []
for k in sorted(yields, key=openmc.data.zam):
nuclide.yield_data[E].append((k, yields[k]))
# Display warnings
if missing_daughter:
print('The following decay modes have daughters with no decay data:')
for mode in missing_daughter:
print(' {}'.format(mode))
print('')
if missing_rx_product:
print('The following reaction products have no decay data:')
for vals in missing_rx_product:
print('{} {} -> {}'.format(*vals))
print('')
if missing_fpy:
print('The following fissionable nuclides have no fission product yields:')
for parent in missing_fpy:
print(' ' + parent)
print('')
if missing_fp:
print('The following nuclides have fission products with no decay data:')
for vals in missing_fp:
print(' {}, E={} eV (total yield={})'.format(*vals))
return chain
@classmethod
def from_xml(cls, filename):
"""Reads a depletion chain XML file.
Parameters
----------
filename : str
The path to the depletion chain XML file.
"""
chain = cls()
# Load XML tree
root = ET.parse(str(filename))
for i, nuclide_elem in enumerate(root.findall('nuclide')):
nuc = Nuclide.from_xml(nuclide_elem)
chain.nuclide_dict[nuc.name] = i
# Check for reaction paths
for rx in nuc.reactions:
if rx.type not in chain.reactions:
chain.reactions.append(rx.type)
chain.nuclides.append(nuc)
return chain
def export_to_xml(self, filename):
"""Writes a depletion chain XML file.
Parameters
----------
filename : str
The path to the depletion chain XML file.
"""
root_elem = ET.Element('depletion_chain')
for nuclide in self.nuclides:
root_elem.append(nuclide.to_xml_element())
tree = ET.ElementTree(root_elem)
if _have_lxml:
tree.write(str(filename), encoding='utf-8', pretty_print=True)
else:
clean_xml_indentation(root_elem)
tree.write(str(filename), encoding='utf-8')
def form_matrix(self, rates):
"""Forms depletion matrix.
Parameters
----------
rates : numpy.ndarray
2D array indexed by (nuclide, reaction)
Returns
-------
scipy.sparse.csr_matrix
Sparse matrix representing depletion.
"""
matrix = defaultdict(float)
reactions = set()
for i, nuc in enumerate(self.nuclides):
if nuc.n_decay_modes != 0:
# Decay paths
# Loss
decay_constant = math.log(2) / nuc.half_life
if decay_constant != 0.0:
matrix[i, i] -= decay_constant
# Gain
for _, target, branching_ratio in nuc.decay_modes:
# Allow for total annihilation for debug purposes
if target != 'Nothing':
branch_val = branching_ratio * decay_constant
if branch_val != 0.0:
k = self.nuclide_dict[target]
matrix[k, i] += branch_val
if nuc.name in rates.index_nuc:
# Extract all reactions for this nuclide in this cell
nuc_ind = rates.index_nuc[nuc.name]
nuc_rates = rates[nuc_ind, :]
for r_type, target, _, br in nuc.reactions:
# Extract reaction index, and then final reaction rate
r_id = rates.index_rx[r_type]
path_rate = nuc_rates[r_id]
# Loss term -- make sure we only count loss once for
# reactions with branching ratios
if r_type not in reactions:
reactions.add(r_type)
if path_rate != 0.0:
matrix[i, i] -= path_rate
# Gain term; allow for total annihilation for debug purposes
if target != 'Nothing':
if r_type != 'fission':
if path_rate != 0.0:
k = self.nuclide_dict[target]
matrix[k, i] += path_rate * br
else:
# Assume that we should always use thermal fission
# yields. At some point it would be nice to account
# for the energy-dependence..
energy, data = sorted(nuc.yield_data.items())[0]
for product, y in data:
yield_val = y * path_rate
if yield_val != 0.0:
k = self.nuclide_dict[product]
matrix[k, i] += yield_val
# Clear set of reactions
reactions.clear()
# Use DOK matrix as intermediate representation, then convert to CSR and return
n = len(self)
matrix_dok = sp.dok_matrix((n, n))
dict.update(matrix_dok, matrix)
return matrix_dok.tocsr()

27
openmc/deplete/dummy_comm.py Normal file
View file

@ -0,0 +1,27 @@
class DummyCommunicator(object):
rank = 0
size = 1
def allgather(self, sendobj):
return [sendobj]
def allreduce(self, sendobj, op=None):
return sendobj
def barrier(self):
pass
def bcast(self, obj, root=0):
return obj
def gather(self, sendobj, root=0):
return [sendobj]
def py2f(self):
return 0
def reduce(self, sendobj, op=None, root=0):
return sendobj
def scatter(self, sendobj, root=0):
return sendobj[0]

10
openmc/deplete/integrator/__init__.py Normal file
View file

@ -0,0 +1,10 @@
"""
Integrator
===========
The integrator subcomponents.
"""
from .cecm import *
from .cram import *
from .predictor import *

78
openmc/deplete/integrator/cecm.py Normal file
View file

@ -0,0 +1,78 @@
"""The CE/CM integrator."""
import copy
from collections.abc import Iterable
from .cram import deplete
from ..results import Results
def cecm(operator, timesteps, power, print_out=True):
r"""Deplete using the CE/CM algorithm.
Implements the second order `CE/CM predictor-corrector algorithm
<https://doi.org/10.13182/NSE14-92>`_. This algorithm is mathematically
defined as:
.. math::
y' &= A(y, t) y(t)
A_p &= A(y_n, t_n)
y_m &= \text{expm}(A_p h/2) y_n
A_c &= A(y_m, t_n + h/2)
y_{n+1} &= \text{expm}(A_c h) y_n
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
power : float or iterable of float
Power of the reactor in [W]. A single value indicates that the power is
constant over all timesteps. An iterable indicates potentially different
power levels for each timestep. For a 2D problem, the power can be given
in [W/cm] as long as the "volume" assigned to a depletion material is
actually an area in [cm^2].
print_out : bool, optional
Whether or not to print out time.
"""
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
chain = operator.chain
t = 0.0
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep reaction rates
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
# Deplete for first half of timestep
x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)
# Get middle-of-timestep reaction rates
x.append(x_middle)
op_results.append(operator(x_middle, p))
# Deplete for full timestep using beginning-of-step materials
x_end = deplete(chain, x[0], op_results[1], dt, print_out)
# Create results, write to disk
Results.save(operator, x, op_results, [t, t + dt], i)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], len(timesteps))

227
openmc/deplete/integrator/cram.py Normal file
View file

@ -0,0 +1,227 @@
"""Chebyshev Rational Approximation Method module
Implements two different forms of CRAM for use in openmc.deplete.
"""
from itertools import repeat
from multiprocessing import Pool
import time
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as sla
from .. import comm
def deplete(chain, x, op_result, dt, print_out):
"""Deplete materials using given reaction rates for a specified time
Parameters
----------
chain : openmc.deplete.Chain
Depletion chain
x : list of numpy.ndarray
Atom number vectors for each material
op_result : openmc.deplete.OperatorResult
Result of applying transport operator (contains reaction rates)
dt : float
Time in [s] to deplete for
print_out : bool
Whether to show elapsed time
Returns
-------
x_result : list of numpy.ndarray
Updated atom number vectors for each material
"""
t_start = time.time()
# Set up iterators
n_mats = len(x)
chains = repeat(chain, n_mats)
vecs = (x[i] for i in range(n_mats))
rates = (op_result.rates[i, :, :] for i in range(n_mats))
dts = repeat(dt, n_mats)
# Use multiprocessing pool to distribute work
with Pool() as pool:
iters = zip(chains, vecs, rates, dts)
x_result = list(pool.starmap(_cram_wrapper, iters))
t_end = time.time()
if comm.rank == 0:
if print_out:
print("Time to matexp: ", t_end - t_start)
return x_result
def _cram_wrapper(chain, n0, rates, dt):
"""Wraps depletion matrix creation / CRAM solve for multiprocess execution
Parameters
----------
chain : DepletionChain
Depletion chain used to construct the burnup matrix
n0 : numpy.array
Vector to operate a matrix exponent on.
rates : numpy.ndarray
2D array indexed by nuclide then by cell.
dt : float
Time to integrate to.
Returns
-------
numpy.array
Results of the matrix exponent.
"""
A = chain.form_matrix(rates)
return CRAM48(A, n0, dt)
def CRAM16(A, n0, dt):
"""Chebyshev Rational Approximation Method, order 16
Algorithm is the 16th order Chebyshev Rational Approximation Method,
implemented in the more stable `incomplete partial fraction (IPF)
<https://doi.org/10.13182/NSE15-26>`_ form.
Parameters
----------
A : scipy.linalg.csr_matrix
Matrix to take exponent of.
n0 : numpy.array
Vector to operate a matrix exponent on.
dt : float
Time to integrate to.
Returns
-------
numpy.array
Results of the matrix exponent.
"""
alpha = np.array([+2.124853710495224e-16,
+5.464930576870210e+3 - 3.797983575308356e+4j,
+9.045112476907548e+1 - 1.115537522430261e+3j,
+2.344818070467641e+2 - 4.228020157070496e+2j,
+9.453304067358312e+1 - 2.951294291446048e+2j,
+7.283792954673409e+2 - 1.205646080220011e+5j,
+3.648229059594851e+1 - 1.155509621409682e+2j,
+2.547321630156819e+1 - 2.639500283021502e+1j,
+2.394538338734709e+1 - 5.650522971778156e+0j],
dtype=np.complex128)
theta = np.array([+0.0,
+3.509103608414918 + 8.436198985884374j,
+5.948152268951177 + 3.587457362018322j,
-5.264971343442647 + 16.22022147316793j,
+1.419375897185666 + 10.92536348449672j,
+6.416177699099435 + 1.194122393370139j,
+4.993174737717997 + 5.996881713603942j,
-1.413928462488886 + 13.49772569889275j,
-10.84391707869699 + 19.27744616718165j],
dtype=np.complex128)
n = A.shape[0]
alpha0 = 2.124853710495224e-16
k = 8
y = np.array(n0, dtype=np.float64)
for l in range(1, k+1):
y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
y *= alpha0
return y
def CRAM48(A, n0, dt):
"""Chebyshev Rational Approximation Method, order 48
Algorithm is the 48th order Chebyshev Rational Approximation Method,
implemented in the more stable `incomplete partial fraction (IPF)
<https://doi.org/10.13182/NSE15-26>`_ form.
Parameters
----------
A : scipy.linalg.csr_matrix
Matrix to take exponent of.
n0 : numpy.array
Vector to operate a matrix exponent on.
dt : float
Time to integrate to.
Returns
-------
numpy.array
Results of the matrix exponent.
"""
theta_r = np.array([-4.465731934165702e+1, -5.284616241568964e+0,
-8.867715667624458e+0, +3.493013124279215e+0,
+1.564102508858634e+1, +1.742097597385893e+1,
-2.834466755180654e+1, +1.661569367939544e+1,
+8.011836167974721e+0, -2.056267541998229e+0,
+1.449208170441839e+1, +1.853807176907916e+1,
+9.932562704505182e+0, -2.244223871767187e+1,
+8.590014121680897e-1, -1.286192925744479e+1,
+1.164596909542055e+1, +1.806076684783089e+1,
+5.870672154659249e+0, -3.542938819659747e+1,
+1.901323489060250e+1, +1.885508331552577e+1,
-1.734689708174982e+1, +1.316284237125190e+1])
theta_i = np.array([+6.233225190695437e+1, +4.057499381311059e+1,
+4.325515754166724e+1, +3.281615453173585e+1,
+1.558061616372237e+1, +1.076629305714420e+1,
+5.492841024648724e+1, +1.316994930024688e+1,
+2.780232111309410e+1, +3.794824788914354e+1,
+1.799988210051809e+1, +5.974332563100539e+0,
+2.532823409972962e+1, +5.179633600312162e+1,
+3.536456194294350e+1, +4.600304902833652e+1,
+2.287153304140217e+1, +8.368200580099821e+0,
+3.029700159040121e+1, +5.834381701800013e+1,
+1.194282058271408e+0, +3.583428564427879e+0,
+4.883941101108207e+1, +2.042951874827759e+1])
theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
alpha_r = np.array([+6.387380733878774e+2, +1.909896179065730e+2,
+4.236195226571914e+2, +4.645770595258726e+2,
+7.765163276752433e+2, +1.907115136768522e+3,
+2.909892685603256e+3, +1.944772206620450e+2,
+1.382799786972332e+5, +5.628442079602433e+3,
+2.151681283794220e+2, +1.324720240514420e+3,
+1.617548476343347e+4, +1.112729040439685e+2,
+1.074624783191125e+2, +8.835727765158191e+1,
+9.354078136054179e+1, +9.418142823531573e+1,
+1.040012390717851e+2, +6.861882624343235e+1,
+8.766654491283722e+1, +1.056007619389650e+2,
+7.738987569039419e+1, +1.041366366475571e+2])
alpha_i = np.array([-6.743912502859256e+2, -3.973203432721332e+2,
-2.041233768918671e+3, -1.652917287299683e+3,
-1.783617639907328e+4, -5.887068595142284e+4,
-9.953255345514560e+3, -1.427131226068449e+3,
-3.256885197214938e+6, -2.924284515884309e+4,
-1.121774011188224e+3, -6.370088443140973e+4,
-1.008798413156542e+6, -8.837109731680418e+1,
-1.457246116408180e+2, -6.388286188419360e+1,
-2.195424319460237e+2, -6.719055740098035e+2,
-1.693747595553868e+2, -1.177598523430493e+1,
-4.596464999363902e+3, -1.738294585524067e+3,
-4.311715386228984e+1, -2.777743732451969e+2])
alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
n = A.shape[0]
alpha0 = 2.258038182743983e-47
k = 24
y = np.array(n0, dtype=np.float64)
for l in range(k):
y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
y *= alpha0
return y

66
openmc/deplete/integrator/predictor.py Normal file
View file

@ -0,0 +1,66 @@
"""First-order predictor algorithm."""
import copy
from collections.abc import Iterable
from .cram import deplete
from ..results import Results
def predictor(operator, timesteps, power, print_out=True):
r"""Deplete using a first-order predictor algorithm.
Implements the first-order predictor algorithm. This algorithm is
mathematically defined as:
.. math::
y' &= A(y, t) y(t)
A_p &= A(y_n, t_n)
y_{n+1} &= \text{expm}(A_p h) y_n
Parameters
----------
operator : openmc.deplete.TransportOperator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
power : float or iterable of float
Power of the reactor in [W]. A single value indicates that the power is
constant over all timesteps. An iterable indicates potentially different
power levels for each timestep. For a 2D problem, the power can be given
in [W/cm] as long as the "volume" assigned to a depletion material is
actually an area in [cm^2].
print_out : bool, optional
Whether or not to print out time.
"""
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
chain = operator.chain
t = 0.0
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep reaction rates
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t + dt], i)
# Deplete for full timestep
x_end = deplete(chain, x[0], op_results[0], dt, print_out)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
Results.save(operator, x, op_results, [t, t], len(timesteps))

219
openmc/deplete/nuclide.py Normal file
View file

@ -0,0 +1,219 @@
"""Nuclide module.
Contains the per-nuclide components of a depletion chain.
"""
from collections import namedtuple
try:
import lxml.etree as ET
except ImportError:
import xml.etree.ElementTree as ET
DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio')
DecayTuple.__doc__ = """\
Decay mode information
Parameters
----------
type : str
Type of the decay mode, e.g., 'beta-'
target : str
Nuclide resulting from decay
branching_ratio : float
Branching ratio of the decay mode
"""
try:
DecayTuple.type.__doc__ = None
DecayTuple.target.__doc__ = None
DecayTuple.branching_ratio.__doc__ = None
except AttributeError:
# Can't set __doc__ on properties on Python 3.4
pass
ReactionTuple = namedtuple('ReactionTuple', 'type target Q branching_ratio')
ReactionTuple.__doc__ = """\
Transmutation reaction information
Parameters
----------
type : str
Type of the reaction, e.g., 'fission'
target : str
nuclide resulting from reaction
Q : float
Q value of the reaction in [eV]
branching_ratio : float
Branching ratio of the reaction
"""
try:
ReactionTuple.type.__doc__ = None
ReactionTuple.target.__doc__ = None
ReactionTuple.Q.__doc__ = None
ReactionTuple.branching_ratio.__doc__ = None
except AttributeError:
pass
class Nuclide(object):
"""Decay modes, reactions, and fission yields for a single nuclide.
Attributes
----------
name : str
Name of nuclide.
half_life : float
Half life of nuclide in [s].
decay_energy : float
Energy deposited from decay in [eV].
n_decay_modes : int
Number of decay pathways.
decay_modes : list of openmc.deplete.DecayTuple
Decay mode information. Each element of the list is a named tuple with
attributes 'type', 'target', and 'branching_ratio'.
n_reaction_paths : int
Number of possible reaction pathways.
reactions : list of openmc.deplete.ReactionTuple
Reaction information. Each element of the list is a named tuple with
attribute 'type', 'target', 'Q', and 'branching_ratio'.
yield_data : dict of float to list
Maps tabulated energy to list of (product, yield) for all
neutron-induced fission products.
yield_energies : list of float
Energies at which fission product yiels exist
"""
def __init__(self):
# Information about the nuclide
self.name = None
self.half_life = None
self.decay_energy = 0.0
# Decay paths
self.decay_modes = []
# Reaction paths
self.reactions = []
# Neutron fission yields, if present
self.yield_data = {}
self.yield_energies = []
@property
def n_decay_modes(self):
return len(self.decay_modes)
@property
def n_reaction_paths(self):
return len(self.reactions)
@classmethod
def from_xml(cls, element):
"""Read nuclide from an XML element.
Parameters
----------
element : xml.etree.ElementTree.Element
XML element to write nuclide data to
Returns
-------
nuc : openmc.deplete.Nuclide
Instance of a nuclide
"""
nuc = cls()
nuc.name = element.get('name')
# Check for half-life
if 'half_life' in element.attrib:
nuc.half_life = float(element.get('half_life'))
nuc.decay_energy = float(element.get('decay_energy', '0'))
# Check for decay paths
for decay_elem in element.iter('decay'):
d_type = decay_elem.get('type')
target = decay_elem.get('target')
branching_ratio = float(decay_elem.get('branching_ratio'))
nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
# Check for reaction paths
for reaction_elem in element.iter('reaction'):
r_type = reaction_elem.get('type')
Q = float(reaction_elem.get('Q', '0'))
branching_ratio = float(reaction_elem.get('branching_ratio', '1'))
# If the type is not fission, get target and Q value, otherwise
# just set null values
if r_type != 'fission':
target = reaction_elem.get('target')
else:
target = None
# Append reaction
nuc.reactions.append(ReactionTuple(
r_type, target, Q, branching_ratio))
fpy_elem = element.find('neutron_fission_yields')
if fpy_elem is not None:
for yields_elem in fpy_elem.iter('fission_yields'):
E = float(yields_elem.get('energy'))
products = yields_elem.find('products').text.split()
yields = [float(y) for y in
yields_elem.find('data').text.split()]
nuc.yield_data[E] = list(zip(products, yields))
nuc.yield_energies = list(sorted(nuc.yield_data.keys()))
return nuc
def to_xml_element(self):
"""Write nuclide to XML element.
Returns
-------
elem : xml.etree.ElementTree.Element
XML element to write nuclide data to
"""
elem = ET.Element('nuclide')
elem.set('name', self.name)
if self.half_life is not None:
elem.set('half_life', str(self.half_life))
elem.set('decay_modes', str(len(self.decay_modes)))
elem.set('decay_energy', str(self.decay_energy))
for mode, daughter, br in self.decay_modes:
mode_elem = ET.SubElement(elem, 'decay')
mode_elem.set('type', mode)
mode_elem.set('target', daughter)
mode_elem.set('branching_ratio', str(br))
elem.set('reactions', str(len(self.reactions)))
for rx, daughter, Q, br in self.reactions:
rx_elem = ET.SubElement(elem, 'reaction')
rx_elem.set('type', rx)
rx_elem.set('Q', str(Q))
if rx != 'fission':
rx_elem.set('target', daughter)
if br != 1.0:
rx_elem.set('branching_ratio', str(br))
if self.yield_data:
fpy_elem = ET.SubElement(elem, 'neutron_fission_yields')
energy_elem = ET.SubElement(fpy_elem, 'energies')
energy_elem.text = ' '.join(str(E) for E in self.yield_energies)
for E in self.yield_energies:
yields_elem = ET.SubElement(fpy_elem, 'fission_yields')
yields_elem.set('energy', str(E))
products_elem = ET.SubElement(yields_elem, 'products')
products_elem.text = ' '.join(x[0] for x in self.yield_data[E])
data_elem = ET.SubElement(yields_elem, 'data')
data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E])
return elem

562
openmc/deplete/operator.py Normal file
View file

@ -0,0 +1,562 @@
"""OpenMC transport operator
This module implements a transport operator for OpenMC so that it can be used by
depletion integrators. The implementation makes use of the Python bindings to
OpenMC's C API so that reading tally results and updating material number
densities is all done in-memory instead of through the filesystem.
"""
import copy
from collections import OrderedDict
from itertools import chain
import os
import time
import xml.etree.ElementTree as ET
import h5py
import numpy as np
import openmc
import openmc.capi
from openmc.data import JOULE_PER_EV
from . import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .reaction_rates import ReactionRates
def _distribute(items):
"""Distribute items across MPI communicator
Parameters
----------
items : list
List of items of distribute
Returns
-------
list
Items assigned to process that called
"""
min_size, extra = divmod(len(items), comm.size)
j = 0
for i in range(comm.size):
chunk_size = min_size + int(i < extra)
if comm.rank == i:
return items[j:j + chunk_size]
j += chunk_size
class Operator(TransportOperator):
"""OpenMC transport operator for depletion.
Instances of this class can be used to perform depletion using OpenMC as the
transport operator. Normally, a user needn't call methods of this class
directly. Instead, an instance of this class is passed to an integrator
function, such as :func:`openmc.deplete.integrator.cecm`.
Parameters
----------
geometry : openmc.Geometry
OpenMC geometry object
settings : openmc.Settings
OpenMC Settings object
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
Attributes
----------
geometry : openmc.Geometry
OpenMC geometry object
settings : openmc.Settings
OpenMC settings object
dilute_initial : float
Initial atom density to add for nuclides that are zero in initial
condition to ensure they exist in the decay chain. Only done for
nuclides with reaction rates. Defaults to 1.0e3.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
Whether or not to round output to OpenMC to 8 digits.
Useful in testing, as OpenMC is incredibly sensitive to exact values.
number : openmc.deplete.AtomNumber
Total number of atoms in simulation.
nuclides_with_data : set of str
A set listing all unique nuclides available from cross_sections.xml.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
reaction_rates : openmc.deplete.ReactionRates
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
local_mats : list of str
All burnable material IDs being managed by a single process
"""
def __init__(self, geometry, settings, chain_file=None):
super().__init__(chain_file)
self.round_number = False
self.settings = settings
self.geometry = geometry
# Clear out OpenMC, create task lists, distribute
openmc.reset_auto_ids()
self.burnable_mats, volume, nuclides = self._get_burnable_mats()
self.local_mats = _distribute(self.burnable_mats)
# Determine which nuclides have incident neutron data
self.nuclides_with_data = self._get_nuclides_with_data()
self._burnable_nucs = [nuc for nuc in self.nuclides_with_data
if nuc in self.chain]
# Extract number densities from the geometry
self._extract_number(self.local_mats, volume, nuclides)
# Create reaction rates array
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
def __call__(self, vec, power, print_out=True):
"""Runs a simulation.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
power : float
Power of the reactor in [W]
print_out : bool, optional
Whether or not to print out time.
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
# Prevent OpenMC from complaining about re-creating tallies
openmc.reset_auto_ids()
# Update status
self.number.set_density(vec)
time_start = time.time()
# Update material compositions and tally nuclides
self._update_materials()
self._tally.nuclides = self._get_tally_nuclides()
# Run OpenMC
openmc.capi.reset()
openmc.capi.run()
time_openmc = time.time()
# Extract results
op_result = self._unpack_tallies_and_normalize(power)
if comm.rank == 0:
time_unpack = time.time()
if print_out:
print("Time to openmc: ", time_openmc - time_start)
print("Time to unpack: ", time_unpack - time_openmc)
return copy.deepcopy(op_result)
def _get_burnable_mats(self):
"""Determine depletable materials, volumes, and nuclids
Returns
-------
burnable_mats : list of str
List of burnable material IDs
volume : OrderedDict of str to float
Volume of each material in [cm^3]
nuclides : list of str
Nuclides in order of how they'll appear in the simulation.
"""
burnable_mats = set()
model_nuclides = set()
volume = OrderedDict()
# Iterate once through the geometry to get dictionaries
for mat in self.geometry.get_all_materials().values():
for nuclide in mat.get_nuclides():
model_nuclides.add(nuclide)
if mat.depletable:
burnable_mats.add(str(mat.id))
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
volume[str(mat.id)] = mat.volume
# Make sure there are burnable materials
if not burnable_mats:
raise RuntimeError(
"No depletable materials were found in the model.")
# Sort the sets
burnable_mats = sorted(burnable_mats, key=int)
model_nuclides = sorted(model_nuclides)
# Construct a global nuclide dictionary, burned first
nuclides = list(self.chain.nuclide_dict)
for nuc in model_nuclides:
if nuc not in nuclides:
nuclides.append(nuc)
return burnable_mats, volume, nuclides
def _extract_number(self, local_mats, volume, nuclides):
"""Construct AtomNumber using geometry
Parameters
----------
local_mats : list of str
Material IDs to be managed by this process
volume : OrderedDict of str to float
Volumes for the above materials in [cm^3]
nuclides : list of str
Nuclides to be used in the simulation.
"""
self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain))
if self.dilute_initial != 0.0:
for nuc in self._burnable_nucs:
self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial)
# Now extract the number densities and store
for mat in self.geometry.get_all_materials().values():
if str(mat.id) in local_mats:
self._set_number_from_mat(mat)
def _set_number_from_mat(self, mat):
"""Extracts material and number densities from openmc.Material
Parameters
----------
mat : openmc.Material
The material to read from
"""
mat_id = str(mat.id)
for nuclide, density in mat.get_nuclide_atom_densities().values():
number = density * 1.0e24
self.number.set_atom_density(mat_id, nuclide, number)
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
# Create XML files
if comm.rank == 0:
self.geometry.export_to_xml()
self.settings.export_to_xml()
self._generate_materials_xml()
# Initialize OpenMC library
comm.barrier()
openmc.capi.init(comm)
# Generate tallies in memory
self._generate_tallies()
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
def finalize(self):
"""Finalize a depletion simulation and release resources."""
openmc.capi.finalize()
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
for rank in range(comm.size):
number_i = comm.bcast(self.number, root=rank)
for mat in number_i.materials:
nuclides = []
densities = []
for nuc in number_i.nuclides:
if nuc in self.nuclides_with_data:
val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
# If nuclide is zero, do not add to the problem.
if val > 0.0:
if self.round_number:
val_magnitude = np.floor(np.log10(val))
val_scaled = val / 10**val_magnitude
val_round = round(val_scaled, 8)
val = val_round * 10**val_magnitude
nuclides.append(nuc)
densities.append(val)
else:
# Only output warnings if values are significantly
# negative. CRAM does not guarantee positive values.
if val < -1.0e-21:
print("WARNING: nuclide ", nuc, " in material ", mat,
" is negative (density = ", val, " at/barn-cm)")
number_i[mat, nuc] = 0.0
mat_internal = openmc.capi.materials[int(mat)]
mat_internal.set_densities(nuclides, densities)
def _generate_materials_xml(self):
"""Creates materials.xml from self.number.
Due to uncertainty with how MPI interacts with OpenMC API, this
constructs the XML manually. The long term goal is to do this
through direct memory writing.
"""
materials = openmc.Materials(self.geometry.get_all_materials()
.values())
# Sort nuclides according to order in AtomNumber object
nuclides = list(self.number.nuclides)
for mat in materials:
mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
materials.export_to_xml()
def _get_tally_nuclides(self):
"""Determine nuclides that should be tallied for reaction rates.
This method returns a list of all nuclides that have neutron data and
are listed in the depletion chain. Technically, we should tally nuclides
that may not appear in the depletion chain because we still need to get
the fission reaction rate for these nuclides in order to normalize
power, but that is left as a future exercise.
Returns
-------
list of str
Tally nuclides
"""
nuc_set = set()
# Create the set of all nuclides in the decay chain in materials marked
# for burning in which the number density is greater than zero.
for nuc in self.number.nuclides:
if nuc in self.nuclides_with_data:
if np.sum(self.number[:, nuc]) > 0.0:
nuc_set.add(nuc)
# Communicate which nuclides have nonzeros to rank 0
if comm.rank == 0:
for i in range(1, comm.size):
nuc_newset = comm.recv(source=i, tag=i)
nuc_set |= nuc_newset
else:
comm.send(nuc_set, dest=0, tag=comm.rank)
if comm.rank == 0:
# Sort nuclides in the same order as self.number
nuc_list = [nuc for nuc in self.number.nuclides
if nuc in nuc_set]
else:
nuc_list = None
# Store list of tally nuclides on each process
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _generate_tallies(self):
"""Generates depletion tallies.
Using information from the depletion chain as well as the nuclides
currently in the problem, this function automatically generates a
tally.xml for the simulation.
"""
# Create tallies for depleting regions
materials = [openmc.capi.materials[int(i)]
for i in self.burnable_mats]
mat_filter = openmc.capi.MaterialFilter(materials)
# Set up a tally that has a material filter covering each depletable
# material and scores corresponding to all reactions that cause
# transmutation. The nuclides for the tally are set later when eval() is
# called.
self._tally = openmc.capi.Tally()
self._tally.scores = self.chain.reactions
self._tally.filters = [mat_filter]
def _unpack_tallies_and_normalize(self, power):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
value and reaction rates from the simulation. The reaction rates are
normalized by the user-specified power, summing the product of the
fission reaction rate times the fission Q value for each material.
Parameters
----------
power : float
Power of the reactor in [W]
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
"""
rates = self.reaction_rates
rates[:, :, :] = 0.0
k_combined = openmc.capi.keff()[0]
# Extract tally bins
materials = self.burnable_mats
nuclides = self._tally.nuclides
# Form fast map
nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Compute fission power
# TODO : improve this calculation
# Keep track of energy produced from all reactions in eV per source
# particle
energy = 0.0
# Create arrays to store fission Q values, reaction rates, and nuclide
# numbers
fission_Q = np.zeros(rates.n_nuc)
rates_expanded = np.zeros((rates.n_nuc, rates.n_react))
number = np.zeros(rates.n_nuc)
fission_ind = rates.index_rx["fission"]
for nuclide in self.chain.nuclides:
if nuclide.name in rates.index_nuc:
for rx in nuclide.reactions:
if rx.type == 'fission':
ind = rates.index_nuc[nuclide.name]
fission_Q[ind] = rx.Q
break
# Extract results
for i, mat in enumerate(self.local_mats):
# Get tally index
slab = materials.index(mat)
# Get material results hyperslab
results = self._tally.results[slab, :, 1]
# Zero out reaction rates and nuclide numbers
rates_expanded[:] = 0.0
number[:] = 0.0
# Expand into our memory layout
j = 0
for nuc, i_nuc_results in zip(nuclides, nuc_ind):
number[i_nuc_results] = self.number[mat, nuc]
for react in react_ind:
rates_expanded[i_nuc_results, react] = results[j]
j += 1
# Accumulate energy from fission
energy += np.dot(rates_expanded[:, fission_ind], fission_Q)
# Divide by total number and store
for i_nuc_results in nuc_ind:
if number[i_nuc_results] != 0.0:
for react in react_ind:
rates_expanded[i_nuc_results, react] /= number[i_nuc_results]
rates[i, :, :] = rates_expanded
# Reduce energy produced from all processes
energy = comm.allreduce(energy)
# Determine power in eV/s
power /= JOULE_PER_EV
# Scale reaction rates to obtain units of reactions/sec
rates *= power / energy
return OperatorResult(k_combined, rates)
def _get_nuclides_with_data(self):
"""Loads a cross_sections.xml file to find participating nuclides.
This allows for nuclides that are important in the decay chain but not
important neutronically, or have no cross section data.
"""
# Reads cross_sections.xml to create a dictionary containing
# participating (burning and not just decaying) nuclides.
try:
filename = os.environ["OPENMC_CROSS_SECTIONS"]
except KeyError:
filename = None
nuclides = set()
try:
tree = ET.parse(filename)
except Exception:
if filename is None:
msg = "No cross_sections.xml specified in materials."
else:
msg = 'Cross section file "{}" is invalid.'.format(filename)
raise IOError(msg)
root = tree.getroot()
for nuclide_node in root.findall('library'):
mats = nuclide_node.get('materials')
if not mats:
continue
for name in mats.split():
# Make a burn list of the union of nuclides in cross_sections.xml
# and nuclides in depletion chain.
if name not in nuclides:
nuclides.add(name)
return nuclides
def get_results_info(self):
"""Returns volume list, material lists, and nuc lists.
Returns
-------
volume : dict of str float
Volumes corresponding to materials in full_burn_dict
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all material IDs to be burned. Used for sorting the simulation.
full_burn_list : list
List of all burnable material IDs
"""
nuc_list = self.number.burnable_nuclides
burn_list = self.local_mats
volume = {}
for i, mat in enumerate(burn_list):
volume[mat] = self.number.volume[i]
# Combine volume dictionaries across processes
volume_list = comm.allgather(volume)
volume = {k: v for d in volume_list for k, v in d.items()}
return volume, nuc_list, burn_list, self.burnable_mats

139
openmc/deplete/reaction_rates.py Normal file
View file

@ -0,0 +1,139 @@
"""ReactionRates module.
An ndarray to store reaction rates with string, integer, or slice indexing.
"""
from collections import OrderedDict
import numpy as np
class ReactionRates(np.ndarray):
"""Reaction rates resulting from a transport operator call
This class is a subclass of :class:`numpy.ndarray` with a few custom
attributes that make it easy to determine what index corresponds to a given
material, nuclide, and reaction rate.
Parameters
----------
local_mats : list of str
Material IDs
nuclides : list of str
Depletable nuclides
reactions : list of str
Transmutation reactions being tracked
Attributes
----------
index_mat : OrderedDict of str to int
A dictionary mapping material ID as string to index.
index_nuc : OrderedDict of str to int
A dictionary mapping nuclide name as string to index.
index_rx : OrderedDict of str to int
A dictionary mapping reaction name as string to index.
n_mat : int
Number of materials.
n_nuc : int
Number of nucs.
n_react : int
Number of reactions.
"""
# NumPy arrays can be created 1) explicitly 2) using view casting, and 3) by
# slicing an existing array. Because of these possibilities, it's necessary
# to put initialization logic in __new__ rather than __init__. Additionally,
# subclasses need to handle the multiple ways of creating arrays by using
# the __array_finalize__ method (discussed here:
# https://docs.scipy.org/doc/numpy/user/basics.subclassing.html)
def __new__(cls, local_mats, nuclides, reactions):
# Create appropriately-sized zeroed-out ndarray
shape = (len(local_mats), len(nuclides), len(reactions))
obj = super().__new__(cls, shape)
obj[:] = 0.0
# Add mapping attributes
obj.index_mat = {mat: i for i, mat in enumerate(local_mats)}
obj.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)}
obj.index_rx = {rx: i for i, rx in enumerate(reactions)}
return obj
def __array_finalize__(self, obj):
if obj is None:
return
self.index_mat = getattr(obj, 'index_mat', None)
self.index_nuc = getattr(obj, 'index_nuc', None)
self.index_rx = getattr(obj, 'index_rx', None)
# Reaction rates are distributed to other processes via multiprocessing,
# which entails pickling the objects. In order to preserve the custom
# attributes, we have to modify how the ndarray is pickled as described
# here: https://stackoverflow.com/a/26599346/1572453
def __reduce__(self):
state = super().__reduce__()
new_state = state[2] + (self.index_mat, self.index_nuc, self.index_rx)
return (state[0], state[1], new_state)
def __setstate__(self, state):
self.index_mat = state[-3]
self.index_nuc = state[-2]
self.index_rx = state[-1]
super().__setstate__(state[0:-3])
@property
def n_mat(self):
return len(self.index_mat)
@property
def n_nuc(self):
return len(self.index_nuc)
@property
def n_react(self):
return len(self.index_rx)
def get(self, mat, nuc, rx):
"""Get reaction rate by material/nuclide/reaction
Parameters
----------
mat : str
Material ID as a string
nuc : str
Nuclide name
rx : str
Name of the reaction
Returns
-------
float
Reaction rate corresponding to given material, nuclide, and reaction
"""
mat = self.index_mat[mat]
nuc = self.index_nuc[nuc]
rx = self.index_rx[rx]
return self[mat, nuc, rx]
def set(self, mat, nuc, rx, value):
"""Set reaction rate by material/nuclide/reaction
Parameters
----------
mat : str
Material ID as a string
nuc : str
Nuclide name
rx : str
Name of the reaction
value : float
Corresponding reaction rate to set
"""
mat = self.index_mat[mat]
nuc = self.index_nuc[nuc]
rx = self.index_rx[rx]
self[mat, nuc, rx] = value

397
openmc/deplete/results.py Normal file
View file

@ -0,0 +1,397 @@
"""The results module.
Contains results generation and saving capabilities.
"""
from collections import OrderedDict
import copy
import numpy as np
import h5py
from . import comm, have_mpi
from .reaction_rates import ReactionRates
_VERSION_RESULTS = (1, 0)
class Results(object):
"""Output of a depletion run
Attributes
----------
k : list of float
Eigenvalue for each substep.
time : list of float
Time at beginning, end of step, in seconds.
n_mat : int
Number of mats.
n_nuc : int
Number of nuclides.
rates : list of ReactionRates
The reaction rates for each substep.
volume : OrderedDict of int to float
Dictionary mapping mat id to volume.
mat_to_ind : OrderedDict of str to int
A dictionary mapping mat ID as string to index.
nuc_to_ind : OrderedDict of str to int
A dictionary mapping nuclide name as string to index.
mat_to_hdf5_ind : OrderedDict of str to int
A dictionary mapping mat ID as string to global index.
n_hdf5_mats : int
Number of materials in entire geometry.
n_stages : int
Number of stages in simulation.
data : numpy.ndarray
Atom quantity, stored by stage, mat, then by nuclide.
"""
def __init__(self):
self.k = None
self.time = None
self.rates = None
self.volume = None
self.mat_to_ind = None
self.nuc_to_ind = None
self.mat_to_hdf5_ind = None
self.data = None
def __getitem__(self, pos):
"""Retrieves an item from results.
Parameters
----------
pos : tuple
A three-length tuple containing a stage index, mat index and a nuc
index. All can be integers or slices. The second two can be
strings corresponding to their respective dictionary.
Returns
-------
float
The atoms for stage, mat, nuc
"""
stage, mat, nuc = pos
if isinstance(mat, str):
mat = self.mat_to_ind[mat]
if isinstance(nuc, str):
nuc = self.nuc_to_ind[nuc]
return self.data[stage, mat, nuc]
def __setitem__(self, pos, val):
"""Sets an item from results.
Parameters
----------
pos : tuple
A three-length tuple containing a stage index, mat index and a nuc
index. All can be integers or slices. The second two can be
strings corresponding to their respective dictionary.
val : float
The value to set data to.
"""
stage, mat, nuc = pos
if isinstance(mat, str):
mat = self.mat_to_ind[mat]
if isinstance(nuc, str):
nuc = self.nuc_to_ind[nuc]
self.data[stage, mat, nuc] = val
@property
def n_mat(self):
return len(self.mat_to_ind)
@property
def n_nuc(self):
return len(self.nuc_to_ind)
@property
def n_hdf5_mats(self):
return len(self.mat_to_hdf5_ind)
@property
def n_stages(self):
return self.data.shape[0]
def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
"""Allocates memory of Results.
Parameters
----------
volume : dict of str float
Volumes corresponding to materials in full_burn_dict
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all mat IDs to be burned. Used for sorting the simulation.
full_burn_list : list of str
List of all burnable material IDs
stages : int
Number of stages in simulation.
"""
self.volume = copy.deepcopy(volume)
self.nuc_to_ind = {nuc: i for i, nuc in enumerate(nuc_list)}
self.mat_to_ind = {mat: i for i, mat in enumerate(burn_list)}
self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
# Create storage array
self.data = np.zeros((stages, self.n_mat, self.n_nuc))
def export_to_hdf5(self, filename, step):
"""Export results to an HDF5 file
Parameters
----------
filename : str
The filename to write to
step : int
What step is this?
"""
if have_mpi and h5py.get_config().mpi:
kwargs = {'driver': 'mpio', 'comm': comm}
else:
kwargs = {}
kwargs['mode'] = "w" if step == 0 else "a"
with h5py.File(filename, **kwargs) as handle:
self._to_hdf5(handle, step)
def _write_hdf5_metadata(self, handle):
"""Writes result metadata in HDF5 file
Parameters
----------
handle : h5py.File or h5py.Group
An hdf5 file or group type to store this in.
"""
# Create and save the 5 dictionaries:
# quantities
# self.mat_to_ind -> self.volume (TODO: support for changing volumes)
# self.nuc_to_ind
# reactions
# self.rates[0].nuc_to_ind (can be different from above, above is superset)
# self.rates[0].react_to_ind
# these are shared by every step of the simulation, and should be deduplicated.
# Store concentration mat and nuclide dictionaries (along with volumes)
handle.attrs['version'] = np.array(_VERSION_RESULTS)
handle.attrs['filetype'] = np.string_('depletion results')
mat_list = sorted(self.mat_to_hdf5_ind, key=int)
nuc_list = sorted(self.nuc_to_ind)
rxn_list = sorted(self.rates[0].index_rx)
n_mats = self.n_hdf5_mats
n_nuc_number = len(nuc_list)
n_nuc_rxn = len(self.rates[0].index_nuc)
n_rxn = len(rxn_list)
n_stages = self.n_stages
mat_group = handle.create_group("materials")
for mat in mat_list:
mat_single_group = mat_group.create_group(mat)
mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat]
mat_single_group.attrs["volume"] = self.volume[mat]
nuc_group = handle.create_group("nuclides")
for nuc in nuc_list:
nuc_single_group = nuc_group.create_group(nuc)
nuc_single_group.attrs["atom number index"] = self.nuc_to_ind[nuc]
if nuc in self.rates[0].index_nuc:
nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc]
rxn_group = handle.create_group("reactions")
for rxn in rxn_list:
rxn_single_group = rxn_group.create_group(rxn)
rxn_single_group.attrs["index"] = self.rates[0].index_rx[rxn]
# Construct array storage
handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number),
maxshape=(None, n_stages, n_mats, n_nuc_number),
chunks=(1, 1, n_mats, n_nuc_number),
dtype='float64')
handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
dtype='float64')
handle.create_dataset("eigenvalues", (1, n_stages),
maxshape=(None, n_stages), dtype='float64')
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
def _to_hdf5(self, handle, index):
"""Converts results object into an hdf5 object.
Parameters
----------
handle : h5py.File or h5py.Group
An HDF5 file or group type to store this in.
index : int
What step is this?
"""
if "/number" not in handle:
comm.barrier()
self._write_hdf5_metadata(handle)
comm.barrier()
# Grab handles
number_dset = handle["/number"]
rxn_dset = handle["/reaction rates"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
# Get number of results stored
number_shape = list(number_dset.shape)
number_results = number_shape[0]
new_shape = index + 1
if number_results < new_shape:
# Extend first dimension by 1
number_shape[0] = new_shape
number_dset.resize(number_shape)
rxn_shape = list(rxn_dset.shape)
rxn_shape[0] = new_shape
rxn_dset.resize(rxn_shape)
eigenvalues_shape = list(eigenvalues_dset.shape)
eigenvalues_shape[0] = new_shape
eigenvalues_dset.resize(eigenvalues_shape)
time_shape = list(time_dset.shape)
time_shape[0] = new_shape
time_dset.resize(time_shape)
# If nothing to write, just return
if len(self.mat_to_ind) == 0:
return
# Add data
# Note, for the last step, self.n_stages = 1, even if n_stages != 1.
n_stages = self.n_stages
inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind]
low = min(inds)
high = max(inds)
for i in range(n_stages):
number_dset[index, i, low:high+1, :] = self.data[i, :, :]
rxn_dset[index, i, low:high+1, :, :] = self.rates[i][:, :, :]
if comm.rank == 0:
eigenvalues_dset[index, i] = self.k[i]
if comm.rank == 0:
time_dset[index, :] = self.time
@classmethod
def from_hdf5(cls, handle, step):
"""Loads results object from HDF5.
Parameters
----------
handle : h5py.File or h5py.Group
An HDF5 file or group type to load from.
step : int
What step is this?
"""
results = cls()
# Grab handles
number_dset = handle["/number"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
results.data = number_dset[step, :, :, :]
results.k = eigenvalues_dset[step, :]
results.time = time_dset[step, :]
# Reconstruct dictionaries
results.volume = OrderedDict()
results.mat_to_ind = OrderedDict()
results.nuc_to_ind = OrderedDict()
rxn_nuc_to_ind = OrderedDict()
rxn_to_ind = OrderedDict()
for mat, mat_handle in handle["/materials"].items():
vol = mat_handle.attrs["volume"]
ind = mat_handle.attrs["index"]
results.volume[mat] = vol
results.mat_to_ind[mat] = ind
for nuc, nuc_handle in handle["/nuclides"].items():
ind_atom = nuc_handle.attrs["atom number index"]
results.nuc_to_ind[nuc] = ind_atom
if "reaction rate index" in nuc_handle.attrs:
rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"]
for rxn, rxn_handle in handle["/reactions"].items():
rxn_to_ind[rxn] = rxn_handle.attrs["index"]
results.rates = []
# Reconstruct reactions
for i in range(results.n_stages):
rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind)
rate[:] = handle["/reaction rates"][step, i, :, :, :]
results.rates.append(rate)
return results
@staticmethod
def save(op, x, op_results, t, step_ind):
"""Creates and writes depletion results to disk
Parameters
----------
op : openmc.deplete.TransportOperator
The operator used to generate these results.
x : list of list of numpy.array
The prior x vectors. Indexed [i][cell] using the above equation.
op_results : list of openmc.deplete.OperatorResult
Results of applying transport operator
t : list of float
Time indices.
step_ind : int
Step index.
"""
# Get indexing terms
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
# Create results
stages = len(x)
results = Results()
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
n_mat = len(burn_list)
for i in range(stages):
for mat_i in range(n_mat):
results[i, mat_i, :] = x[i][mat_i][:]
results.k = [r.k for r in op_results]
results.rates = [r.rates for r in op_results]
results.time = t
results.export_to_hdf5("depletion_results.h5", step_ind)

105
openmc/deplete/results_list.py Normal file
View file

@ -0,0 +1,105 @@
import h5py
import numpy as np
from .results import Results, _VERSION_RESULTS
from openmc.checkvalue import check_filetype_version
class ResultsList(list):
"""A list of openmc.deplete.Results objects
Parameters
----------
filename : str
The filename to read from.
"""
def __init__(self, filename):
super().__init__()
with h5py.File(str(filename), "r") as fh:
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
# Get number of results stored
n = fh["number"].value.shape[0]
for i in range(n):
self.append(Results.from_hdf5(fh, i))
def get_atoms(self, mat, nuc):
"""Get nuclide concentration over time from a single material
Parameters
----------
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
Returns
-------
time : numpy.ndarray
Array of times in [s]
concentration : numpy.ndarray
Total number of atoms for specified nuclide
"""
time = np.empty_like(self)
concentration = np.empty_like(self)
# Evaluate value in each region
for i, result in enumerate(self):
time[i] = result.time[0]
concentration[i] = result[0, mat, nuc]
return time, concentration
def get_reaction_rate(self, mat, nuc, rx):
"""Get reaction rate in a single material/nuclide over time
Parameters
----------
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
rx : str
Reaction rate to evaluate
Returns
-------
time : numpy.ndarray
Array of times in [s]
rate : numpy.ndarray
Array of reaction rates
"""
time = np.empty_like(self)
rate = np.empty_like(self)
# Evaluate value in each region
for i, result in enumerate(self):
time[i] = result.time[0]
rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
return time, rate
def get_eigenvalue(self):
"""Evaluates the eigenvalue from a results list.
Returns
-------
time : numpy.ndarray
Array of times in [s]
eigenvalue : numpy.ndarray
k-eigenvalue at each time
"""
time = np.empty_like(self)
eigenvalue = np.empty_like(self)
# Get time/eigenvalue at each point
for i, result in enumerate(self):
time[i] = result.time[0]
eigenvalue[i] = result.k[0]
return time, eigenvalue

6
openmc/element.py
View file

@ -1,8 +1,6 @@
from collections import OrderedDict
import re
import os
from six import string_types
from xml.etree import ElementTree as ET
import openmc
@ -29,9 +27,9 @@ class Element(str):
"""
def __new__(cls, name):
cv.check_type('element name', name, string_types)
cv.check_type('element name', name, str)
cv.check_length('element name', name, 1, 2)
return super(Element, cls).__new__(cls, name)
return super().__new__(cls, name)
@property
def name(self):

2
openmc/examples.py
View file

@ -587,7 +587,7 @@ def slab_mg(reps=None, as_macro=True):
# Define the materials file
model.xs_data = xs
model.materials.cross_sections = "../1d_mgxs.h5"
model.materials.cross_sections = "../../1d_mgxs.h5"
# Define surfaces.
# Assembly/Problem Boundary

42
openmc/executor.py
View file

@ -1,10 +1,7 @@
from __future__ import print_function
from collections import Iterable
from collections.abc import Iterable
import subprocess
from numbers import Integral
from six import string_types
import openmc
from openmc import VolumeCalculation
@ -15,18 +12,22 @@ def _run(args, output, cwd):
stderr=subprocess.STDOUT, universal_newlines=True)
# Capture and re-print OpenMC output in real-time
lines = []
while True:
# If OpenMC is finished, break loop
line = p.stdout.readline()
if not line and p.poll() is not None:
break
lines.append(line)
if output:
# If user requested output, print to screen
print(line, end='')
# Return the returncode (integer, zero if no problems encountered)
return p.returncode
# Raise an exception if return status is non-zero
if p.returncode != 0:
raise subprocess.CalledProcessError(p.returncode, ' '.join(args),
''.join(lines))
def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
@ -41,8 +42,13 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
cwd : str, optional
Path to working directory to run in
Raises
------
subprocess.CalledProcessError
If the `openmc` executable returns a non-zero status
"""
return _run([openmc_exec, '-p'], output, cwd)
_run([openmc_exec, '-p'], output, cwd)
def plot_inline(plots, openmc_exec='openmc', cwd='.', convert_exec='convert'):
@ -63,6 +69,11 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.', convert_exec='convert'):
convert_exec : str, optional
Command that can convert PPM files into PNG files
Raises
------
subprocess.CalledProcessError
If the `openmc` executable returns a non-zero status
"""
from IPython.display import Image, display
@ -121,6 +132,11 @@ def calculate_volumes(threads=None, output=True, cwd='.',
Path to working directory to run in. Defaults to the current working
directory.
Raises
------
subprocess.CalledProcessError
If the `openmc` executable returns a non-zero status
See Also
--------
openmc.VolumeCalculation
@ -133,7 +149,7 @@ def calculate_volumes(threads=None, output=True, cwd='.',
if mpi_args is not None:
args = mpi_args + args
return _run(args, output, cwd)
_run(args, output, cwd)
def run(particles=None, threads=None, geometry_debug=False,
@ -167,8 +183,12 @@ def run(particles=None, threads=None, geometry_debug=False,
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
"""
Raises
------
subprocess.CalledProcessError
If the `openmc` executable returns a non-zero status
"""
args = [openmc_exec]
if isinstance(particles, Integral) and particles > 0:
@ -180,7 +200,7 @@ def run(particles=None, threads=None, geometry_debug=False,
if geometry_debug:
args.append('-g')
if isinstance(restart_file, string_types):
if isinstance(restart_file, str):
args += ['-r', restart_file]
if tracks:
@ -189,4 +209,4 @@ def run(particles=None, threads=None, geometry_debug=False,
if mpi_args is not None:
args = mpi_args + args
return _run(args, output, cwd)
_run(args, output, cwd)

296
openmc/filter.py
View file

@ -1,4 +1,3 @@
from __future__ import division
from abc import ABCMeta
from collections import Iterable, OrderedDict
import copy
@ -8,7 +7,6 @@ from numbers import Real, Integral
import operator
from xml.etree import ElementTree as ET
from six import add_metaclass
import numpy as np
import pandas as pd
@ -17,6 +15,7 @@ import openmc.checkvalue as cv
from .cell import Cell
from .material import Material
from .mixin import IDManagerMixin
from .surface import Surface
from .universe import Universe
@ -24,12 +23,11 @@ _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface',
'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal',
'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom']
_CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in',
3: 'x-max out', 4: 'x-max in',
5: 'y-min out', 6: 'y-min in',
7: 'y-max out', 8: 'y-max in',
9: 'z-min out', 10: 'z-min in',
11: 'z-max out', 12: 'z-max in'}
_CURRENT_NAMES = (
'x-min out', 'x-min in', 'x-max out', 'x-max in',
'y-min out', 'y-min in', 'y-max out', 'y-max in',
'z-min out', 'z-min in', 'z-max out', 'z-max in'
)
class FilterMeta(ABCMeta):
@ -66,12 +64,10 @@ class FilterMeta(ABCMeta):
namespace[func_name].__doc__ = old_doc
# Make the class.
return super(FilterMeta, cls).__new__(cls, name, bases, namespace,
**kwargs)
return super().__new__(cls, name, bases, namespace, **kwargs)
@add_metaclass(FilterMeta)
class Filter(IDManagerMixin):
class Filter(IDManagerMixin, metaclass=FilterMeta):
"""Tally modifier that describes phase-space and other characteristics.
Parameters
@ -501,9 +497,9 @@ class CellFilter(WithIDFilter):
Parameters
----------
bins : openmc.Cell, Integral, or iterable thereof
The Cells to tally. Either openmc.Cell objects or their
Integral ID numbers can be used.
bins : openmc.Cell, int, or iterable thereof
The cells to tally. Either openmc.Cell objects or their ID numbers can
be used.
filter_id : int
Unique identifier for the filter
@ -568,86 +564,29 @@ class CellbornFilter(WithIDFilter):
expected_type = Cell
class SurfaceFilter(Filter):
"""Bins particle currents on Mesh surfaces.
class SurfaceFilter(WithIDFilter):
"""Filters particles by surface crossing
Parameters
----------
bins : Iterable of Integral
Indices corresponding to which face of a mesh cell the current is
crossing.
bins : openmc.Surface, int, or iterable of Integral
The surfaces to tally over. Either openmc.Surface objects or their ID
numbers can be used.
filter_id : int
Unique identifier for the filter
Attributes
----------
bins : Iterable of Integral
Indices corresponding to which face of a mesh cell the current is
crossing.
The surfaces to tally over. Either openmc.Surface objects or their ID
numbers can be used.
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
"""
@Filter.bins.setter
def bins(self, bins):
# Format the bins as a 1D numpy array.
bins = np.atleast_1d(bins)
# Check the bin values.
cv.check_iterable_type('filter bins', bins, Integral)
for edge in bins:
cv.check_greater_than('filter bin', edge, 0, equality=True)
self._bins = bins
@property
def num_bins(self):
# Need to handle number of bins carefully -- for surface current
# tallies, the number of bins depends on the mesh, which we don't have a
# reference to in this filter
return self._num_bins
def get_pandas_dataframe(self, data_size, stride, **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 : int
The total number of bins in the tally corresponding to this filter
stride : int
Stride in memory for the filter
Returns
-------
pandas.DataFrame
A Pandas DataFrame with a column of strings describing which surface
the current is crossing and which direction it points. 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.
See also
--------
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
df = pd.DataFrame()
filter_bins = np.repeat(self.bins, stride)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_bins = [_CURRENT_NAMES[x] for x in filter_bins]
df = pd.concat([df, pd.DataFrame(
{self.short_name.lower(): filter_bins})])
return df
expected_type = Surface
class MeshFilter(Filter):
@ -675,7 +614,7 @@ class MeshFilter(Filter):
def __init__(self, mesh, filter_id=None):
self.mesh = mesh
super(MeshFilter, self).__init__(mesh.id, filter_id)
super().__init__(mesh.id, filter_id)
@classmethod
def from_hdf5(cls, group, **kwargs):
@ -693,7 +632,6 @@ class MeshFilter(Filter):
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(mesh_obj, filter_id=filter_id)
out._num_bins = group['n_bins'].value
return out
@ -709,10 +647,7 @@ class MeshFilter(Filter):
@property
def num_bins(self):
try:
return self._num_bins
except AttributeError:
return reduce(operator.mul, self.mesh.dimension)
return reduce(operator.mul, self.mesh.dimension)
def check_bins(self, bins):
if not len(bins) == 1:
@ -736,37 +671,40 @@ class MeshFilter(Filter):
# Filter bins for a mesh are an (x,y,z) tuple. Convert (x,y,z) to a
# single bin -- this is similar to subroutine mesh_indices_to_bin in
# openmc/src/mesh.F90.
if len(self.mesh.dimension) == 3:
n_dim = len(self.mesh.dimension)
if n_dim == 3:
i, j, k = filter_bin
nx, ny, nz = self.mesh.dimension
val = (filter_bin[0] - 1) * ny * nz + \
(filter_bin[1] - 1) * nz + \
(filter_bin[2] - 1)
else:
return (i - 1) + (j - 1)*nx + (k - 1)*nx*ny
elif n_dim == 2:
i, j, *_ = filter_bin
nx, ny = self.mesh.dimension
val = (filter_bin[0] - 1) * ny + \
(filter_bin[1] - 1)
return val
return (i - 1) + (j - 1)*nx
else:
return filter_bin[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 3-tuple of x,y,z cell indices for a 3D mesh
if len(self.mesh.dimension) == 3:
n_dim = len(self.mesh.dimension)
if n_dim == 3:
# Construct 3-tuple of x,y,z cell indices for a 3D mesh
nx, ny, nz = self.mesh.dimension
x = bin_index / (ny * nz)
y = (bin_index - (x * ny * nz)) / nz
z = bin_index - (x * ny * nz) - (y * nz)
x = (bin_index % nx) + 1
y = (bin_index % (nx * ny)) // nx + 1
z = bin_index // (nx * ny) + 1
return (x, y, z)
# Construct 2-tuple of x,y cell indices for a 2D mesh
else:
elif n_dim == 2:
# Construct 2-tuple of x,y cell indices for a 2D mesh
nx, ny = self.mesh.dimension
x = bin_index / ny
y = bin_index - (x * ny)
x = (bin_index % nx) + 1
y = bin_index // nx + 1
return (x, y)
else:
return (bin_index + 1,)
def get_pandas_dataframe(self, data_size, stride, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
@ -803,7 +741,134 @@ class MeshFilter(Filter):
filter_dict = {}
# Append Mesh ID as outermost index of multi-index
mesh_key = 'mesh {0}'.format(self.mesh.id)
mesh_key = 'mesh {}'.format(self.mesh.id)
# Find mesh dimensions - use 3D indices for simplicity
n_dim = len(self.mesh.dimension)
if n_dim == 3:
nx, ny, nz = self.mesh.dimension
elif n_dim == 2:
nx, ny = self.mesh.dimension
nz = 1
else:
nx = self.mesh.dimension
ny = nz = 1
# Generate multi-index sub-column for x-axis
filter_bins = np.arange(1, nx + 1)
repeat_factor = stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'x')] = filter_bins
# Generate multi-index sub-column for y-axis
filter_bins = np.arange(1, ny + 1)
repeat_factor = nx * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'y')] = filter_bins
# Generate multi-index sub-column for z-axis
filter_bins = np.arange(1, nz + 1)
repeat_factor = nx * ny * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'z')] = filter_bins
# Initialize a Pandas DataFrame from the mesh dictionary
df = pd.concat([df, pd.DataFrame(filter_dict)])
return df
class MeshSurfaceFilter(MeshFilter):
"""Filter events by surface crossings on a regular, rectangular mesh.
Parameters
----------
mesh : openmc.Mesh
The Mesh object that events will be tallied onto
filter_id : int
Unique identifier for the filter
Attributes
----------
bins : Integral
The Mesh ID
mesh : openmc.Mesh
The Mesh object that events will be tallied onto
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
"""
@property
def num_bins(self):
n_dim = len(self.mesh.dimension)
return 4*n_dim*reduce(operator.mul, self.mesh.dimension)
def get_bin_index(self, filter_bin):
# Split bin into mesh/surface parts
*mesh_tuple, surf = filter_bin
# Get index for mesh alone
mesh_index = super().get_bin_index(mesh_tuple)
# Combine surface and mesh index
n_dim = len(self.mesh.dimension)
for surf_index, name in enumerate(_CURRENT_NAMES):
if surf == name:
return 4*n_dim*mesh_index + surf_index
else:
raise ValueError("'{}' is not a valid mesh surface.".format(surf))
def get_bin(self, bin_index):
n_dim = len(self.mesh.dimension)
mesh_index, surf_index = divmod(bin_index, 4*n_dim)
mesh_bin = super().get_bin(mesh_index)
return mesh_bin + (_CURRENT_NAMES[surf_index],)
def get_pandas_dataframe(self, data_size, stride, **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 : int
The total number of bins in the tally corresponding to this filter
stride : int
Stride in memory for the filter
Returns
-------
pandas.DataFrame
A Pandas DataFrame with three columns describing the x,y,z mesh
cell indices corresponding to 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.
See also
--------
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
df = pd.DataFrame()
# Initialize dictionary to build Pandas Multi-index column
filter_dict = {}
# Append Mesh ID as outermost index of multi-index
mesh_key = 'mesh {}'.format(self.mesh.id)
# Find mesh dimensions - use 3D indices for simplicity
if len(self.mesh.dimension) == 3:
@ -817,7 +882,7 @@ class MeshFilter(Filter):
# Generate multi-index sub-column for x-axis
filter_bins = np.arange(1, nx + 1)
repeat_factor = ny * nz * stride
repeat_factor = 12 * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
@ -825,7 +890,7 @@ class MeshFilter(Filter):
# Generate multi-index sub-column for y-axis
filter_bins = np.arange(1, ny + 1)
repeat_factor = nz * stride
repeat_factor = 12 * nx * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
@ -833,16 +898,21 @@ class MeshFilter(Filter):
# Generate multi-index sub-column for z-axis
filter_bins = np.arange(1, nz + 1)
repeat_factor = stride
repeat_factor = 12 * nx * ny * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'z')] = filter_bins
# Initialize a Pandas DataFrame from the mesh dictionary
df = pd.concat([df, pd.DataFrame(filter_dict)])
# Generate multi-index sub-column for surface
repeat_factor = stride
filter_bins = np.repeat(_CURRENT_NAMES, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'surf')] = filter_bins
return df
# Initialize a Pandas DataFrame from the mesh dictionary
return pd.concat([df, pd.DataFrame(filter_dict)])
class RealFilter(Filter):
@ -872,7 +942,7 @@ class RealFilter(Filter):
# This logic is used when merging tallies with real filters
return self.bins[0] >= other.bins[-1]
else:
return super(RealFilter, self).__gt__(other)
return super().__gt__(other)
@property
def num_bins(self):
@ -1130,7 +1200,7 @@ class DistribcellFilter(Filter):
def __init__(self, cell, filter_id=None):
self._paths = None
super(DistribcellFilter, self).__init__(cell, filter_id)
super().__init__(cell, filter_id)
@classmethod
def from_hdf5(cls, group, **kwargs):

46
openmc/geometry.py
View file

@ -1,9 +1,8 @@
from collections import OrderedDict, Iterable
from collections import OrderedDict
from collections.abc import Iterable
from copy import deepcopy
from xml.etree import ElementTree as ET
from six import string_types
import openmc
from openmc.clean_xml import clean_xml_indentation
from openmc.checkvalue import check_type
@ -14,8 +13,9 @@ class Geometry(object):
Parameters
----------
root_universe : openmc.Universe, optional
Root universe which contains all others
root : openmc.Universe or Iterable of openmc.Cell, optional
Root universe which contains all others, or an iterable of cells that
should be used to create a root universe.
Attributes
----------
@ -27,11 +27,17 @@ class Geometry(object):
"""
def __init__(self, root_universe=None):
def __init__(self, root=None):
self._root_universe = None
self._offsets = {}
if root_universe is not None:
self.root_universe = root_universe
if root is not None:
if isinstance(root, openmc.Universe):
self.root_universe = root
else:
univ = openmc.Universe()
for cell in root:
univ.add_cell(cell)
self._root_universe = univ
@property
def root_universe(self):
@ -132,7 +138,7 @@ class Geometry(object):
"""
# Make sure we are working with an iterable
return_list = (isinstance(paths, Iterable) and
not isinstance(paths, string_types))
not isinstance(paths, str))
path_list = paths if return_list else [paths]
indices = []
@ -249,7 +255,7 @@ class Geometry(object):
for cell in self.get_all_cells().values():
if cell.fill_type == 'lattice':
if cell.fill not in lattices:
if cell.fill.id not in lattices:
lattices[cell.fill.id] = cell.fill
return lattices
@ -306,9 +312,7 @@ class Geometry(object):
elif not matching and name in material_name:
materials.add(material)
materials = list(materials)
materials.sort(key=lambda x: x.id)
return materials
return sorted(materials, key=lambda x: x.id)
def get_cells_by_name(self, name, case_sensitive=False, matching=False):
"""Return a list of cells with matching names.
@ -346,9 +350,7 @@ class Geometry(object):
elif not matching and name in cell_name:
cells.add(cell)
cells = list(cells)
cells.sort(key=lambda x: x.id)
return cells
return sorted(cells, key=lambda x: x.id)
def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False):
"""Return a list of cells with fills with matching names.
@ -393,9 +395,7 @@ class Geometry(object):
elif not matching and name in fill_name:
cells.add(cell)
cells = list(cells)
cells.sort(key=lambda x: x.id)
return cells
return sorted(cells, key=lambda x: x.id)
def get_universes_by_name(self, name, case_sensitive=False, matching=False):
"""Return a list of universes with matching names.
@ -433,9 +433,7 @@ class Geometry(object):
elif not matching and name in universe_name:
universes.add(universe)
universes = list(universes)
universes.sort(key=lambda x: x.id)
return universes
return sorted(universes, key=lambda x: x.id)
def get_lattices_by_name(self, name, case_sensitive=False, matching=False):
"""Return a list of lattices with matching names.
@ -473,9 +471,7 @@ class Geometry(object):
elif not matching and name in lattice_name:
lattices.add(lattice)
lattices = list(lattices)
lattices.sort(key=lambda x: x.id)
return lattices
return sorted(lattices, key=lambda x: x.id)
def determine_paths(self, instances_only=False):
"""Determine paths through CSG tree for cells and materials.

40
openmc/lattice.py
View file

@ -1,13 +1,11 @@
from __future__ import division
from abc import ABCMeta
from collections import OrderedDict, Iterable
from collections import OrderedDict
from collections.abc import Iterable
from copy import deepcopy
from math import sqrt, floor
from numbers import Real, Integral
from xml.etree import ElementTree as ET
from six import add_metaclass, string_types
import numpy as np
import openmc.checkvalue as cv
@ -15,8 +13,7 @@ import openmc
from openmc.mixin import IDManagerMixin
@add_metaclass(ABCMeta)
class Lattice(IDManagerMixin):
class Lattice(IDManagerMixin, metaclass=ABCMeta):
"""A repeating structure wherein each element is a universe.
Parameters
@ -54,25 +51,6 @@ class Lattice(IDManagerMixin):
self._outer = None
self._universes = None
def __eq__(self, other):
if not isinstance(other, Lattice):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
elif np.any(self.pitch != other.pitch):
return False
elif self.outer != other.outer:
return False
elif np.any(self.universes != other.universes):
return False
else:
return True
def __ne__(self, other):
return not self == other
@property
def name(self):
return self._name
@ -92,7 +70,7 @@ class Lattice(IDManagerMixin):
@name.setter
def name(self, name):
if name is not None:
cv.check_type('lattice name', name, string_types)
cv.check_type('lattice name', name, str)
self._name = name
else:
self._name = ''
@ -512,7 +490,7 @@ class RectLattice(Lattice):
"""
def __init__(self, lattice_id=None, name=''):
super(RectLattice, self).__init__(lattice_id, name)
super().__init__(lattice_id, name)
# Initialize Lattice class attributes
self._lower_left = None
@ -579,7 +557,11 @@ class RectLattice(Lattice):
@property
def ndim(self):
return len(self.pitch)
if self.pitch is not None:
return len(self.pitch)
else:
raise ValueError('Number of dimensions cannot be determined until '
'the lattice pitch has been set.')
@property
def shape(self):
@ -839,7 +821,7 @@ class HexLattice(Lattice):
"""
def __init__(self, lattice_id=None, name=''):
super(HexLattice, self).__init__(lattice_id, name)
super().__init__(lattice_id, name)
# Initialize Lattice class attributes
self._num_rings = None

6
openmc/macroscopic.py
View file

@ -1,5 +1,3 @@
from six import string_types
from openmc.checkvalue import check_type
@ -19,8 +17,8 @@ class Macroscopic(str):
"""
def __new__(cls, name):
check_type('name', name, string_types)
return super(Macroscopic, cls).__new__(cls, name)
check_type('name', name, str)
return super().__new__(cls, name)
@property
def name(self):

207
openmc/material.py
View file

@ -4,7 +4,6 @@ from numbers import Real, Integral
import warnings
from xml.etree import ElementTree as ET
from six import string_types
import numpy as np
import openmc
@ -15,7 +14,7 @@ from .mixin import IDManagerMixin
# Units for density supported by OpenMC
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum',
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
'macro']
@ -25,7 +24,7 @@ class Material(IDManagerMixin):
To create a material, one should create an instance of this class, add
nuclides or elements with :meth:`Material.add_nuclide` or
`Material.add_element`, respectively, and set the total material density
with `Material.export_to_xml()`. The material can then be assigned to a cell
with `Material.set_density()`. The material can then be assigned to a cell
using the :attr:`Cell.fill` attribute.
Parameters
@ -49,12 +48,11 @@ class Material(IDManagerMixin):
density : float
Density of the material (units defined separately)
density_units : str
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/cm3',
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/m3',
'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only
applies in the case of a multi-group calculation.
depletable : bool
Indicate whether the material is depletable. This attribute can be used
by downstream depletion applications.
Indicate whether the material is depletable.
nuclides : list of tuple
List in which each item is a 3-tuple consisting of a nuclide string, the
percent density, and the percent type ('ao' or 'wo').
@ -75,6 +73,9 @@ class Material(IDManagerMixin):
:meth:`Geometry.determine_paths` method.
num_instances : int
The number of instances of this material throughout the geometry.
fissionable_mass : float
Mass of fissionable nuclides in the material in [g]. Requires that the
:attr:`volume` attribute is set.
"""
@ -87,7 +88,7 @@ class Material(IDManagerMixin):
self.name = name
self.temperature = temperature
self._density = None
self._density_units = ''
self._density_units = 'sum'
self._depletable = False
self._paths = None
self._num_instances = None
@ -217,7 +218,7 @@ class Material(IDManagerMixin):
def name(self, name):
if name is not None:
cv.check_type('name for Material ID="{}"'.format(self._id),
name, string_types)
name, str)
self._name = name
else:
self._name = ''
@ -243,9 +244,21 @@ class Material(IDManagerMixin):
@isotropic.setter
def isotropic(self, isotropic):
cv.check_iterable_type('Isotropic scattering nuclides', isotropic,
string_types)
str)
self._isotropic = list(isotropic)
@property
def fissionable_mass(self):
if self.volume is None:
raise ValueError("Volume must be set in order to determine mass.")
density = 0.0
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
Z = openmc.data.zam(nuc)[0]
if Z >= 90:
density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
/ openmc.data.AVOGADRO
return density*self.volume
@classmethod
def from_hdf5(cls, group):
"""Create material from HDF5 group
@ -300,7 +313,7 @@ class Material(IDManagerMixin):
"""
if volume_calc.domain_type == 'material':
if self.id in volume_calc.volumes:
self._volume = volume_calc.volumes[self.id][0]
self._volume = volume_calc.volumes[self.id].n
self._atoms = volume_calc.atoms[self.id]
else:
raise ValueError('No volume information found for this material.')
@ -312,7 +325,7 @@ class Material(IDManagerMixin):
Parameters
----------
units : {'g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro'}
units : {'g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro'}
Physical units of density.
density : float, optional
Value of the density. Must be specified unless units is given as
@ -345,7 +358,7 @@ class Material(IDManagerMixin):
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
if not isinstance(filename, string_types) and filename is not None:
if not isinstance(filename, str) and filename is not None:
msg = 'Unable to add OTF material file to Material ID="{}" with a ' \
'non-string name "{}"'.format(self._id, filename)
raise ValueError(msg)
@ -366,33 +379,31 @@ class Material(IDManagerMixin):
Parameters
----------
nuclide : str
Nuclide to add
Nuclide to add, e.g., 'Mo95'
percent : float
Atom or weight percent
percent_type : {'ao', 'wo'}
'ao' for atom percent and 'wo' for weight percent
"""
cv.check_type('nuclide', nuclide, str)
cv.check_type('percent', percent, Real)
cv.check_value('percent type', percent_type, {'ao', 'wo'})
if self._macroscopic is not None:
msg = 'Unable to add a Nuclide to Material ID="{}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(nuclide, string_types):
msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \
'non-string value "{}"'.format(self._id, nuclide)
raise ValueError(msg)
elif not isinstance(percent, Real):
msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \
'non-floating point value "{}"'.format(self._id, percent)
raise ValueError(msg)
elif percent_type not in ('ao', 'wo'):
msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \
'percent type "{}"'.format(self._id, percent_type)
raise ValueError(msg)
# If nuclide name doesn't look valid, give a warning
try:
Z, _, _ = openmc.data.zam(nuclide)
except ValueError as e:
warnings.warn(str(e))
else:
# For actinides, have the material be depletable by default
if Z >= 89:
self.depletable = True
self._nuclides.append((nuclide, percent, percent_type))
@ -405,7 +416,7 @@ class Material(IDManagerMixin):
Nuclide to remove
"""
cv.check_type('nuclide', nuclide, string_types)
cv.check_type('nuclide', nuclide, str)
# If the Material contains the Nuclide, delete it
for nuc in self._nuclides:
@ -434,7 +445,7 @@ class Material(IDManagerMixin):
'has already been added'.format(self._id, macroscopic)
raise ValueError(msg)
if not isinstance(macroscopic, string_types):
if not isinstance(macroscopic, str):
msg = 'Unable to add a Macroscopic to Material ID="{}" with a ' \
'non-string value "{}"'.format(self._id, macroscopic)
raise ValueError(msg)
@ -465,7 +476,7 @@ class Material(IDManagerMixin):
"""
if not isinstance(macroscopic, string_types):
if not isinstance(macroscopic, str):
msg = 'Unable to remove a Macroscopic "{}" in Material ID="{}" ' \
'since it is not a string'.format(self._id, macroscopic)
raise ValueError(msg)
@ -480,7 +491,7 @@ class Material(IDManagerMixin):
Parameters
----------
element : str
Element to add
Element to add, e.g., 'Zr'
percent : float
Atom or weight percent
percent_type : {'ao', 'wo'}, optional
@ -492,27 +503,15 @@ class Material(IDManagerMixin):
(natural composition).
"""
cv.check_type('nuclide', element, str)
cv.check_type('percent', percent, Real)
cv.check_value('percent type', percent_type, {'ao', 'wo'})
if self._macroscopic is not None:
msg = 'Unable to add an Element to Material ID="{}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(element, string_types):
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'non-string value "{}"'.format(self._id, element)
raise ValueError(msg)
if not isinstance(percent, Real):
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'non-floating point value "{}"'.format(self._id, percent)
raise ValueError(msg)
if percent_type not in ['ao', 'wo']:
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'percent type "{}"'.format(self._id, percent_type)
raise ValueError(msg)
if enrichment is not None:
if not isinstance(enrichment, Real):
msg = 'Unable to add an Element to Material ID="{}" with a ' \
@ -538,10 +537,15 @@ class Material(IDManagerMixin):
format(enrichment, self._id)
warnings.warn(msg)
# Make sure element name is just that
if not element.isalpha():
raise ValueError("Element name should be given by the "
"element's symbol, e.g., 'Zr'")
# Add naturally-occuring isotopes
element = openmc.Element(element)
for nuclide in element.expand(percent, percent_type, enrichment):
self._nuclides.append(nuclide)
self.add_nuclide(*nuclide)
def add_s_alpha_beta(self, name, fraction=1.0):
r"""Add an :math:`S(\alpha,\beta)` table to the material
@ -563,7 +567,7 @@ class Material(IDManagerMixin):
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(name, string_types):
if not isinstance(name, str):
msg = 'Unable to add an S(a,b) table to Material ID="{}" with a ' \
'non-string table name "{}"'.format(self._id, name)
raise ValueError(msg)
@ -688,6 +692,51 @@ class Material(IDManagerMixin):
return nuclides
def get_mass_density(self, nuclide=None):
"""Return mass density of one or all nuclides
Parameters
----------
nuclides : str, optional
Nuclide for which density is desired. If not specified, the density
for the entire material is given.
Returns
-------
float
Density of the nuclide/material in [g/cm^3]
"""
mass_density = 0.0
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
/ openmc.data.AVOGADRO
if nuclide is None or nuclide == nuc:
mass_density += density_i
return mass_density
def get_mass(self, nuclide=None):
"""Return mass of one or all nuclides.
Note that this method requires that the :attr:`Material.volume` has
already been set.
Parameters
----------
nuclides : str, optional
Nuclide for which mass is desired. If not specified, the density
for the entire material is given.
Returns
-------
float
Mass of the nuclide/material in [g]
"""
if self.volume is None:
raise ValueError("Volume must be set in order to determine mass.")
return self.volume*self.get_mass_density(nuclide)
def clone(self, memo=None):
"""Create a copy of this material with a new unique ID.
@ -886,7 +935,7 @@ class Materials(cv.CheckedList):
"""
def __init__(self, materials=None):
super(Materials, self).__init__(Material, 'materials collection')
super().__init__(Material, 'materials collection')
self._cross_sections = None
self._multipole_library = None
@ -903,49 +952,14 @@ class Materials(cv.CheckedList):
@cross_sections.setter
def cross_sections(self, cross_sections):
cv.check_type('cross sections', cross_sections, string_types)
cv.check_type('cross sections', cross_sections, str)
self._cross_sections = cross_sections
@multipole_library.setter
def multipole_library(self, multipole_library):
cv.check_type('cross sections', multipole_library, string_types)
cv.check_type('cross sections', multipole_library, str)
self._multipole_library = multipole_library
def add_material(self, material):
"""Append material to collection
.. deprecated:: 0.8
Use :meth:`Materials.append` instead.
Parameters
----------
material : openmc.Material
Material to add
"""
warnings.warn("Materials.add_material(...) has been deprecated and may be "
"removed in a future version. Use Material.append(...) "
"instead.", DeprecationWarning)
self.append(material)
def add_materials(self, materials):
"""Add multiple materials to the collection
.. deprecated:: 0.8
Use compound assignment instead.
Parameters
----------
materials : Iterable of openmc.Material
Materials to add
"""
warnings.warn("Materials.add_materials(...) has been deprecated and may be "
"removed in a future version. Use compound assignment "
"instead.", DeprecationWarning)
for material in materials:
self.append(material)
def append(self, material):
"""Append material to collection
@ -955,7 +969,7 @@ class Materials(cv.CheckedList):
Material to append
"""
super(Materials, self).append(material)
super().append(material)
def insert(self, index, material):
"""Insert material before index
@ -968,24 +982,7 @@ class Materials(cv.CheckedList):
Material to insert
"""
super(Materials, self).insert(index, material)
def remove_material(self, material):
"""Remove a material from the file
.. deprecated:: 0.8
Use :meth:`Materials.remove` instead.
Parameters
----------
material : openmc.Material
Material to remove
"""
warnings.warn("Materials.remove_material(...) has been deprecated and "
"may be removed in a future version. Use "
"Materials.remove(...) instead.", DeprecationWarning)
self.remove(material)
super().insert(index, material)
def make_isotropic_in_lab(self):
for material in self:

23
openmc/mesh.py
View file

@ -3,7 +3,6 @@ from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
from six import string_types
import numpy as np
import openmc.checkvalue as cv
@ -83,11 +82,29 @@ class Mesh(IDManagerMixin):
def num_mesh_cells(self):
return np.prod(self._dimension)
@property
def indices(self):
ndim = len(self._dimension)
if ndim == 3:
nx, ny, nz = self.dimension
return ((x, y, z)
for z in range(1, nz + 1)
for y in range(1, ny + 1)
for x in range(1, nx + 1))
elif ndim == 2:
nx, ny = self.dimension
return ((x, y)
for y in range(1, ny + 1)
for x in range(1, nx + 1))
else:
nx, = self.dimension
return ((x,) for x in range(1, nx + 1))
@name.setter
def name(self, name):
if name is not None:
cv.check_type('name for mesh ID="{0}"'.format(self._id),
name, string_types)
name, str)
self._name = name
else:
self._name = ''
@ -95,7 +112,7 @@ class Mesh(IDManagerMixin):
@type.setter
def type(self, meshtype):
cv.check_type('type for mesh ID="{0}"'.format(self._id),
meshtype, string_types)
meshtype, str)
cv.check_value('type for mesh ID="{0}"'.format(self._id),
meshtype, ['regular'])
self._type = meshtype

28
openmc/mgxs/library.py
View file

@ -3,10 +3,10 @@ import os
import copy
import pickle
from numbers import Integral
from collections import OrderedDict, Iterable
from collections import OrderedDict
from collections.abc import Iterable
from warnings import warn
from six import string_types
import numpy as np
import openmc
@ -271,7 +271,7 @@ class Library(object):
@name.setter
def name(self, name):
cv.check_type('name', name, string_types)
cv.check_type('name', name, str)
self._name = name
@mgxs_types.setter
@ -280,7 +280,7 @@ class Library(object):
if mgxs_types == 'all':
self._mgxs_types = all_mgxs_types
else:
cv.check_iterable_type('mgxs_types', mgxs_types, string_types)
cv.check_iterable_type('mgxs_types', mgxs_types, str)
for mgxs_type in mgxs_types:
cv.check_value('mgxs_type', mgxs_type, all_mgxs_types)
self._mgxs_types = mgxs_types
@ -587,7 +587,7 @@ class Library(object):
self._nuclides = statepoint.summary.nuclides
if statepoint.run_mode == 'eigenvalue':
self._keff = statepoint.k_combined[0]
self._keff = statepoint.k_combined.n
# Load tallies for each MGXS for each domain and mgxs type
for domain in self.domains:
@ -814,8 +814,8 @@ class Library(object):
'since a statepoint has not yet been loaded'
raise ValueError(msg)
cv.check_type('filename', filename, string_types)
cv.check_type('directory', directory, string_types)
cv.check_type('filename', filename, str)
cv.check_type('directory', directory, str)
import h5py
@ -857,8 +857,8 @@ class Library(object):
"""
cv.check_type('filename', filename, string_types)
cv.check_type('directory', directory, string_types)
cv.check_type('filename', filename, str)
cv.check_type('directory', directory, str)
# Make directory if it does not exist
if not os.path.exists(directory):
@ -892,8 +892,8 @@ class Library(object):
"""
cv.check_type('filename', filename, string_types)
cv.check_type('directory', directory, string_types)
cv.check_type('filename', filename, str)
cv.check_type('directory', directory, str)
# Make directory if it does not exist
if not os.path.exists(directory):
@ -953,8 +953,8 @@ class Library(object):
cv.check_type('domain', domain, (openmc.Material, openmc.Cell,
openmc.Universe, openmc.Mesh))
cv.check_type('xsdata_name', xsdata_name, string_types)
cv.check_type('nuclide', nuclide, string_types)
cv.check_type('xsdata_name', xsdata_name, str)
cv.check_type('nuclide', nuclide, str)
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
if subdomain is not None:
cv.check_iterable_type('subdomain', subdomain, Integral,
@ -1213,7 +1213,7 @@ class Library(object):
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
if xsdata_names is not None:
cv.check_iterable_type('xsdata_names', xsdata_names, string_types)
cv.check_iterable_type('xsdata_names', xsdata_names, str)
# If gathering material-specific data, set the xs_type to macro
if not self.by_nuclide:

101
openmc/mgxs/mdgxs.py
View file

@ -1,6 +1,5 @@
from __future__ import division
from collections import Iterable, OrderedDict
from collections import OrderedDict
from collections.abc import Iterable
import itertools
from numbers import Integral
import warnings
@ -9,7 +8,6 @@ import sys
import copy
from abc import ABCMeta
from six import add_metaclass, string_types
import numpy as np
import openmc
@ -29,7 +27,6 @@ MDGXS_TYPES = ['delayed-nu-fission',
MAX_DELAYED_GROUPS = 8
@add_metaclass(ABCMeta)
class MDGXS(MGXS):
"""An abstract multi-delayed-group cross section for some energy and delayed
group structures within some spatial domain.
@ -133,8 +130,8 @@ class MDGXS(MGXS):
def __init__(self, domain=None, domain_type=None, energy_groups=None,
delayed_groups=None, by_nuclide=False, name='',
num_polar=1, num_azimuthal=1):
super(MDGXS, self).__init__(domain, domain_type, energy_groups,
by_nuclide, name, num_polar, num_azimuthal)
super().__init__(domain, domain_type, energy_groups, by_nuclide, name,
num_polar, num_azimuthal)
self._delayed_groups = None
@ -355,7 +352,7 @@ class MDGXS(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral,
max_depth=3)
for subdomain in subdomains:
@ -363,7 +360,7 @@ class MDGXS(MGXS):
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(groups, string_types):
if not isinstance(groups, str):
cv.check_iterable_type('groups', groups, Integral)
for group in groups:
filters.append(openmc.EnergyFilter)
@ -371,7 +368,7 @@ class MDGXS(MGXS):
(self.energy_groups.get_group_bounds(group),))
# Construct list of delayed group tuples for all requested groups
if not isinstance(delayed_groups, string_types):
if not isinstance(delayed_groups, str):
cv.check_type('delayed groups', delayed_groups, list, int)
for delayed_group in delayed_groups:
filters.append(openmc.DelayedGroupFilter)
@ -475,7 +472,7 @@ class MDGXS(MGXS):
"""
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
cv.check_iterable_type('energy_groups', groups, Integral)
cv.check_type('delayed groups', delayed_groups, list, int)
@ -551,7 +548,7 @@ class MDGXS(MGXS):
"""
merged_mdgxs = super(MDGXS, self).merge(other)
merged_mdgxs = super().merge(other)
# Merge delayed groups
if self.delayed_groups != other.delayed_groups:
@ -581,11 +578,11 @@ class MDGXS(MGXS):
"""
if self.delayed_groups is None:
super(MDGXS, self).print_xs(subdomains, nuclides, xs_type)
super().print_xs(subdomains, nuclides, xs_type)
return
# Construct a collection of the subdomains to report
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -602,7 +599,7 @@ class MDGXS(MGXS):
elif nuclides == 'sum':
nuclides = ['sum']
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
else:
nuclides = ['sum']
@ -725,8 +722,8 @@ class MDGXS(MGXS):
"""
cv.check_type('filename', filename, string_types)
cv.check_type('directory', directory, string_types)
cv.check_type('filename', filename, str)
cv.check_type('directory', directory, str)
cv.check_value('format', format, ['csv', 'excel', 'pickle', 'latex'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@ -816,11 +813,11 @@ class MDGXS(MGXS):
"""
if not isinstance(groups, string_types):
if not isinstance(groups, str):
cv.check_iterable_type('groups', groups, Integral)
if nuclides != 'all' and nuclides != 'sum':
cv.check_iterable_type('nuclides', nuclides, string_types)
if not isinstance(delayed_groups, string_types):
cv.check_iterable_type('nuclides', nuclides, str)
if not isinstance(delayed_groups, str):
cv.check_type('delayed groups', delayed_groups, list, int)
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@ -858,7 +855,7 @@ class MDGXS(MGXS):
columns = self._df_convert_columns_to_bins(df)
# Select out those groups the user requested
if not isinstance(groups, string_types):
if not isinstance(groups, str):
if 'group in' in df:
df = df[df['group in'].isin(groups)]
if 'group out' in df:
@ -1011,9 +1008,8 @@ class ChiDelayed(MDGXS):
def __init__(self, domain=None, domain_type=None, energy_groups=None,
delayed_groups=None, by_nuclide=False, name='',
num_polar=1, num_azimuthal=1):
super(ChiDelayed, self).__init__(domain, domain_type, energy_groups,
delayed_groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, energy_groups, delayed_groups,
by_nuclide, name, num_polar, num_azimuthal)
self._rxn_type = 'chi-delayed'
self._estimator = 'analog'
@ -1059,7 +1055,7 @@ class ChiDelayed(MDGXS):
# Compute chi
self._xs_tally = self.rxn_rate_tally / delayed_nu_fission_in
super(ChiDelayed, self)._compute_xs()
super()._compute_xs()
# Add the coarse energy filter back to the nu-fission tally
delayed_nu_fission_in.filters.append(energy_filter)
@ -1131,8 +1127,7 @@ class ChiDelayed(MDGXS):
delayed_nu_fission_in.remove_filter(energy_filter)
# Call super class method and null out derived tallies
slice_xs = super(ChiDelayed, self).get_slice(nuclides, groups,
delayed_groups)
slice_xs = super().get_slice(nuclides, groups, delayed_groups)
slice_xs._rxn_rate_tally = None
slice_xs._xs_tally = None
@ -1288,7 +1283,7 @@ class ChiDelayed(MDGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral,
max_depth=3)
for subdomain in subdomains:
@ -1296,7 +1291,7 @@ class ChiDelayed(MDGXS):
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(groups, string_types):
if not isinstance(groups, str):
cv.check_iterable_type('groups', groups, Integral)
for group in groups:
filters.append(openmc.EnergyoutFilter)
@ -1304,7 +1299,7 @@ class ChiDelayed(MDGXS):
(self.energy_groups.get_group_bounds(group),))
# Construct list of delayed group tuples for all requested groups
if not isinstance(delayed_groups, string_types):
if not isinstance(delayed_groups, str):
cv.check_type('delayed groups', delayed_groups, list, int)
for delayed_group in delayed_groups:
filters.append(openmc.DelayedGroupFilter)
@ -1352,7 +1347,7 @@ class ChiDelayed(MDGXS):
# Get chi delayed for user-specified nuclides in the domain
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
xs = self.xs_tally.get_values(filters=filters,
filter_bins=filter_bins,
nuclides=nuclides, value=value)
@ -1525,10 +1520,8 @@ class DelayedNuFissionXS(MDGXS):
def __init__(self, domain=None, domain_type=None, energy_groups=None,
delayed_groups=None, by_nuclide=False, name='',
num_polar=1, num_azimuthal=1):
super(DelayedNuFissionXS, self).__init__(domain, domain_type,
energy_groups, delayed_groups,
by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, energy_groups, delayed_groups,
by_nuclide, name, num_polar, num_azimuthal)
self._rxn_type = 'delayed-nu-fission'
@ -1661,9 +1654,8 @@ class Beta(MDGXS):
def __init__(self, domain=None, domain_type=None, energy_groups=None,
delayed_groups=None, by_nuclide=False, name='',
num_polar=1, num_azimuthal=1):
super(Beta, self).__init__(domain, domain_type, energy_groups,
delayed_groups, by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, energy_groups, delayed_groups,
by_nuclide, name, num_polar, num_azimuthal)
self._rxn_type = 'beta'
@property
@ -1689,7 +1681,7 @@ class Beta(MDGXS):
# Compute beta
self._xs_tally = self.rxn_rate_tally / nu_fission
super(Beta, self)._compute_xs()
super()._compute_xs()
return self._xs_tally
@ -1845,9 +1837,8 @@ class DecayRate(MDGXS):
def __init__(self, domain=None, domain_type=None, energy_groups=None,
delayed_groups=None, by_nuclide=False, name='',
num_polar=1, num_azimuthal=1):
super(DecayRate, self).__init__(domain, domain_type, energy_groups,
delayed_groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, energy_groups, delayed_groups,
by_nuclide, name, num_polar, num_azimuthal)
self._rxn_type = 'decay-rate'
@property
@ -1882,7 +1873,7 @@ class DecayRate(MDGXS):
# Compute the decay rate
self._xs_tally = self.rxn_rate_tally / delayed_nu_fission
super(DecayRate, self)._compute_xs()
super()._compute_xs()
return self._xs_tally
@ -1914,7 +1905,6 @@ class DecayRate(MDGXS):
return self._get_homogenized_mgxs(other_mgxs, 'delayed-nu-fission')
@add_metaclass(ABCMeta)
class MatrixMDGXS(MDGXS):
"""An abstract multi-delayed-group cross section for some energy group and
delayed group structure within some spatial domain. This class is
@ -2117,7 +2107,7 @@ class MatrixMDGXS(MDGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral,
max_depth=3)
for subdomain in subdomains:
@ -2125,7 +2115,7 @@ class MatrixMDGXS(MDGXS):
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(in_groups, string_types):
if not isinstance(in_groups, str):
cv.check_iterable_type('groups', in_groups, Integral)
for group in in_groups:
filters.append(openmc.EnergyFilter)
@ -2133,7 +2123,7 @@ class MatrixMDGXS(MDGXS):
self.energy_groups.get_group_bounds(group),))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(out_groups, string_types):
if not isinstance(out_groups, str):
cv.check_iterable_type('groups', out_groups, Integral)
for group in out_groups:
filters.append(openmc.EnergyoutFilter)
@ -2141,7 +2131,7 @@ class MatrixMDGXS(MDGXS):
self.energy_groups.get_group_bounds(group),))
# Construct list of delayed group tuples for all requested groups
if not isinstance(delayed_groups, string_types):
if not isinstance(delayed_groups, str):
cv.check_type('delayed groups', delayed_groups, list, int)
for delayed_group in delayed_groups:
filters.append(openmc.DelayedGroupFilter)
@ -2266,8 +2256,7 @@ class MatrixMDGXS(MDGXS):
"""
# Call super class method and null out derived tallies
slice_xs = super(MatrixMDGXS, self).get_slice(nuclides, in_groups,
delayed_groups)
slice_xs = super().get_slice(nuclides, in_groups, delayed_groups)
slice_xs._rxn_rate_tally = None
slice_xs._xs_tally = None
@ -2312,7 +2301,7 @@ class MatrixMDGXS(MDGXS):
"""
# Construct a collection of the subdomains to report
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -2329,7 +2318,7 @@ class MatrixMDGXS(MDGXS):
if nuclides == 'sum':
nuclides = ['sum']
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
else:
nuclides = ['sum']
@ -2614,12 +2603,8 @@ class DelayedNuFissionMatrixXS(MatrixMDGXS):
def __init__(self, domain=None, domain_type=None, energy_groups=None,
delayed_groups=None, by_nuclide=False, name='',
num_polar=1, num_azimuthal=1):
super(DelayedNuFissionMatrixXS, self).__init__(domain, domain_type,
energy_groups,
delayed_groups,
by_nuclide, name,
num_polar,
num_azimuthal)
super().__init__(domain, domain_type, energy_groups, delayed_groups,
by_nuclide, name, num_polar, num_azimuthal)
self._rxn_type = 'delayed-nu-fission'
self._hdf5_key = 'delayed-nu-fission matrix'
self._estimator = 'analog'

183
openmc/mgxs/mgxs.py
View file

@ -1,15 +1,12 @@
from __future__ import division
from collections import OrderedDict
from numbers import Integral
import warnings
import os
import copy
from abc import ABCMeta
import itertools
from six import add_metaclass, string_types
import numpy as np
import h5py
import openmc
import openmc.checkvalue as cv
@ -115,8 +112,7 @@ def _df_column_convert_to_bin(df, current_name, new_name, values_to_bin,
df.rename(columns={current_name: new_name}, inplace=True)
@add_metaclass(ABCMeta)
class MGXS(object):
class MGXS(metaclass=ABCMeta):
"""An abstract multi-group cross section for some energy group structure
within some spatial domain.
@ -579,7 +575,7 @@ class MGXS(object):
@name.setter
def name(self, name):
cv.check_type('name', name, string_types)
cv.check_type('name', name, str)
self._name = name
@by_nuclide.setter
@ -589,7 +585,7 @@ class MGXS(object):
@nuclides.setter
def nuclides(self, nuclides):
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
self._nuclides = nuclides
@estimator.setter
@ -805,7 +801,7 @@ class MGXS(object):
"""
cv.check_type('nuclide', nuclide, string_types)
cv.check_type('nuclide', nuclide, str)
# Get list of all nuclides in the spatial domain
nuclides = self.domain.get_nuclide_densities()
@ -940,8 +936,7 @@ class MGXS(object):
# NOTE: This is important if tally merging was used
if self.domain_type == 'mesh':
filters = [_DOMAIN_TO_FILTER[self.domain_type]]
xyz = [range(1, x + 1) for x in self.domain.dimension]
filter_bins = [tuple(itertools.product(*xyz))]
filter_bins = [tuple(self.domain.indices)]
elif self.domain_type != 'distribcell':
filters = [_DOMAIN_TO_FILTER[self.domain_type]]
filter_bins = [(self.domain.id,)]
@ -1032,7 +1027,7 @@ class MGXS(object):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral,
max_depth=3)
@ -1043,7 +1038,7 @@ class MGXS(object):
filter_bins.append(tuple(subdomain_bins))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(groups, string_types):
if not isinstance(groups, str):
cv.check_iterable_type('groups', groups, Integral)
filters.append(openmc.EnergyFilter)
energy_bins = []
@ -1218,7 +1213,7 @@ class MGXS(object):
"""
# Construct a collection of the subdomain filter bins to average across
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral)
subdomains = [(subdomain,) for subdomain in subdomains]
subdomains = [tuple(subdomains)]
@ -1375,7 +1370,7 @@ class MGXS(object):
"""
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
cv.check_iterable_type('energy_groups', groups, Integral)
# Build lists of filters and filter bins to slice
@ -1529,13 +1524,12 @@ class MGXS(object):
"""
# Construct a collection of the subdomains to report
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
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 == 'mesh':
xyz = [range(1, x + 1) for x in self.domain.dimension]
subdomains = list(itertools.product(*xyz))
subdomains = list(self.domain.indices)
else:
subdomains = [self.domain.id]
@ -1546,7 +1540,7 @@ class MGXS(object):
elif nuclides == 'sum':
nuclides = ['sum']
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
else:
nuclides = ['sum']
@ -1682,13 +1676,8 @@ class MGXS(object):
ValueError
When this method is called before the multi-group cross section is
computed from tally data.
ImportError
When h5py is not installed.
"""
import h5py
# Make directory if it does not exist
if not os.path.exists(directory):
os.makedirs(directory)
@ -1702,7 +1691,7 @@ class MGXS(object):
xs_results = h5py.File(filename, 'w', libver=libver)
# Construct a collection of the subdomains to report
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1710,8 +1699,7 @@ class MGXS(object):
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]
subdomains = list(itertools.product(*xyz))
subdomains = list(self.domain.indices)
else:
subdomains = [self.domain.id]
@ -1723,7 +1711,7 @@ class MGXS(object):
elif nuclides == 'sum':
nuclides = ['sum']
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
else:
nuclides = ['sum']
@ -1801,8 +1789,8 @@ class MGXS(object):
"""
cv.check_type('filename', filename, string_types)
cv.check_type('directory', directory, string_types)
cv.check_type('filename', filename, str)
cv.check_type('directory', directory, str)
cv.check_value('format', format, ['csv', 'excel', 'pickle', 'latex'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@ -1888,10 +1876,10 @@ class MGXS(object):
"""
if not isinstance(groups, string_types):
if not isinstance(groups, str):
cv.check_iterable_type('groups', groups, Integral)
if nuclides != 'all' and nuclides != 'sum':
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
# Get a Pandas DataFrame from the derived xs tally
@ -1927,7 +1915,7 @@ class MGXS(object):
columns = self._df_convert_columns_to_bins(df)
# Select out those groups the user requested
if not isinstance(groups, string_types):
if not isinstance(groups, str):
if 'group in' in df:
df = df[df['group in'].isin(groups)]
if 'group out' in df:
@ -1980,7 +1968,6 @@ class MGXS(object):
return 'cm^-1' if xs_type == 'macro' else 'barns'
@add_metaclass(ABCMeta)
class MatrixMGXS(MGXS):
"""An abstract multi-group cross section for some energy group structure
within some spatial domain. This class is specifically intended for
@ -2168,7 +2155,7 @@ class MatrixMGXS(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral,
max_depth=3)
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
@ -2178,7 +2165,7 @@ class MatrixMGXS(MGXS):
filter_bins.append(tuple(subdomain_bins))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(in_groups, string_types):
if not isinstance(in_groups, str):
cv.check_iterable_type('groups', in_groups, Integral)
filters.append(openmc.EnergyFilter)
for group in in_groups:
@ -2186,7 +2173,7 @@ class MatrixMGXS(MGXS):
filter_bins.append(tuple(energy_bins))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(out_groups, string_types):
if not isinstance(out_groups, str):
cv.check_iterable_type('groups', out_groups, Integral)
for group in out_groups:
filters.append(openmc.EnergyoutFilter)
@ -2301,7 +2288,7 @@ class MatrixMGXS(MGXS):
"""
# Call super class method and null out derived tallies
slice_xs = super(MatrixMGXS, self).get_slice(nuclides, in_groups)
slice_xs = super().get_slice(nuclides, in_groups)
slice_xs._rxn_rate_tally = None
slice_xs._xs_tally = None
@ -2346,13 +2333,12 @@ class MatrixMGXS(MGXS):
"""
# Construct a collection of the subdomains to report
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
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 == 'mesh':
xyz = [range(1, x + 1) for x in self.domain.dimension]
subdomains = list(itertools.product(*xyz))
subdomains = list(self.domain.indices)
else:
subdomains = [self.domain.id]
@ -2363,7 +2349,7 @@ class MatrixMGXS(MGXS):
if nuclides == 'sum':
nuclides = ['sum']
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
else:
nuclides = ['sum']
@ -2576,9 +2562,8 @@ class TotalXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(TotalXS, self).__init__(domain, domain_type,
groups, by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._rxn_type = 'total'
@ -2713,9 +2698,8 @@ class TransportXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None, nu=False,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(TransportXS, self).__init__(domain, domain_type,
groups, by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
# Use tracklength estimators for the total MGXS term, and
# analog estimators for the transport correction term
@ -2724,7 +2708,7 @@ class TransportXS(MGXS):
self.nu = nu
def __deepcopy__(self, memo):
clone = super(TransportXS, self).__deepcopy__(memo)
clone = super().__deepcopy__(memo)
clone._nu = self.nu
return clone
@ -2921,9 +2905,8 @@ class AbsorptionXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(AbsorptionXS, self).__init__(domain, domain_type,
groups, by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._rxn_type = 'absorption'
@ -3048,9 +3031,8 @@ class CaptureXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(CaptureXS, self).__init__(domain, domain_type,
groups, by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._rxn_type = 'capture'
@property
@ -3203,16 +3185,15 @@ class FissionXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None, nu=False,
prompt=False, by_nuclide=False, name='', num_polar=1,
num_azimuthal=1):
super(FissionXS, self).__init__(domain, domain_type,
groups, by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._nu = False
self._prompt = False
self.nu = nu
self.prompt = prompt
def __deepcopy__(self, memo):
clone = super(FissionXS, self).__deepcopy__(memo)
clone = super().__deepcopy__(memo)
clone._nu = self.nu
clone._prompt = self.prompt
return clone
@ -3371,9 +3352,8 @@ class KappaFissionXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(KappaFissionXS, self).__init__(domain, domain_type,
groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._rxn_type = 'kappa-fission'
@ -3504,13 +3484,12 @@ class ScatterXS(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1,
num_azimuthal=1, nu=False):
super(ScatterXS, self).__init__(domain, domain_type,
groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self.nu = nu
def __deepcopy__(self, memo):
clone = super(ScatterXS, self).__deepcopy__(memo)
clone = super().__deepcopy__(memo)
clone._nu = self.nu
return clone
@ -3721,9 +3700,8 @@ class ScatterMatrixXS(MatrixMGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1,
num_azimuthal=1, nu=False):
super(ScatterMatrixXS, self).__init__(domain, domain_type,
groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._formulation = 'simple'
self._correction = 'P0'
self._scatter_format = 'legendre'
@ -3734,7 +3712,7 @@ class ScatterMatrixXS(MatrixMGXS):
self.nu = nu
def __deepcopy__(self, memo):
clone = super(ScatterMatrixXS, self).__deepcopy__(memo)
clone = super().__deepcopy__(memo)
clone._formulation = self.formulation
clone._correction = self.correction
clone._scatter_format = self.scatter_format
@ -3825,7 +3803,7 @@ class ScatterMatrixXS(MatrixMGXS):
@property
def tally_keys(self):
if self.formulation == 'simple':
return super(ScatterMatrixXS, self).tally_keys
return super().tally_keys
else:
# Add keys for groupwise scattering cross section
tally_keys = ['flux (tracklength)', 'scatter']
@ -4155,7 +4133,7 @@ class ScatterMatrixXS(MatrixMGXS):
[score_prefix + '{}'.format(i)
for i in range(self.legendre_order + 1)]
super(ScatterMatrixXS, self).load_from_statepoint(statepoint)
super().load_from_statepoint(statepoint)
def get_slice(self, nuclides=[], in_groups=[], out_groups=[],
legendre_order='same'):
@ -4195,7 +4173,7 @@ class ScatterMatrixXS(MatrixMGXS):
"""
# Call super class method and null out derived tallies
slice_xs = super(ScatterMatrixXS, self).get_slice(nuclides, in_groups)
slice_xs = super().get_slice(nuclides, in_groups)
slice_xs._rxn_rate_tally = None
slice_xs._xs_tally = None
@ -4311,7 +4289,7 @@ class ScatterMatrixXS(MatrixMGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3)
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
subdomain_bins = []
@ -4320,7 +4298,7 @@ class ScatterMatrixXS(MatrixMGXS):
filter_bins.append(tuple(subdomain_bins))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(in_groups, string_types):
if not isinstance(in_groups, str):
cv.check_iterable_type('groups', in_groups, Integral)
filters.append(openmc.EnergyFilter)
energy_bins = []
@ -4330,7 +4308,7 @@ class ScatterMatrixXS(MatrixMGXS):
filter_bins.append(tuple(energy_bins))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(out_groups, string_types):
if not isinstance(out_groups, str):
cv.check_iterable_type('groups', out_groups, Integral)
for group in out_groups:
filters.append(openmc.EnergyoutFilter)
@ -4486,8 +4464,7 @@ class ScatterMatrixXS(MatrixMGXS):
"""
df = super(ScatterMatrixXS, self).get_pandas_dataframe(
groups, nuclides, xs_type, paths)
df = super().get_pandas_dataframe(groups, nuclides, xs_type, paths)
if self.scatter_format == 'legendre':
# Add a moment column to dataframe
@ -4543,13 +4520,12 @@ class ScatterMatrixXS(MatrixMGXS):
"""
# Construct a collection of the subdomains to report
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
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 == 'mesh':
xyz = [range(1, x + 1) for x in self.domain.dimension]
subdomains = list(itertools.product(*xyz))
subdomains = list(self.domain.indices)
else:
subdomains = [self.domain.id]
@ -4560,7 +4536,7 @@ class ScatterMatrixXS(MatrixMGXS):
if nuclides == 'sum':
nuclides = ['sum']
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
else:
nuclides = ['sum']
@ -4811,9 +4787,8 @@ class MultiplicityMatrixXS(MatrixMGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(MultiplicityMatrixXS, self).__init__(domain, domain_type, groups,
by_nuclide, name, num_polar,
num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._rxn_type = 'multiplicity matrix'
self._estimator = 'analog'
self._valid_estimators = ['analog']
@ -4848,7 +4823,7 @@ class MultiplicityMatrixXS(MatrixMGXS):
# Compute the multiplicity
self._xs_tally = self.rxn_rate_tally / scatter
super(MultiplicityMatrixXS, self)._compute_xs()
super()._compute_xs()
return self._xs_tally
@ -4977,9 +4952,8 @@ class ScatterProbabilityMatrix(MatrixMGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(ScatterProbabilityMatrix, self).__init__(
domain, domain_type, groups, by_nuclide,
name, num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide,
name, num_polar, num_azimuthal)
self._rxn_type = 'scatter'
self._hdf5_key = 'scatter probability matrix'
@ -5021,7 +4995,7 @@ class ScatterProbabilityMatrix(MatrixMGXS):
# Compute the group-to-group probabilities
self._xs_tally = self.tallies[self.rxn_type] / norm
super(ScatterProbabilityMatrix, self)._compute_xs()
super()._compute_xs()
return self._xs_tally
@ -5151,9 +5125,8 @@ class NuFissionMatrixXS(MatrixMGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1,
num_azimuthal=1, prompt=False):
super(NuFissionMatrixXS, self).__init__(domain, domain_type,
groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
if not prompt:
self._rxn_type = 'nu-fission'
self._hdf5_key = 'nu-fission matrix'
@ -5174,7 +5147,7 @@ class NuFissionMatrixXS(MatrixMGXS):
self._prompt = prompt
def __deepcopy__(self, memo):
clone = super(NuFissionMatrixXS, self).__deepcopy__(memo)
clone = super().__deepcopy__(memo)
clone._prompt = self.prompt
return clone
@ -5308,8 +5281,8 @@ class Chi(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
prompt=False, by_nuclide=False, name='', num_polar=1,
num_azimuthal=1):
super(Chi, self).__init__(domain, domain_type, groups, by_nuclide,
name, num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
if not prompt:
self._rxn_type = 'chi'
else:
@ -5319,7 +5292,7 @@ class Chi(MGXS):
self.prompt = prompt
def __deepcopy__(self, memo):
clone = super(Chi, self).__deepcopy__(memo)
clone = super().__deepcopy__(memo)
clone._prompt = self.prompt
return clone
@ -5449,7 +5422,7 @@ class Chi(MGXS):
nu_fission_in.remove_filter(energy_filter)
# Call super class method and null out derived tallies
slice_xs = super(Chi, self).get_slice(nuclides, groups)
slice_xs = super().get_slice(nuclides, groups)
slice_xs._rxn_rate_tally = None
slice_xs._xs_tally = None
@ -5586,7 +5559,7 @@ class Chi(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if not isinstance(subdomains, string_types):
if not isinstance(subdomains, str):
cv.check_iterable_type('subdomains', subdomains, Integral,
max_depth=3)
filters.append(_DOMAIN_TO_FILTER[self.domain_type])
@ -5596,7 +5569,7 @@ class Chi(MGXS):
filter_bins.append(tuple(subdomain_bins))
# Construct list of energy group bounds tuples for all requested groups
if not isinstance(groups, string_types):
if not isinstance(groups, str):
cv.check_iterable_type('groups', groups, Integral)
filters.append(openmc.EnergyoutFilter)
energy_bins = []
@ -5644,7 +5617,7 @@ class Chi(MGXS):
# Get chi for user-specified nuclides in the domain
else:
cv.check_iterable_type('nuclides', nuclides, string_types)
cv.check_iterable_type('nuclides', nuclides, str)
xs = self.xs_tally.get_values(filters=filters,
filter_bins=filter_bins,
nuclides=nuclides, value=value)
@ -5727,8 +5700,7 @@ class Chi(MGXS):
"""
# Build the dataframe using the parent class method
df = super(Chi, self).get_pandas_dataframe(
groups, nuclides, xs_type, paths=paths)
df = super().get_pandas_dataframe(groups, nuclides, xs_type, paths=paths)
# If user requested micro cross sections, multiply by the atom
# densities to cancel out division made by the parent class method
@ -5886,9 +5858,8 @@ class InverseVelocity(MGXS):
def __init__(self, domain=None, domain_type=None, groups=None,
by_nuclide=False, name='', num_polar=1, num_azimuthal=1):
super(InverseVelocity, self).__init__(domain, domain_type,
groups, by_nuclide, name,
num_polar, num_azimuthal)
super().__init__(domain, domain_type, groups, by_nuclide, name,
num_polar, num_azimuthal)
self._rxn_type = 'inverse-velocity'
def get_units(self, xs_type='macro'):

5
openmc/mgxs_library.py
View file

@ -2,7 +2,6 @@ import copy
from numbers import Real, Integral
import os
from six import string_types
import numpy as np
import h5py
from scipy.interpolate import interp1d
@ -381,7 +380,7 @@ class XSdata(object):
@name.setter
def name(self, name):
check_type('name for XSdata', name, string_types)
check_type('name for XSdata', name, str)
self._name = name
@energy_groups.setter
@ -2517,7 +2516,7 @@ class MGXSLibrary(object):
"""
check_type('filename', filename, string_types)
check_type('filename', filename, str)
# Create and write to the HDF5 file
file = h5py.File(filename, "w", libver=libver)

14
openmc/mixin.py
View file

@ -45,14 +45,24 @@ class IDManagerMixin(object):
@id.setter
def id(self, uid):
cls = type(self)
name = cls.__name__
# The first time this is called for a class, we search through the MRO
# to determine which class actually holds next_id and used_ids. Since
# next_id is an integer (immutable), we can't modify it directly through
# the instance without just creating a new attribute
try:
cls = self._id_class
except AttributeError:
for cls in self.__class__.__mro__:
if 'next_id' in cls.__dict__:
break
if uid is None:
while cls.next_id in cls.used_ids:
cls.next_id += 1
self._id = cls.next_id
cls.used_ids.add(cls.next_id)
else:
name = cls.__name__
cv.check_type('{} ID'.format(name), uid, Integral)
cv.check_greater_than('{} ID'.format(name), uid, 0, equality=True)
if uid in cls.used_ids:

97
openmc/model/funcs.py
View file

@ -1,10 +1,103 @@
from __future__ import division
from collections import Iterable, OrderedDict
from collections import OrderedDict
from collections.abc import Iterable
from math import sqrt
from numbers import Real
from openmc import XPlane, YPlane, Plane, ZCylinder
from openmc.checkvalue import check_type, check_value
import openmc.data
def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K',
press_unit='MPa', density=None, **kwargs):
"""Return a Material with the composition of boron dissolved in water.
The water density can be determined from a temperature and pressure, or it
can be set directly.
The concentration of boron has no effect on the stoichiometric ratio of H
and O---they are fixed at 2-1.
Parameters
----------
boron_ppm : float
The weight fraction in parts-per-million of elemental boron in the
water.
temperature : float
Temperature in [K] used to compute water density.
pressure : float
Pressure in [MPa] used to compute water density.
temp_unit : {'K', 'C', 'F'}
The units used for the `temperature` argument.
press_unit : {'MPa', 'psi'}
The units used for the `pressure` argument.
density : float
Water density in [g / cm^3]. If specified, this value overrides the
temperature and pressure arguments.
**kwargs
All keyword arguments are passed to the created Material object.
Returns
-------
openmc.Material
"""
# Perform any necessary unit conversions.
check_value('temperature unit', temp_unit, ('K', 'C', 'F'))
if temp_unit == 'K':
T = temperature
elif temp_unit == 'C':
T = temperature + 273.15
elif temp_unit == 'F':
T = (temperature + 459.67) * 5.0 / 9.0
check_value('pressure unit', press_unit, ('MPa', 'psi'))
if press_unit == 'MPa':
P = pressure
elif press_unit == 'psi':
P = pressure * 0.006895
# Set the density of water, either from an explicitly given density or from
# temperature and pressure.
if density is not None:
water_density = density
else:
water_density = openmc.data.water_density(T, P)
# Compute the density of the solution.
solution_density = water_density / (1 - boron_ppm * 1e-6)
# Compute the molar mass of pure water.
hydrogen = openmc.Element('H')
oxygen = openmc.Element('O')
M_H2O = 0.0
for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'):
M_H2O += frac * openmc.data.atomic_mass(iso_name)
for iso_name, frac, junk in oxygen.expand(1.0, 'ao'):
M_H2O += frac * openmc.data.atomic_mass(iso_name)
# Compute the molar mass of boron.
boron = openmc.Element('B')
M_B = 0.0
for iso_name, frac, junk in boron.expand(1.0, 'ao'):
M_B += frac * openmc.data.atomic_mass(iso_name)
# Compute the number fractions of each element.
frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O
frac_H = 2 * frac_H2O
frac_O = frac_H2O
frac_B = boron_ppm * 1e-6 / M_B
# Build the material.
if density is None:
out = openmc.Material(temperature=T, **kwargs)
else:
out = openmc.Material(**kwargs)
out.add_element('H', frac_H, 'ao')
out.add_element('O', frac_O, 'ao')
out.add_element('B', frac_B, 'ao')
out.set_density('g/cc', solution_density)
out.add_s_alpha_beta('c_H_in_H2O')
return out
def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),

56
openmc/model/model.py
View file

@ -1,7 +1,7 @@
from collections import Iterable
from collections.abc import Iterable
import openmc
from openmc.checkvalue import check_type
from openmc.checkvalue import check_type, check_value
class Model(object):
@ -136,9 +136,49 @@ class Model(object):
for plot in plots:
self._plots.append(plot)
def export_to_xml(self):
"""Export model to XML files.
def deplete(self, timesteps, power, chain_file=None, method='cecm',
**kwargs):
"""Deplete model using specified timesteps/power
Parameters
----------
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
power : float or iterable of float
Power of the reactor in [W]. A single value indicates that the power
is constant over all timesteps. An iterable indicates potentially
different power levels for each timestep. For a 2D problem, the
power can be given in [W/cm] as long as the "volume" assigned to a
depletion material is actually an area in [cm^2].
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
method : {'cecm', 'predictor'}
Integration method used for depletion
**kwargs
Keyword arguments passed to integration function (e.g.,
:func:`openmc.deplete.integrator.cecm`)
"""
# Import the depletion module. This is done here rather than the module
# header to delay importing openmc.capi (through openmc.deplete) which
# can be tough to install properly.
import openmc.deplete as dep
# Create OpenMC transport operator
op = dep.Operator(self.geometry, self.settings, chain_file)
# Perform depletion
if method == 'predictor':
dep.integrator.predictor(op, timesteps, power, **kwargs)
elif method == 'cecm':
dep.integrator.cecm(op, timesteps, power, **kwargs)
else:
check_value('method', method, ('cecm', 'predictor'))
def export_to_xml(self):
"""Export model to XML files."""
self.settings.export_to_xml()
self.geometry.export_to_xml()
@ -166,19 +206,17 @@ class Model(object):
Parameters
----------
**kwargs
All keyword arguments are passed to openmc.run
All keyword arguments are passed to :func:`openmc.run`
Returns
-------
2-tuple of float
uncertainties.UFloat
Combined estimator of k-effective from the statepoint
"""
self.export_to_xml()
return_code = openmc.run(**kwargs)
assert (return_code == 0), "OpenMC did not execute successfully"
openmc.run(**kwargs)
n = self.settings.batches
if self.settings.statepoint is not None:

26
openmc/model/triso.py
View file

@ -1,22 +1,17 @@
from __future__ import division
import copy
import warnings
import itertools
import random
from collections import Iterable, defaultdict
from collections import defaultdict
from collections.abc import Iterable
from numbers import Real
from random import uniform, gauss
from heapq import heappush, heappop
from math import pi, sin, cos, floor, log10, sqrt
from abc import ABCMeta, abstractproperty, abstractmethod
from six import add_metaclass
import numpy as np
try:
import scipy.spatial
_SCIPY_AVAILABLE = True
except ImportError:
_SCIPY_AVAILABLE = False
import scipy.spatial
import openmc
import openmc.checkvalue as cv
@ -51,7 +46,7 @@ class TRISO(openmc.Cell):
def __init__(self, outer_radius, fill, center=(0., 0., 0.)):
self._surface = openmc.Sphere(R=outer_radius)
super(TRISO, self).__init__(fill=fill, region=-self._surface)
super().__init__(fill=fill, region=-self._surface)
self.center = np.asarray(center)
@property
@ -96,8 +91,7 @@ class TRISO(openmc.Cell):
k_min:k_max+1, j_min:j_max+1, i_min:i_max+1]))
@add_metaclass(ABCMeta)
class _Domain(object):
class _Domain(metaclass=ABCMeta):
"""Container in which to pack particles.
Parameters
@ -252,7 +246,7 @@ class _CubicDomain(_Domain):
"""
def __init__(self, length, particle_radius, center=[0., 0., 0.]):
super(_CubicDomain, self).__init__(particle_radius, center)
super().__init__(particle_radius, center)
self.length = length
@property
@ -331,7 +325,7 @@ class _CylindricalDomain(_Domain):
"""
def __init__(self, length, radius, particle_radius, center=[0., 0., 0.]):
super(_CylindricalDomain, self).__init__(particle_radius, center)
super().__init__(particle_radius, center)
self.length = length
self.radius = radius
@ -421,7 +415,7 @@ class _SphericalDomain(_Domain):
"""
def __init__(self, radius, particle_radius, center=[0., 0., 0.]):
super(_SphericalDomain, self).__init__(particle_radius, center)
super().__init__(particle_radius, center)
self.radius = radius
@property
@ -837,10 +831,6 @@ def _close_random_pack(domain, particles, contraction_rate):
if rods:
inner_diameter[0] = rods[0][0]
if not _SCIPY_AVAILABLE:
raise ImportError('SciPy must be installed to perform '
'close random packing.')
n_particles = len(particles)
diameter = 2*domain.particle_radius

4
openmc/nuclide.py
View file

@ -1,7 +1,5 @@
import warnings
from six import string_types
import openmc.checkvalue as cv
@ -34,7 +32,7 @@ class Nuclide(str):
'"{}" is being renamed as "{}".'.format(orig_name, name)
warnings.warn(msg)
return super(Nuclide, cls).__new__(cls, name)
return super().__new__(cls, name)
@property
def name(self):

107
openmc/plots.py
View file

@ -1,10 +1,10 @@
from collections import Iterable, Mapping
from collections.abc import Iterable, Mapping
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import subprocess
import sys
import warnings
from six import string_types
import numpy as np
import openmc
@ -297,7 +297,7 @@ class Plot(IDManagerMixin):
@name.setter
def name(self, name):
cv.check_type('plot name', name, string_types)
cv.check_type('plot name', name, str)
self._name = name
@width.setter
@ -322,7 +322,7 @@ class Plot(IDManagerMixin):
@filename.setter
def filename(self, filename):
cv.check_type('filename', filename, string_types)
cv.check_type('filename', filename, str)
self._filename = filename
@color_by.setter
@ -343,7 +343,7 @@ class Plot(IDManagerMixin):
@background.setter
def background(self, background):
cv.check_type('plot background', background, Iterable)
if isinstance(background, string_types):
if isinstance(background, str):
if background.lower() not in _SVG_COLORS:
raise ValueError("'{}' is not a valid color.".format(background))
else:
@ -359,7 +359,7 @@ class Plot(IDManagerMixin):
for key, value in colors.items():
cv.check_type('plot color key', key, (openmc.Cell, openmc.Material))
cv.check_type('plot color value', value, Iterable)
if isinstance(value, string_types):
if isinstance(value, str):
if value.lower() not in _SVG_COLORS:
raise ValueError("'{}' is not a valid color.".format(value))
else:
@ -380,7 +380,7 @@ class Plot(IDManagerMixin):
@mask_background.setter
def mask_background(self, mask_background):
cv.check_type('plot mask background', mask_background, Iterable)
if isinstance(mask_background, string_types):
if isinstance(mask_background, str):
if mask_background.lower() not in _SVG_COLORS:
raise ValueError("'{}' is not a valid color.".format(mask_background))
else:
@ -558,7 +558,7 @@ class Plot(IDManagerMixin):
cv.check_type('background', background, Iterable)
# Get a background (R,G,B) tuple to apply in alpha compositing
if isinstance(background, string_types):
if isinstance(background, str):
if background.lower() not in _SVG_COLORS:
raise ValueError("'{}' is not a valid color.".format(background))
background = _SVG_COLORS[background.lower()]
@ -570,7 +570,7 @@ class Plot(IDManagerMixin):
# other than those the user wishes to highlight
for domain, color in self.colors.items():
if domain not in domains:
if isinstance(color, string_types):
if isinstance(color, str):
color = _SVG_COLORS[color.lower()]
r, g, b = color
r = int(((1-alpha) * background[0]) + (alpha * r))
@ -610,7 +610,7 @@ class Plot(IDManagerMixin):
if self._background is not None:
subelement = ET.SubElement(element, "background")
color = self._background
if isinstance(color, string_types):
if isinstance(color, str):
color = _SVG_COLORS[color.lower()]
subelement.text = ' '.join(str(x) for x in color)
@ -619,7 +619,7 @@ class Plot(IDManagerMixin):
key=lambda x: x[0].id):
subelement = ET.SubElement(element, "color")
subelement.set("id", str(domain.id))
if isinstance(color, string_types):
if isinstance(color, str):
color = _SVG_COLORS[color.lower()]
subelement.set("rgb", ' '.join(str(x) for x in color))
@ -629,7 +629,7 @@ class Plot(IDManagerMixin):
str(d.id) for d in self._mask_components))
color = self._mask_background
if color is not None:
if isinstance(color, string_types):
if isinstance(color, str):
color = _SVG_COLORS[color.lower()]
subelement.set("background", ' '.join(
str(x) for x in color))
@ -651,6 +651,49 @@ class Plot(IDManagerMixin):
return element
def to_ipython_image(self, openmc_exec='openmc', cwd='.',
convert_exec='convert'):
"""Render plot as an image
This method runs OpenMC in plotting mode to produce a bitmap image which
is then converted to a .png file and loaded in as an
:class:`IPython.display.Image` object. As such, it requires that your
model geometry, materials, and settings have already been exported to
XML.
Parameters
----------
openmc_exec : str
Path to OpenMC executable
cwd : str, optional
Path to working directory to run in
convert_exec : str, optional
Command that can convert PPM files into PNG files
Returns
-------
IPython.display.Image
Image generated
"""
from IPython.display import Image
# Create plots.xml
Plots([self]).export_to_xml()
# Run OpenMC in geometry plotting mode
openmc.plot_geometry(False, openmc_exec, cwd)
# Convert to .png
if self.filename is not None:
ppm_file = '{}.ppm'.format(self.filename)
else:
ppm_file = 'plot_{}.ppm'.format(self.id)
png_file = ppm_file.replace('.ppm', '.png')
subprocess.check_call([convert_exec, ppm_file, png_file])
return Image(png_file)
class Plots(cv.CheckedList):
"""Collection of Plots used for an OpenMC simulation.
@ -674,28 +717,11 @@ class Plots(cv.CheckedList):
"""
def __init__(self, plots=None):
super(Plots, self).__init__(Plot, 'plots collection')
super().__init__(Plot, 'plots collection')
self._plots_file = ET.Element("plots")
if plots is not None:
self += plots
def add_plot(self, plot):
"""Add a plot to the file.
.. deprecated:: 0.8
Use :meth:`Plots.append` instead.
Parameters
----------
plot : openmc.Plot
Plot to add
"""
warnings.warn("Plots.add_plot(...) has been deprecated and may be "
"removed in a future version. Use Plots.append(...) "
"instead.", DeprecationWarning)
self.append(plot)
def append(self, plot):
"""Append plot to collection
@ -705,7 +731,7 @@ class Plots(cv.CheckedList):
Plot to append
"""
super(Plots, self).append(plot)
super().append(plot)
def insert(self, index, plot):
"""Insert plot before index
@ -718,24 +744,7 @@ class Plots(cv.CheckedList):
Plot to insert
"""
super(Plots, self).insert(index, plot)
def remove_plot(self, plot):
"""Remove a plot from the file.
.. deprecated:: 0.8
Use :meth:`Plots.remove` instead.
Parameters
----------
plot : openmc.Plot
Plot to remove
"""
warnings.warn("Plots.remove_plot(...) has been deprecated and may be "
"removed in a future version. Use Plots.remove(...) "
"instead.", DeprecationWarning)
self.remove(plot)
super().insert(index, plot)
def colorize(self, geometry, seed=1):
"""Generate a consistent color scheme for each domain in each plot.

10
openmc/plotter.py
View file

@ -2,8 +2,6 @@ from numbers import Integral, Real
from itertools import chain
import string
from six import string_types
import matplotlib.pyplot as plt
import numpy as np
import openmc.checkvalue as cv
@ -126,6 +124,8 @@ def plot_xs(this, types, divisor_types=None, temperature=294., data_type=None,
generated.
"""
import matplotlib.pyplot as plt
cv.check_type("plot_CE", plot_CE, bool)
if data_type is None:
@ -137,7 +137,7 @@ def plot_xs(this, types, divisor_types=None, temperature=294., data_type=None,
data_type = 'material'
elif isinstance(this, openmc.Macroscopic):
data_type = 'macroscopic'
elif isinstance(this, string_types):
elif isinstance(this, str):
if this[-1] in string.digits:
data_type = 'nuclide'
else:
@ -275,7 +275,7 @@ def calculate_cexs(this, data_type, types, temperature=294., sab_name=None,
# Check types
cv.check_type('temperature', temperature, Real)
if sab_name:
cv.check_type('sab_name', sab_name, string_types)
cv.check_type('sab_name', sab_name, str)
if enrichment:
cv.check_type('enrichment', enrichment, Real)
@ -648,7 +648,7 @@ def calculate_mgxs(this, data_type, types, orders=None, temperature=294.,
cv.check_type('temperature', temperature, Real)
if enrichment:
cv.check_type('enrichment', enrichment, Real)
cv.check_iterable_type('types', types, string_types)
cv.check_iterable_type('types', types, str)
cv.check_type("cross_sections", cross_sections, str)
library = openmc.MGXSLibrary.from_hdf5(cross_sections)

15
openmc/region.py
View file

@ -1,15 +1,14 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable, OrderedDict, MutableSequence
from collections import OrderedDict
from collections.abc import Iterable, MutableSequence
from copy import deepcopy
from six import add_metaclass
import numpy as np
from openmc.checkvalue import check_type
@add_metaclass(ABCMeta)
class Region(object):
class Region(metaclass=ABCMeta):
"""Region of space that can be assigned to a cell.
Region is an abstract base class that is inherited by
@ -39,10 +38,8 @@ class Region(object):
def __eq__(self, other):
if not isinstance(other, type(self)):
return False
elif str(self) != str(other):
return False
else:
return True
return str(self) == str(other)
def __ne__(self, other):
return not self == other
@ -463,7 +460,7 @@ class Union(Region, MutableSequence):
if memo is None:
memo = {}
clone = copy.deepcopy(self)
clone = deepcopy(self)
clone[:] = [n.clone(memo) for n in self]
return clone
@ -584,6 +581,6 @@ class Complement(Region):
if memo is None:
memo = {}
clone = copy.deepcopy(self)
clone = deepcopy(self)
clone.node = self.node.clone(memo)
return clone

6
openmc/search.py
View file

@ -1,4 +1,4 @@
from collections import Callable
from collections.abc import Callable
from numbers import Real
import scipy.optimize as sopt
@ -60,9 +60,9 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
if print_iterations:
text = 'Iteration: {}; Guess of {:.2e} produced a keff of ' + \
'{:1.5f} +/- {:1.5f}'
print(text.format(len(guesses), guess, keff[0], keff[1]))
print(text.format(len(guesses), guess, keff.n, keff.s))
return (keff[0] - target)
return keff.n - target
def search_for_keff(model_builder, initial_guess=None, target=1.0,

214
openmc/settings.py
View file

@ -1,10 +1,9 @@
from collections import Iterable, MutableSequence, Mapping
from collections.abc import Iterable, MutableSequence, Mapping
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
from openmc.clean_xml import clean_xml_indentation
@ -29,13 +28,6 @@ class Settings(object):
deviation.
create_fission_neutrons : bool
Indicate whether fission neutrons should be created or not.
cross_sections : str
Indicates the path to an XML cross section listing file (usually named
cross_sections.xml). If it is not set, the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable will be used for
continuous-energy calculations and
:envvar:`OPENMC_MG_CROSS_SECTIONS` will be used for multi-group
calculations to find the path to the XML cross section file.
cutoff : dict
Dictionary defining weight cutoff and energy cutoff. The dictionary may
have three keys, 'weight', 'weight_avg' and 'energy'. Value for 'weight'
@ -60,13 +52,10 @@ class Settings(object):
type are 'variance', 'std_dev', and 'rel_err'. The threshold value
should be a float indicating the variance, standard deviation, or
relative error used.
log_grid_bins : int
Number of bins for logarithmic energy grid search
max_order : None or int
Maximum scattering order to apply globally when in multi-group mode.
multipole_library : str
Indicates the path to a directory containing a windowed multipole
cross section library. If it is not set, the
:envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. A
multipole library is optional.
no_reduce : bool
Indicate that all user-defined and global tallies should not be reduced
across processes in a parallel calculation.
@ -220,19 +209,12 @@ class Settings(object):
# Uniform fission source subelement
self._ufs_mesh = None
# Domain decomposition subelement
self._dd_mesh_dimension = None
self._dd_mesh_lower_left = None
self._dd_mesh_upper_right = None
self._dd_nodemap = None
self._dd_allow_leakage = False
self._dd_count_interactions = False
self._resonance_scattering = {}
self._volume_calculations = cv.CheckedList(
VolumeCalculation, 'volume calculations')
self._create_fission_neutrons = None
self._log_grid_bins = None
@property
def run_mode(self):
@ -274,14 +256,6 @@ class Settings(object):
def confidence_intervals(self):
return self._confidence_intervals
@property
def cross_sections(self):
return self._cross_sections
@property
def multipole_library(self):
return self._multipole_library
@property
def ptables(self):
return self._ptables
@ -362,30 +336,6 @@ class Settings(object):
def ufs_mesh(self):
return self._ufs_mesh
@property
def dd_mesh_dimension(self):
return self._dd_mesh_dimension
@property
def dd_mesh_lower_left(self):
return self._dd_mesh_lower_left
@property
def dd_mesh_upper_right(self):
return self._dd_mesh_upper_right
@property
def dd_nodemap(self):
return self._dd_nodemap
@property
def dd_allow_leakage(self):
return self._dd_allow_leakage
@property
def dd_count_interactions(self):
return self._dd_count_interactions
@property
def resonance_scattering(self):
return self._resonance_scattering
@ -398,6 +348,10 @@ class Settings(object):
def create_fission_neutrons(self):
return self._create_fission_neutrons
@property
def log_grid_bins(self):
return self._log_grid_bins
@run_mode.setter
def run_mode(self, run_mode):
cv.check_value('run mode', run_mode, _RUN_MODES)
@ -484,7 +438,7 @@ class Settings(object):
if key in ('summary', 'tallies'):
cv.check_type("output['{}']".format(key), value, bool)
else:
cv.check_type("output['path']", value, string_types)
cv.check_type("output['path']", value, str)
self._output = output
@verbosity.setter
@ -531,23 +485,6 @@ class Settings(object):
cv.check_type('confidence interval', confidence_intervals, bool)
self._confidence_intervals = confidence_intervals
@cross_sections.setter
def cross_sections(self, cross_sections):
warnings.warn('Settings.cross_sections has been deprecated and will be '
'removed in a future version. Materials.cross_sections '
'should defined instead.', DeprecationWarning)
cv.check_type('cross sections', cross_sections, string_types)
self._cross_sections = cross_sections
@multipole_library.setter
def multipole_library(self, multipole_library):
warnings.warn('Settings.multipole_library has been deprecated and will '
'be removed in a future version. '
'Materials.multipole_library should defined instead.',
DeprecationWarning)
cv.check_type('multipole library', multipole_library, string_types)
self._multipole_library = multipole_library
@ptables.setter
def ptables(self, ptables):
cv.check_type('probability tables', ptables, bool)
@ -696,85 +633,6 @@ class Settings(object):
cv.check_length('UFS mesh upper-right corner', ufs_mesh.upper_right, 3)
self._ufs_mesh = ufs_mesh
@dd_mesh_dimension.setter
def dd_mesh_dimension(self, dimension):
# TODO: remove this when domain decomposition is merged
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
cv.check_type('DD mesh dimension', dimension, Iterable, Integral)
cv.check_length('DD mesh dimension', dimension, 3)
self._dd_mesh_dimension = dimension
@dd_mesh_lower_left.setter
def dd_mesh_lower_left(self, lower_left):
# TODO: remove this when domain decomposition is merged
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
cv.check_type('DD mesh lower left corner', lower_left, Iterable, Real)
cv.check_length('DD mesh lower left corner', lower_left, 3)
self._dd_mesh_lower_left = lower_left
@dd_mesh_upper_right.setter
def dd_mesh_upper_right(self, upper_right):
# TODO: remove this when domain decomposition is merged
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
cv.check_type('DD mesh upper right corner', upper_right, Iterable, Real)
cv.check_length('DD mesh upper right corner', upper_right, 3)
self._dd_mesh_upper_right = upper_right
@dd_nodemap.setter
def dd_nodemap(self, nodemap):
# TODO: remove this when domain decomposition is merged
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
cv.check_type('DD nodemap', nodemap, Iterable)
nodemap = np.array(nodemap).flatten()
if self._dd_mesh_dimension is None:
msg = 'Must set DD mesh dimension before setting the nodemap'
raise ValueError(msg)
else:
len_nodemap = np.prod(self._dd_mesh_dimension)
if len(nodemap) < len_nodemap or len(nodemap) > len_nodemap:
msg = 'Unable to set DD nodemap with length "{0}" which ' \
'does not have the same dimensionality as the domain ' \
'mesh'.format(len(nodemap))
raise ValueError(msg)
self._dd_nodemap = nodemap
@dd_allow_leakage.setter
def dd_allow_leakage(self, allow):
# TODO: remove this when domain decomposition is merged
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
cv.check_type('DD allow leakage', allow, bool)
self._dd_allow_leakage = allow
@dd_count_interactions.setter
def dd_count_interactions(self, interactions):
# TODO: remove this when domain decomposition is merged
warnings.warn('This feature is not yet implemented in a release '
'version of openmc')
cv.check_type('DD count interactions', interactions, bool)
self._dd_count_interactions = interactions
@resonance_scattering.setter
def resonance_scattering(self, res):
cv.check_type('resonance scattering settings', res, Mapping)
@ -796,7 +654,7 @@ class Settings(object):
cv.check_greater_than(name, value, 0)
elif key == 'nuclides':
cv.check_type('resonance scattering nuclides', value,
Iterable, string_types)
Iterable, str)
self._resonance_scattering = res
@volume_calculations.setter
@ -812,6 +670,12 @@ class Settings(object):
create_fission_neutrons, bool)
self._create_fission_neutrons = create_fission_neutrons
@log_grid_bins.setter
def log_grid_bins(self, log_grid_bins):
cv.check_type('log grid bins', log_grid_bins, Real)
cv.check_greater_than('log grid bins', log_grid_bins, 0)
self._log_grid_bins = log_grid_bins
def _create_run_mode_subelement(self, root):
elem = ET.SubElement(root, "run_mode")
elem.text = self._run_mode
@ -913,16 +777,6 @@ class Settings(object):
element = ET.SubElement(root, "confidence_intervals")
element.text = str(self._confidence_intervals).lower()
def _create_cross_sections_subelement(self, root):
if self._cross_sections is not None:
element = ET.SubElement(root, "cross_sections")
element.text = str(self._cross_sections)
def _create_multipole_library_subelement(self, root):
if self._multipole_library is not None:
element = ET.SubElement(root, "multipole_library")
element.text = str(self._multipole_library)
def _create_ptables_subelement(self, root):
if self._ptables is not None:
element = ET.SubElement(root, "ptables")
@ -1033,33 +887,6 @@ class Settings(object):
subelement = ET.SubElement(root, "ufs_mesh")
subelement.text = str(self.ufs_mesh.id)
def _create_dd_subelement(self, root):
if self._dd_mesh_lower_left is not None and \
self._dd_mesh_upper_right is not None and \
self._dd_mesh_dimension is not None:
element = ET.SubElement(root, "domain_decomposition")
subelement = ET.SubElement(element, "mesh")
subsubelement = ET.SubElement(subelement, "dimension")
subsubelement.text = ' '.join(map(str, self._dd_mesh_dimension))
subsubelement = ET.SubElement(subelement, "lower_left")
subsubelement.text = ' '.join(map(str, self._dd_mesh_lower_left))
subsubelement = ET.SubElement(subelement, "upper_right")
subsubelement.text = ' '.join(map(str, self._dd_mesh_upper_right))
if self._dd_nodemap is not None:
subelement = ET.SubElement(element, "nodemap")
subelement.text = ' '.join(map(str, self._dd_nodemap))
subelement = ET.SubElement(element, "allow_leakage")
subelement.text = str(self._dd_allow_leakage).lower()
subelement = ET.SubElement(element, "count_interactions")
subelement.text = str(self._dd_count_interactions).lower()
def _create_resonance_scattering_subelement(self, root):
res = self.resonance_scattering
if res:
@ -1085,6 +912,11 @@ class Settings(object):
elem = ET.SubElement(root, "create_fission_neutrons")
elem.text = str(self._create_fission_neutrons).lower()
def _create_log_grid_bins_subelement(self, root):
if self._log_grid_bins is not None:
elem = ET.SubElement(root, "log_grid_bins")
elem.text = str(self._log_grid_bins)
def export_to_xml(self, path='settings.xml'):
"""Export simulation settings to an XML file.
@ -1109,8 +941,6 @@ class Settings(object):
self._create_statepoint_subelement(root_element)
self._create_sourcepoint_subelement(root_element)
self._create_confidence_intervals(root_element)
self._create_cross_sections_subelement(root_element)
self._create_multipole_library_subelement(root_element)
self._create_energy_mode_subelement(root_element)
self._create_max_order_subelement(root_element)
self._create_ptables_subelement(root_element)
@ -1128,10 +958,10 @@ class Settings(object):
self._create_trace_subelement(root_element)
self._create_track_subelement(root_element)
self._create_ufs_mesh_subelement(root_element)
self._create_dd_subelement(root_element)
self._create_resonance_scattering_subelement(root_element)
self._create_volume_calcs_subelement(root_element)
self._create_create_fission_neutrons_subelement(root_element)
self._create_log_grid_bins_subelement(root_element)
# Clean the indentation in the file to be user-readable
clean_xml_indentation(root_element)

4
openmc/source.py
View file

@ -2,8 +2,6 @@ from numbers import Real
import sys
from xml.etree import ElementTree as ET
from six import string_types
from openmc.stats.univariate import Univariate
from openmc.stats.multivariate import UnitSphere, Spatial
import openmc.checkvalue as cv
@ -78,7 +76,7 @@ class Source(object):
@file.setter
def file(self, filename):
cv.check_type('source file', filename, string_types)
cv.check_type('source file', filename, str)
self._file = filename
@space.setter

17
openmc/statepoint.py
View file

@ -1,3 +1,4 @@
from datetime import datetime
import re
import os
import warnings
@ -5,6 +6,7 @@ import glob
import numpy as np
import h5py
from uncertainties import ufloat
import openmc
import openmc.checkvalue as cv
@ -47,8 +49,8 @@ class StatePoint(object):
CMFD fission source distribution over all mesh cells and energy groups.
current_batch : int
Number of batches simulated
date_and_time : str
Date and time when simulation began
date_and_time : datetime.datetime
Date and time at which statepoint was written
entropy : numpy.ndarray
Shannon entropy of fission source at each batch
filters : dict
@ -59,8 +61,8 @@ class StatePoint(object):
global_tallies : numpy.ndarray of compound datatype
Global tallies for k-effective estimates and leakage. The compound
datatype has fields 'name', 'sum', 'sum_sq', 'mean', and 'std_dev'.
k_combined : list
Combined estimator for k-effective and its uncertainty
k_combined : uncertainties.UFloat
Combined estimator for k-effective
k_col_abs : float
Cross-product of collision and absorption estimates of k-effective
k_col_tra : float
@ -187,7 +189,8 @@ class StatePoint(object):
@property
def date_and_time(self):
return self._f.attrs['date_and_time'].decode()
s = self._f.attrs['date_and_time'].decode()
return datetime.strptime(s, '%Y-%m-%d %H:%M:%S')
@property
def entropy(self):
@ -255,7 +258,7 @@ class StatePoint(object):
@property
def k_combined(self):
if self.run_mode == 'eigenvalue':
return self._f['k_combined'].value
return ufloat(*self._f['k_combined'].value)
else:
return None
@ -457,7 +460,7 @@ class StatePoint(object):
@property
def version(self):
return tuple(self._f.attrs['version'])
return tuple(self._f.attrs['openmc_version'])
@property
def summary(self):

21
openmc/stats/multivariate.py
View file

@ -1,19 +1,17 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
from collections.abc import Iterable
from math import pi
from numbers import Real
import sys
from xml.etree import ElementTree as ET
from six import add_metaclass
import numpy as np
import openmc.checkvalue as cv
from openmc.stats.univariate import Univariate, Uniform
@add_metaclass(ABCMeta)
class UnitSphere(object):
class UnitSphere(metaclass=ABCMeta):
"""Distribution of points on the unit sphere.
This abstract class is used for angular distributions, since a direction is
@ -77,7 +75,7 @@ class PolarAzimuthal(UnitSphere):
"""
def __init__(self, mu=None, phi=None, reference_uvw=[0., 0., 1.]):
super(PolarAzimuthal, self).__init__(reference_uvw)
super().__init__(reference_uvw)
if mu is not None:
self.mu = mu
else:
@ -130,7 +128,7 @@ class Isotropic(UnitSphere):
"""
def __init__(self):
super(Isotropic, self).__init__()
super().__init__()
def to_xml_element(self):
"""Return XML representation of the isotropic distribution
@ -163,7 +161,7 @@ class Monodirectional(UnitSphere):
def __init__(self, reference_uvw=[1., 0., 0.]):
super(Monodirectional, self).__init__(reference_uvw)
super().__init__(reference_uvw)
def to_xml_element(self):
"""Return XML representation of the monodirectional distribution
@ -181,8 +179,7 @@ class Monodirectional(UnitSphere):
return element
@add_metaclass(ABCMeta)
class Spatial(object):
class Spatial(metaclass=ABCMeta):
"""Distribution of locations in three-dimensional Euclidean space.
Classes derived from this abstract class can be used for spatial
@ -225,7 +222,7 @@ class CartesianIndependent(Spatial):
def __init__(self, x, y, z):
super(CartesianIndependent, self).__init__()
super().__init__()
self.x = x
self.y = y
self.z = z
@ -301,7 +298,7 @@ class Box(Spatial):
def __init__(self, lower_left, upper_right, only_fissionable=False):
super(Box, self).__init__()
super().__init__()
self.lower_left = lower_left
self.upper_right = upper_right
self.only_fissionable = only_fissionable
@ -374,7 +371,7 @@ class Point(Spatial):
"""
def __init__(self, xyz=(0., 0., 0.)):
super(Point, self).__init__()
super().__init__()
self.xyz = xyz
@property

37
openmc/stats/univariate.py
View file

@ -1,10 +1,9 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
from collections.abc import Iterable
from numbers import Real
import sys
from xml.etree import ElementTree as ET
from six import add_metaclass
import numpy as np
import openmc.checkvalue as cv
@ -15,8 +14,7 @@ _INTERPOLATION_SCHEMES = ['histogram', 'linear-linear', 'linear-log',
'log-linear', 'log-log']
@add_metaclass(ABCMeta)
class Univariate(EqualityMixin):
class Univariate(EqualityMixin, metaclass=ABCMeta):
"""Probability distribution of a single random variable.
The Univariate class is an abstract class that can be derived to implement a
@ -59,7 +57,7 @@ class Discrete(Univariate):
"""
def __init__(self, x, p):
super(Discrete, self).__init__()
super().__init__()
self.x = x
self.p = p
@ -133,7 +131,7 @@ class Uniform(Univariate):
"""
def __init__(self, a=0.0, b=1.0):
super(Uniform, self).__init__()
super().__init__()
self.a = a
self.b = b
@ -194,17 +192,17 @@ class Maxwell(Univariate):
Parameters
----------
theta : float
Effective temperature for distribution
Effective temperature for distribution in eV
Attributes
----------
theta : float
Effective temperature for distribution
Effective temperature for distribution in eV
"""
def __init__(self, theta):
super(Maxwell, self).__init__()
super().__init__()
self.theta = theta
def __len__(self):
@ -250,21 +248,21 @@ class Watt(Univariate):
Parameters
----------
a : float
First parameter of distribution
First parameter of distribution in units of eV
b : float
Second parameter of distribution
Second parameter of distribution in units of 1/eV
Attributes
----------
a : float
First parameter of distribution
First parameter of distribution in units of eV
b : float
Second parameter of distribution
Second parameter of distribution in units of 1/eV
"""
def __init__(self, a=0.988e6, b=2.249e-6):
super(Watt, self).__init__()
super().__init__()
self.a = a
self.b = b
@ -344,7 +342,7 @@ class Tabular(Univariate):
def __init__(self, x, p, interpolation='linear-linear',
ignore_negative=False):
super(Tabular, self).__init__()
super().__init__()
self._ignore_negative = ignore_negative
self.x = x
self.p = p
@ -444,10 +442,9 @@ class Legendre(Univariate):
def coefficients(self, coefficients):
cv.check_type('Legendre expansion coefficients', coefficients,
Iterable, Real)
for l in range(len(coefficients)):
coefficients[l] *= (2.*l + 1.)/2.
self._legendre_polynomial = np.polynomial.legendre.Legendre(
coefficients)
l = np.arange(len(coefficients))
coeffs = (2.*l + 1.)/2. * np.array(coefficients)
self._legendre_polynomial = np.polynomial.Legendre(coeffs)
def to_xml_element(self, element_name):
raise NotImplementedError
@ -473,7 +470,7 @@ class Mixture(Univariate):
"""
def __init__(self, probability, distribution):
super(Mixture, self).__init__()
super().__init__()
self.probability = probability
self.distribution = distribution

2
openmc/summary.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
import re
import warnings

44
openmc/surface.py
View file

@ -1,4 +1,3 @@
from __future__ import division
from abc import ABCMeta
from collections import OrderedDict
from copy import deepcopy
@ -6,7 +5,6 @@ from functools import partial
from numbers import Real, Integral
from xml.etree import ElementTree as ET
from six import add_metaclass, string_types
import numpy as np
from openmc.checkvalue import check_type, check_value
@ -115,14 +113,14 @@ class Surface(IDManagerMixin):
@name.setter
def name(self, name):
if name is not None:
check_type('surface name', name, string_types)
check_type('surface name', name, str)
self._name = name
else:
self._name = ''
@boundary_type.setter
def boundary_type(self, boundary_type):
check_type('boundary type', boundary_type, string_types)
check_type('boundary type', boundary_type, str)
check_value('boundary type', boundary_type, _BOUNDARY_TYPES)
self._boundary_type = boundary_type
@ -328,7 +326,7 @@ class Plane(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
A=1., B=0., C=0., D=0., name=''):
super(Plane, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._type = 'plane'
self._coeff_keys = ['A', 'B', 'C', 'D']
@ -412,7 +410,7 @@ class Plane(Surface):
XML element containing source data
"""
element = super(Plane, self).to_xml_element()
element = super().to_xml_element()
# Add periodic surface pair information
if self.boundary_type == 'periodic':
@ -462,7 +460,7 @@ class XPlane(Plane):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., name=''):
super(XPlane, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._type = 'x-plane'
self._coeff_keys = ['x0']
@ -568,7 +566,7 @@ class YPlane(Plane):
def __init__(self, surface_id=None, boundary_type='transmission',
y0=0., name=''):
# Initialize YPlane class attributes
super(YPlane, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._type = 'y-plane'
self._coeff_keys = ['y0']
@ -674,7 +672,7 @@ class ZPlane(Plane):
def __init__(self, surface_id=None, boundary_type='transmission',
z0=0., name=''):
# Initialize ZPlane class attributes
super(ZPlane, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._type = 'z-plane'
self._coeff_keys = ['z0']
@ -738,8 +736,7 @@ class ZPlane(Plane):
return point[2] - self.z0
@add_metaclass(ABCMeta)
class Cylinder(Surface):
class Cylinder(Surface, metaclass=ABCMeta):
"""A cylinder whose length is parallel to the x-, y-, or z-axis.
Parameters
@ -776,7 +773,7 @@ class Cylinder(Surface):
"""
def __init__(self, surface_id=None, boundary_type='transmission',
R=1., name=''):
super(Cylinder, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['R']
self.r = R
@ -836,7 +833,7 @@ class XCylinder(Cylinder):
def __init__(self, surface_id=None, boundary_type='transmission',
y0=0., z0=0., R=1., name=''):
super(XCylinder, self).__init__(surface_id, boundary_type, R, name=name)
super().__init__(surface_id, boundary_type, R, name=name)
self._type = 'x-cylinder'
self._coeff_keys = ['y0', 'z0', 'R']
@ -958,7 +955,7 @@ class YCylinder(Cylinder):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., z0=0., R=1., name=''):
super(YCylinder, self).__init__(surface_id, boundary_type, R, name=name)
super().__init__(surface_id, boundary_type, R, name=name)
self._type = 'y-cylinder'
self._coeff_keys = ['x0', 'z0', 'R']
@ -1080,7 +1077,7 @@ class ZCylinder(Cylinder):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., R=1., name=''):
super(ZCylinder, self).__init__(surface_id, boundary_type, R, name=name)
super().__init__(surface_id, boundary_type, R, name=name)
self._type = 'z-cylinder'
self._coeff_keys = ['x0', 'y0', 'R']
@ -1206,7 +1203,7 @@ class Sphere(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R=1., name=''):
super(Sphere, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._type = 'sphere'
self._coeff_keys = ['x0', 'y0', 'z0', 'R']
@ -1305,8 +1302,7 @@ class Sphere(Surface):
return x**2 + y**2 + z**2 - self.r**2
@add_metaclass(ABCMeta)
class Cone(Surface):
class Cone(Surface, metaclass=ABCMeta):
"""A conical surface parallel to the x-, y-, or z-axis.
Parameters
@ -1354,7 +1350,7 @@ class Cone(Surface):
"""
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super(Cone, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
self.x0 = x0
@ -1376,7 +1372,7 @@ class Cone(Surface):
@property
def r2(self):
return self.coefficients['r2']
return self.coefficients['R2']
@x0.setter
def x0(self, x0):
@ -1449,7 +1445,7 @@ class XCone(Cone):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super(XCone, self).__init__(surface_id, boundary_type, x0, y0,
super().__init__(surface_id, boundary_type, x0, y0,
z0, R2, name=name)
self._type = 'x-cone'
@ -1525,7 +1521,7 @@ class YCone(Cone):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super(YCone, self).__init__(surface_id, boundary_type, x0, y0, z0,
super().__init__(surface_id, boundary_type, x0, y0, z0,
R2, name=name)
self._type = 'y-cone'
@ -1601,7 +1597,7 @@ class ZCone(Cone):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super(ZCone, self).__init__(surface_id, boundary_type, x0, y0, z0,
super().__init__(surface_id, boundary_type, x0, y0, z0,
R2, name=name)
self._type = 'z-cone'
@ -1666,7 +1662,7 @@ class Quadric(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
a=0., b=0., c=0., d=0., e=0., f=0., g=0.,
h=0., j=0., k=0., name=''):
super(Quadric, self).__init__(surface_id, boundary_type, name=name)
super().__init__(surface_id, boundary_type, name=name)
self._type = 'quadric'
self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k']

Some files were not shown because too many files have changed in this diff Show more