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

Merge branch 'develop' into cmfd-added-funcs

Browse source
This commit is contained in:
Shikhar Kumar 2019-04-07 22:37:45 -04:00
commit 196d24f0db
166 changed files with 4802 additions and 3166 deletions
Download patch file Download diff file Expand all files Collapse all files

33
CMakeLists.txt
View file

@ -1,5 +1,5 @@
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc CXX)
project(openmc C CXX)
# Setup output directories
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
@ -20,16 +20,13 @@ option(optimize "Turn on all compiler optimization flags" OFF)
option(coverage "Compile with coverage analysis flags" OFF)
option(dagmc "Enable support for DAGMC (CAD) geometry" OFF)
# Maximum number of nested coordinates levels
set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
#===============================================================================
# MPI for distributed-memory parallelism
#===============================================================================
set(MPI_ENABLED FALSE)
if($ENV{CXX} MATCHES "(mpi[^/]*|CC)$")
message("-- Detected MPI wrapper: $ENV{CXX}")
message(STATUS "Detected MPI wrapper: $ENV{CXX}")
set(MPI_ENABLED TRUE)
endif()
@ -38,9 +35,6 @@ endif()
#===============================================================================
if(dagmc)
find_package(DAGMC REQUIRED)
if(NOT DAGMC_FOUND)
message(FATAL_ERROR "Could not find DAGMC installation")
endif()
link_directories(${DAGMC_LIBRARY_DIRS})
endif()
@ -66,15 +60,12 @@ if(NOT DEFINED HDF5_PREFER_PARALLEL)
endif()
endif()
find_package(HDF5 COMPONENTS HL)
if(NOT HDF5_FOUND)
message(FATAL_ERROR "Could not find HDF5")
endif()
find_package(HDF5 REQUIRED COMPONENTS C HL)
if(HDF5_IS_PARALLEL)
if(NOT MPI_ENABLED)
message(FATAL_ERROR "Parallel HDF5 must be used with MPI.")
endif()
message("-- Using parallel HDF5")
message(STATUS "Using parallel HDF5")
endif()
#===============================================================================
@ -98,7 +89,7 @@ if(debug)
list(APPEND cxxflags -g -O0)
endif()
if(profile)
list(APPEND cxxflags -pg)
list(APPEND cxxflags -g -fno-omit-frame-pointer)
endif()
if(optimize)
list(REMOVE_ITEM cxxflags -O2)
@ -106,8 +97,8 @@ if(optimize)
endif()
# Show flags being used
message(STATUS "C++ flags: ${cxxflags}")
message(STATUS "Linker flags: ${ldflags}")
message(STATUS "OpenMC C++ flags: ${cxxflags}")
message(STATUS "OpenMC Linker flags: ${ldflags}")
#===============================================================================
# pugixml library
@ -149,7 +140,7 @@ set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
# the RPATH to be used when installing, but only if it's not a system directory
list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir)
if("${isSystemDir}" STREQUAL "-1")
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
endif()
#===============================================================================
@ -255,13 +246,11 @@ set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc)
target_include_directories(libopenmc
PUBLIC include
PRIVATE ${HDF5_INCLUDE_DIRS})
PUBLIC include ${HDF5_INCLUDE_DIRS})
# Set compile flags
target_compile_options(libopenmc PRIVATE ${cxxflags})
target_compile_definitions(libopenmc PRIVATE -DMAX_COORD=${maxcoord})
if (HDF5_IS_PARALLEL)
target_compile_definitions(libopenmc PRIVATE -DPHDF5)
endif()
@ -281,8 +270,8 @@ endif()
# target_link_libraries treats any arguments starting with - but not -l as
# linker flags. Thus, we can pass both linker flags and libraries together.
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml
faddeeva xtensor)
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}
pugixml faddeeva xtensor)
if(dagmc)
target_compile_definitions(libopenmc PRIVATE DAGMC)

399
cmake/Modules/FindHDF5.cmake
View file

@ -1,399 +0,0 @@
#.rst:
# FindHDF5
# --------
#
# Find HDF5, a library for reading and writing self describing array data.
#
#
#
# This module invokes the HDF5 wrapper compiler that should be installed
# alongside HDF5. Depending upon the HDF5 Configuration, the wrapper
# compiler is called either h5cc or h5pcc. If this succeeds, the module
# will then call the compiler with the -show argument to see what flags
# are used when compiling an HDF5 client application.
#
# The module will optionally accept the COMPONENTS argument. If no
# COMPONENTS are specified, then the find module will default to finding
# only the HDF5 C library. If one or more COMPONENTS are specified, the
# module will attempt to find the language bindings for the specified
# components. The only valid components are C, CXX, Fortran, HL, and
# Fortran_HL. If the COMPONENTS argument is not given, the module will
# attempt to find only the C bindings.
#
# On UNIX systems, this module will read the variable
# HDF5_USE_STATIC_LIBRARIES to determine whether or not to prefer a
# static link to a dynamic link for HDF5 and all of it's dependencies.
# To use this feature, make sure that the HDF5_USE_STATIC_LIBRARIES
# variable is set before the call to find_package.
#
# To provide the module with a hint about where to find your HDF5
# installation, you can set the environment variable HDF5_ROOT. The
# Find module will then look in this path when searching for HDF5
# executables, paths, and libraries.
#
# In addition to finding the includes and libraries required to compile
# an HDF5 client application, this module also makes an effort to find
# tools that come with the HDF5 distribution that may be useful for
# regression testing.
#
# This module will define the following variables:
#
# ::
#
# HDF5_INCLUDE_DIRS - Location of the hdf5 includes
# HDF5_INCLUDE_DIR - Location of the hdf5 includes (deprecated)
# HDF5_DEFINITIONS - Required compiler definitions for HDF5
# HDF5_C_LIBRARIES - Required libraries for the HDF5 C bindings.
# HDF5_CXX_LIBRARIES - Required libraries for the HDF5 C++ bindings
# HDF5_Fortran_LIBRARIES - Required libraries for the HDF5 Fortran bindings
# HDF5_HL_LIBRARIES - Required libraries for the HDF5 high level API
# HDF5_Fortran_HL_LIBRARIES - Required libraries for the high level Fortran
# bindings.
# HDF5_LIBRARIES - Required libraries for all requested bindings
# HDF5_FOUND - true if HDF5 was found on the system
# HDF5_VERSION - HDF5 version in format Major.Minor.Release
# HDF5_LIBRARY_DIRS - the full set of library directories
# HDF5_IS_PARALLEL - Whether or not HDF5 was found with parallel IO support
# HDF5_C_COMPILER_EXECUTABLE - the path to the HDF5 C wrapper compiler
# HDF5_CXX_COMPILER_EXECUTABLE - the path to the HDF5 C++ wrapper compiler
# HDF5_Fortran_COMPILER_EXECUTABLE - the path to the HDF5 Fortran wrapper compiler
# HDF5_DIFF_EXECUTABLE - the path to the HDF5 dataset comparison tool
#=============================================================================
# Copyright 2015 Axel Huebl, Helmholtz-Zentrum Dresden - Rossendorf
# Copyright 2009 Kitware, Inc.
#
# Distributed under the OSI-approved BSD License (the "License");
# see accompanying file Copyright.txt for details.
#
# This software is distributed WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the License for more information.
#=============================================================================
# (To distribute this file outside of CMake, substitute the full
# License text for the above reference.)
# This module is maintained by Will Dicharry <wdicharry@stellarscience.com>.
include(SelectLibraryConfigurations)
include(FindPackageHandleStandardArgs)
# List of the valid HDF5 components
set( HDF5_VALID_COMPONENTS
C
CXX
Fortran
HL
Fortran_HL
)
# Validate the list of find components.
if( NOT HDF5_FIND_COMPONENTS )
set( HDF5_LANGUAGE_BINDINGS "C" )
else()
# add the extra specified components, ensuring that they are valid.
foreach( component ${HDF5_FIND_COMPONENTS} )
list( FIND HDF5_VALID_COMPONENTS ${component} component_location )
if( ${component_location} EQUAL -1 )
message( FATAL_ERROR
"\"${component}\" is not a valid HDF5 component." )
else()
list( APPEND HDF5_LANGUAGE_BINDINGS ${component} )
endif()
endforeach()
endif()
# Determine whether to search for serial or parallel executable first
if(HDF5_PREFER_PARALLEL)
set(HDF5_C_COMPILER_NAMES h5pcc h5cc)
set(HDF5_CXX_COMPILER_NAMES h5pc++ h5c++)
set(HDF5_Fortran_COMPILER_NAMES h5pfc h5fc)
else()
set(HDF5_C_COMPILER_NAMES h5cc h5pcc)
set(HDF5_CXX_COMPILER_NAMES h5c++ h5pc++)
set(HDF5_Fortran_COMPILER_NAMES h5fc h5pfc)
endif()
# try to find the HDF5 wrapper compilers
find_program( HDF5_C_COMPILER_EXECUTABLE
NAMES ${HDF5_C_COMPILER_NAMES}
HINTS ENV HDF5_ROOT
PATH_SUFFIXES bin Bin
DOC "HDF5 Wrapper compiler. Used only to detect HDF5 compile flags." )
mark_as_advanced( HDF5_C_COMPILER_EXECUTABLE )
find_program( HDF5_CXX_COMPILER_EXECUTABLE
NAMES ${HDF5_CXX_COMPILER_NAMES}
HINTS ENV HDF5_ROOT
PATH_SUFFIXES bin Bin
DOC "HDF5 C++ Wrapper compiler. Used only to detect HDF5 compile flags." )
mark_as_advanced( HDF5_CXX_COMPILER_EXECUTABLE )
find_program( HDF5_Fortran_COMPILER_EXECUTABLE
NAMES ${HDF5_Fortran_COMPILER_NAMES}
HINTS ENV HDF5_ROOT
PATH_SUFFIXES bin Bin
DOC "HDF5 Fortran Wrapper compiler. Used only to detect HDF5 compile flags." )
mark_as_advanced( HDF5_Fortran_COMPILER_EXECUTABLE )
unset(HDF5_C_COMPILER_NAMES)
unset(HDF5_CXX_COMPILER_NAMES)
unset(HDF5_Fortran_COMPILER_NAMES)
find_program( HDF5_DIFF_EXECUTABLE
NAMES h5diff
HINTS ENV HDF5_ROOT
PATH_SUFFIXES bin Bin
DOC "HDF5 file differencing tool." )
mark_as_advanced( HDF5_DIFF_EXECUTABLE )
# Invoke the HDF5 wrapper compiler. The compiler return value is stored to the
# return_value argument, the text output is stored to the output variable.
macro( _HDF5_invoke_compiler language output return_value )
if( HDF5_${language}_COMPILER_EXECUTABLE )
exec_program( ${HDF5_${language}_COMPILER_EXECUTABLE}
ARGS -show
OUTPUT_VARIABLE ${output}
RETURN_VALUE ${return_value}
)
if( ${${return_value}} EQUAL 0 )
# do nothing
else()
message( STATUS
"Unable to determine HDF5 ${language} flags from HDF5 wrapper." )
endif()
endif()
endmacro()
# Parse a compile line for definitions, includes, library paths, and libraries.
macro( _HDF5_parse_compile_line
compile_line_var
include_paths
definitions
library_paths
libraries )
# Match the include paths
string( REGEX MATCHALL "-I([^\" ]+)" include_path_flags
"${${compile_line_var}}"
)
foreach( IPATH ${include_path_flags} )
string( REGEX REPLACE "^-I" "" IPATH ${IPATH} )
string( REPLACE "//" "/" IPATH ${IPATH} )
list( APPEND ${include_paths} ${IPATH} )
endforeach()
# Match the definitions
string( REGEX MATCHALL "-D[^ ]*" definition_flags "${${compile_line_var}}" )
foreach( DEF ${definition_flags} )
list( APPEND ${definitions} ${DEF} )
endforeach()
# Match the library paths
string( REGEX MATCHALL "-L([^\" ]+|\"[^\"]+\")" library_path_flags
"${${compile_line_var}}"
)
foreach( LPATH ${library_path_flags} )
string( REGEX REPLACE "^-L" "" LPATH ${LPATH} )
string( REPLACE "//" "/" LPATH ${LPATH} )
list( APPEND ${library_paths} ${LPATH} )
endforeach()
# now search for the library names specified in the compile line (match -l...)
# match only -l's preceded by a space or comma
# this is to exclude directory names like xxx-linux/
string( REGEX MATCHALL "[, ]-l([^\", ]+)" library_name_flags
"${${compile_line_var}}" )
# strip the -l from all of the library flags and add to the search list
foreach( LIB ${library_name_flags} )
string( REGEX REPLACE "^[, ]-l" "" LIB ${LIB} )
list( APPEND ${libraries} ${LIB} )
endforeach()
endmacro()
# Try to find HDF5 using an installed hdf5-config.cmake
if( NOT HDF5_FOUND )
find_package( HDF5 QUIET NO_MODULE )
if( HDF5_FOUND )
set( HDF5_INCLUDE_DIRS ${HDF5_INCLUDE_DIR} )
set( HDF5_LIBRARIES )
set( HDF5_C_TARGET hdf5 )
set( HDF5_CXX_TARGET hdf5_cpp )
set( HDF5_HL_TARGET hdf5_hl )
set( HDF5_Fortran_TARGET hdf5_fortran )
set( HDF5_Fortran_HL_TARGET hdf5_hl_fortran )
foreach( _component ${HDF5_LANGUAGE_BINDINGS} )
list( FIND HDF5_VALID_COMPONENTS ${_component} _component_location )
get_target_property( _comp_location ${HDF5_${_component}_TARGET} LOCATION )
if( _comp_location )
set( HDF5_${_component}_LIBRARY ${_comp_location} CACHE PATH
"HDF5 ${_component} library" )
mark_as_advanced( HDF5_${_component}_LIBRARY )
list( APPEND HDF5_LIBRARIES ${HDF5_${_component}_LIBRARY} )
endif()
endforeach()
endif()
endif()
if( NOT HDF5_FOUND )
_HDF5_invoke_compiler( C HDF5_C_COMPILE_LINE HDF5_C_RETURN_VALUE )
_HDF5_invoke_compiler( CXX HDF5_CXX_COMPILE_LINE HDF5_CXX_RETURN_VALUE )
_HDF5_invoke_compiler( Fortran HDF5_Fortran_COMPILE_LINE HDF5_Fortran_RETURN_VALUE )
set(HDF5_HL_COMPILE_LINE ${HDF5_C_COMPILE_LINE})
set(HDF5_Fortran_HL_COMPILE_LINE ${HDF5_Fortran_COMPILE_LINE})
# seed the initial lists of libraries to find with items we know we need
set( HDF5_C_LIBRARY_NAMES_INIT hdf5 )
set( HDF5_HL_LIBRARY_NAMES_INIT hdf5_hl ${HDF5_C_LIBRARY_NAMES_INIT} )
set( HDF5_CXX_LIBRARY_NAMES_INIT hdf5_cpp ${HDF5_C_LIBRARY_NAMES_INIT} )
set( HDF5_Fortran_LIBRARY_NAMES_INIT hdf5_fortran
${HDF5_C_LIBRARY_NAMES_INIT} )
set( HDF5_Fortran_HL_LIBRARY_NAMES_INIT hdf5hl_fortran hdf5_hl
${HDF5_Fortran_LIBRARY_NAMES_INIT} )
foreach( LANGUAGE ${HDF5_LANGUAGE_BINDINGS} )
if( HDF5_${LANGUAGE}_COMPILE_LINE )
_HDF5_parse_compile_line( HDF5_${LANGUAGE}_COMPILE_LINE
HDF5_${LANGUAGE}_INCLUDE_FLAGS
HDF5_${LANGUAGE}_DEFINITIONS
HDF5_${LANGUAGE}_LIBRARY_DIRS
HDF5_${LANGUAGE}_LIBRARY_NAMES
)
# take a guess that the includes may be in the 'include' sibling
# directory of a library directory.
foreach( dir ${HDF5_${LANGUAGE}_LIBRARY_DIRS} )
list( APPEND HDF5_${LANGUAGE}_INCLUDE_FLAGS ${dir}/../include )
endforeach()
endif()
# set the definitions for the language bindings.
list( APPEND HDF5_DEFINITIONS ${HDF5_${LANGUAGE}_DEFINITIONS} )
# find the HDF5 include directories
if(${LANGUAGE} MATCHES "Fortran")
set(HDF5_INCLUDE_FILENAME hdf5.mod)
else()
set(HDF5_INCLUDE_FILENAME hdf5.h)
endif()
find_path( HDF5_${LANGUAGE}_INCLUDE_DIR ${HDF5_INCLUDE_FILENAME}
HINTS
${HDF5_${LANGUAGE}_INCLUDE_FLAGS}
ENV
HDF5_ROOT
PATHS
$ENV{HOME}/.local/include
PATH_SUFFIXES
include
Include
)
mark_as_advanced( HDF5_${LANGUAGE}_INCLUDE_DIR )
list( APPEND HDF5_INCLUDE_DIRS ${HDF5_${LANGUAGE}_INCLUDE_DIR} )
# find the HDF5 libraries
foreach( LIB ${HDF5_${LANGUAGE}_LIBRARY_NAMES_INIT} )
if( UNIX AND HDF5_USE_STATIC_LIBRARIES )
# According to bug 1643 on the CMake bug tracker, this is the
# preferred method for searching for a static library.
# See http://www.cmake.org/Bug/view.php?id=1643. We search
# first for the full static library name, but fall back to a
# generic search on the name if the static search fails.
set( THIS_LIBRARY_SEARCH_DEBUG lib${LIB}d.a ${LIB}d )
set( THIS_LIBRARY_SEARCH_RELEASE lib${LIB}.a ${LIB} )
else()
set( THIS_LIBRARY_SEARCH_DEBUG ${LIB}d )
set( THIS_LIBRARY_SEARCH_RELEASE ${LIB} )
endif()
find_library( HDF5_${LIB}_LIBRARY_DEBUG
NAMES ${THIS_LIBRARY_SEARCH_DEBUG}
HINTS ${HDF5_${LANGUAGE}_LIBRARY_DIRS}
ENV HDF5_ROOT
PATH_SUFFIXES lib Lib )
find_library( HDF5_${LIB}_LIBRARY_RELEASE
NAMES ${THIS_LIBRARY_SEARCH_RELEASE}
HINTS ${HDF5_${LANGUAGE}_LIBRARY_DIRS}
ENV HDF5_ROOT
PATH_SUFFIXES lib Lib )
select_library_configurations( HDF5_${LIB} )
list(APPEND HDF5_${LANGUAGE}_LIBRARIES ${HDF5_${LIB}_LIBRARY})
endforeach()
list( APPEND HDF5_LIBRARY_DIRS ${HDF5_${LANGUAGE}_LIBRARY_DIRS} )
# When the wrapper lists a library with -l, e.g. -lz, simply use it as
# is. If find_library is called for these libraries, you end up with
# local libraries that will not be suitable when cross-compiling for the
# Intel Xeon Phi.
foreach(LIBNAME ${HDF5_${LANGUAGE}_LIBRARY_NAMES})
list(APPEND HDF5_${LANGUAGE}_LIBRARIES "-l${LIBNAME}")
endforeach()
# Append the libraries for this language binding to the list of all
# required libraries.
list(APPEND HDF5_LIBRARIES ${HDF5_${LANGUAGE}_LIBRARIES})
endforeach()
# We may have picked up some duplicates in various lists during the above
# process for the language bindings (both the C and C++ bindings depend on
# libz for example). Remove the duplicates. It appears that the default
# CMake behavior is to remove duplicates from the end of a list. However,
# for link lines, this is incorrect since unresolved symbols are searched
# for down the link line. Therefore, we reverse the list, remove the
# duplicates, and then reverse it again to get the duplicates removed from
# the beginning.
macro( _remove_duplicates_from_beginning _list_name )
list( REVERSE ${_list_name} )
list( REMOVE_DUPLICATES ${_list_name} )
list( REVERSE ${_list_name} )
endmacro()
if( HDF5_INCLUDE_DIRS )
_remove_duplicates_from_beginning( HDF5_INCLUDE_DIRS )
endif()
if( HDF5_LIBRARY_DIRS )
_remove_duplicates_from_beginning( HDF5_LIBRARY_DIRS )
endif()
# If the HDF5 include directory was found, open H5pubconf.h to determine if
# HDF5 was compiled with parallel IO support
set( HDF5_IS_PARALLEL FALSE )
set( HDF5_VERSION "" )
foreach( _dir IN LISTS HDF5_INCLUDE_DIRS )
foreach(_hdr "${_dir}/H5pubconf.h" "${_dir}/H5pubconf-64.h" "${_dir}/H5pubconf-32.h")
if( EXISTS "${_hdr}" )
file( STRINGS "${_hdr}"
HDF5_HAVE_PARALLEL_DEFINE
REGEX "HAVE_PARALLEL 1" )
if( HDF5_HAVE_PARALLEL_DEFINE )
set( HDF5_IS_PARALLEL TRUE )
endif()
unset(HDF5_HAVE_PARALLEL_DEFINE)
file( STRINGS "${_hdr}"
HDF5_VERSION_DEFINE
REGEX "^[ \t]*#[ \t]*define[ \t]+H5_VERSION[ \t]+" )
if( "${HDF5_VERSION_DEFINE}" MATCHES
"H5_VERSION[ \t]+\"([0-9]+\\.[0-9]+\\.[0-9]+).*\"" )
set( HDF5_VERSION "${CMAKE_MATCH_1}" )
endif()
unset(HDF5_VERSION_DEFINE)
endif()
endforeach()
endforeach()
set( HDF5_IS_PARALLEL ${HDF5_IS_PARALLEL} CACHE BOOL
"HDF5 library compiled with parallel IO support" )
mark_as_advanced( HDF5_IS_PARALLEL )
# For backwards compatibility we set HDF5_INCLUDE_DIR to the value of
# HDF5_INCLUDE_DIRS
if( HDF5_INCLUDE_DIRS )
set( HDF5_INCLUDE_DIR "${HDF5_INCLUDE_DIRS}" )
endif()
endif()
find_package_handle_standard_args( HDF5
REQUIRED_VARS HDF5_LIBRARIES HDF5_INCLUDE_DIRS
VERSION_VAR HDF5_VERSION
)

3
docs/source/conf.py
View file

@ -53,9 +53,10 @@ extensions = ['sphinx.ext.autodoc',
'sphinx.ext.autosummary',
'sphinx.ext.intersphinx',
'sphinx.ext.viewcode',
'sphinx.ext.imgconverter',
'sphinx_numfig',
'notebook_sphinxext']
if not on_rtd:
extensions.append('sphinx.ext.imgconverter')
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']

25
docs/source/devguide/styleguide.rst
View file

@ -12,6 +12,11 @@ adding new code in OpenMC.
C++
---
Indentation
-----------
Use two spaces per indentation level.
Miscellaneous
-------------
@ -126,6 +131,15 @@ single declaration to avoid confusion:
Curly braces
------------
For a class declaration, the opening brace should be on the same line that
lists the name of the class.
.. code-block:: C++
class Matrix {
...
};
For a function definition, the opening and closing braces should each be on
their own lines. This helps distinguish function code from the argument list.
If the entire function fits on one or two lines, then the braces can be on the
@ -210,11 +224,18 @@ Use of third-party Python packages should be limited to numpy_, scipy_,
matplotlib_, pandas_, and h5py_. Use of other third-party packages must be
implemented as optional dependencies rather than required dependencies.
Prefer pathlib_ when working with filesystem paths over functions in the os_
module or other standard-library modules. Functions that accept arguments that
represent a filesystem path should work with both strings and Path_ objects.
.. _C++ Core Guidelines: http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines
.. _PEP8: https://www.python.org/dev/peps/pep-0008/
.. _numpydoc: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt
.. _numpydoc: https://numpydoc.readthedocs.io/en/latest/format.html
.. _numpy: http://www.numpy.org/
.. _scipy: https://www.scipy.org/
.. _matplotlib: https://matplotlib.org/
.. _pandas: https://pandas.pydata.org/
.. _h5py: http://www.h5py.org/
.. _h5py: https://www.h5py.org/
.. _pathlib: https://docs.python.org/3/library/pathlib.html
.. _os: https://docs.python.org/3/library/os.html
.. _Path: https://docs.python.org/3/library/pathlib.html#pathlib.Path

13
docs/source/io_formats/nuclear_data.rst
View file

@ -133,11 +133,17 @@ Incident Photon Data
**/<element>/bremsstrahlung/**
:Attributes: - **I** (*double*) -- Mean excitation energy in [eV]
:Datasets: - **electron_energy** (*double[]*) -- Incident electron energy in [eV]
- **photon_energy** (*double[]*) -- Outgoing photon energy as
fraction of incident electron energy
- **dcs** (*double[][]*) -- Bremsstrahlung differential cross section
at each incident energy in [mb/eV]
- **ionization_energy** (*double[]*) -- Ionization potential of each
subshell in [eV]
- **num_electrons** (*int[]*) -- Number of electrons per subshell,
with conduction electrons indicated by a negative value
**/<element>/coherent/**
@ -176,13 +182,6 @@ Incident Photon Data
:Datasets: - **xs** (*double[]*) -- Total photoionization cross section in [b]
**/<element>/stopping_powers/**
:Datasets: - **I** (*double*) -- Mean excitation energy in [eV]
- **energy** (*double[]*) -- Energies in [eV]
- **s_collision** (*double[]*) -- Collision stopping power in [eV-cm\ :sup:`2`\ /g]
- **s_radiative** (*double[]*) -- Radiative stopping power in [eV-cm\ :sup:`2`\ /g]
**/<element>/subshells/**
:Attributes: - **designators** (*char[][]*) -- Designator for each shell, e.g. 'M2'

9
docs/source/io_formats/settings.rst
View file

@ -730,15 +730,6 @@ sections.
*Default*: 10 K
---------------------
``<threads>`` Element
---------------------
The ``<threads>`` element indicates the number of OpenMP threads to be used for
a simulation. It has no attributes and accepts a positive integer value.
*Default*: None (Determined by environment variable :envvar:`OMP_NUM_THREADS`)
.. _trace:
-------------------

168
docs/source/methods/photon_physics.rst
View file

@ -5,7 +5,7 @@ Photon Physics
==============
Photons, being neutral particles, behave much in the same manner as neutrons,
traveling in straight lines and experiencing occasional collisions which change
traveling in straight lines and experiencing occasional collisions that change
their energy and direction. Photons undergo four basic interactions as they pass
through matter: coherent (Rayleigh) scattering, incoherent (Compton) scattering,
photoelectric effect, and pair/triplet production. Photons with energy in the
@ -728,10 +728,18 @@ the cross section differential in energy loss. The total stopping power
power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to
bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`,
which refers to the energy loss due to inelastic collisions with bound
electrons in the material that result in ionization and excitation. To obtain
the radiative stopping power for positrons, the radiative stopping power for
electrons is multiplied by :eq:`positron-factor`. Currently, the collision
stopping power for electrons is also used for positrons.
electrons in the material that result in ionization and excitation. The
radiative stopping power for electrons is given by
.. math::
:label: radiative-stopping-power
S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa)
d\kappa.
To obtain the radiative stopping power for positrons,
:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`.
While the models for photon interactions with matter described above can safely
assume interactions occur with free atoms, sampling the target atom based on
@ -754,14 +762,97 @@ power is calculated using Bragg's additivity rule as
S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T),
where :math:`w_i` is the mass fraction of the :math:`i`-th element. The
collision stopping power, however, is a function of certain quantities such as
the mean excitation energy :math:`I` and the density effect correction
:math:`\delta_F` that depend on molecular properties. These quantities cannot
simply be summed over constituent elements in a compound, but should instead be
calculated for the material. Currently, we use Bragg's additivity rule to
calculate the collision stopping power as well, but this is not a good
approximation and should be fixed in the future.
where :math:`w_i` is the mass fraction of the :math:`i`-th element and
:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using
:eq:`radiative-stopping-power`. The collision stopping power, however, is a
function of certain quantities such as the mean excitation energy :math:`I` and
the density effect correction :math:`\delta_F` that depend on molecular
properties. These quantities cannot simply be summed over constituent elements
in a compound, but should instead be calculated for the material. The Bethe
formula can be used to find the collision stopping power of the material:
.. math::
:label: material-collision-stopping-power
S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M}
[\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)],
where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass,
:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For
electrons,
.. math::
:label: F-electron
F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2],
while for positrons
.. math::
:label: F-positron
F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 +
4/(\tau + 2)^3].
The density effect correction :math:`\delta_F` takes into account the reduction
of the collision stopping power due to the polarization of the material the
charged particle is passing through by the electric field of the particle.
It can be evaluated using the method described by Sternheimer_, where the
equation for :math:`\delta_F` is
.. math::
:label: density-effect-correction
\delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] -
l^2(1-\beta^2).
Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition,
given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in
the :math:`i`-th subshell. The frequency :math:`l` is the solution of the
equation
.. math::
:label: density-effect-l
\frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2},
where :math:`\bar{v}_i` is defined as
.. math::
:label: density-effect-nubar
\bar{\nu}_i = h\nu_i \rho / h\nu_p.
The plasma energy :math:`h\nu_p` of the medium is given by
.. math::
:label: plasma-frequency
h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}},
where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the
material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator
energy, and :math:`\rho` is an adjustment factor introduced to give agreement
between the experimental values of the oscillator energies and the mean
excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are
defined as
.. math::
:label: density-effect-li
l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\
l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0,
where the second case applies to conduction electrons. For a conductor,
:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective
number of conduction electrons, and :math:`v_n = 0`. The adjustment factor
:math:`\rho` is determined using the equation for the mean excitation energy:
.. math::
:label: mean-excitation-energy
\ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} +
f_n \ln (h\nu_pf_n^{1/2}).
.. _ttb:
@ -891,6 +982,55 @@ direction of the incident charged particle, which is a reasonable approximation
at higher energies when the bremsstrahlung radiation is emitted at small
angles.
-----------------
Photon Production
-----------------
In coupled neutron-photon transport, a source neutron is tracked, and photons
produced from neutron reactions are transported after the neutron's history has
terminated. Since these secondary photons form the photon source for the
problem, it is important to correctly describe their energy and angular
distributions as the accuracy of the calculation relies on the accuracy of this
source. The photon production cross section for a particular reaction :math:`i`
and incident neutron energy :math:`E` is defined as
.. math::
:label: photon-production-xs
\sigma_{\gamma, i}(E) = y_i(E)\sigma_i(E),
where :math:`y_i(E)` is the photon yield corresponding to an incident neutron
reaction having cross section :math:`\sigma_i(E)`.
The yield of photons during neutron transport is determined as the sum of the
photon yields from each individual reaction. In OpenMC, production of photons
is treated in an average sense. That is, the total photon production cross
section is used at a collision site to determine how many photons to produce
rather than the photon production from the reaction that actually took place.
This is partly done for convenience but also because the use of variance
reduction techniques such as implicit capture make it difficult in practice to
directly sample photon production from individual reactions.
In OpenMC, secondary photons are created after a nuclide has been sampled in a
neutron collision. The expected number of photons produced is
.. math::
:label: expected-number-photons
n = w\frac{\sigma_{\gamma}(E)}{\sigma_T(E)},
where :math:`w` is the weight of the neutron, :math:`\sigma_{\gamma}` is the
photon production cross section for the sampled nuclide, and :math:`\sigma_T`
is the total cross section for the nuclide. :math:`\lfloor n \rfloor` photons
are created with an additional photon produced with probability :math:`n -
\lfloor n \rfloor`. Next, a reaction is sampled for each secondary photon. The
probability of sampling the :math:`i`-th reaction is given by
:math:`\sigma_{\gamma, i}(E)/\sum_j\sigma_{\gamma, j}(E)`, where
:math:`\sum_j\sigma_{\gamma, j} = \sigma_{\gamma}` is the total photon
production cross section. The secondary angle and energy distributions
associated with the reaction are used to sample the angle and energy of the
emitted photon.
.. _Koblinger: https://doi.org/10.13182/NSE75-A26663
.. _anomalous scattering: http://pd.chem.ucl.ac.uk/pdnn/diff1/anomscat.htm
@ -906,3 +1046,5 @@ angles.
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: http://www.oecd-nea.org/globalsearch/download.php?doc=77434
.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067

22
docs/source/publications.rst
View file

@ -59,7 +59,7 @@ Coupling and Multi-physics
- Ze-Long Zhao, Yongwei Yang, and Shuang Hong, "`Application of FLUKA and OpenMC
in coupled physics calculation of target and subcritical reactor for ADS
<https://doi.org/10.1007/s41365-018-0539-1>`_," *Nucl. Sci. Tech.*, **30**
<https://doi.org/10.1007/s41365-018-0539-1>`_," *Nucl. Sci. Tech.*, **30**: 10
(2019).
- April Novak, Paul Romano, Brycen Wendt, Ron Rahaman, Elia Merzari, Leslie
@ -134,6 +134,11 @@ Geometry and Visualization
Miscellaneous
-------------
- Faisal Qayyum, Muhammad R. Ali, Awais Zahur, and R. Khan, "`Improvements in
methodology to determine feedback reactivity coefficients
<https://doi.org/10.1007/s41365-019-0588-0>`_," *Nucl. Sci. Tech.*, **30**: 63
(2019).
- Amanda L. Lund and Paul K. Romano, "`Implementation and Validation of Photon
Transport in OpenMC <https://doi.org/10.2172/1490825>`_", Argonne National
Laboratory, Technical Report ANL/MCS-TM-381 (2018).
@ -200,11 +205,21 @@ Miscellaneous
Multigroup Cross Section Generation
-----------------------------------
- William Boyd, Adam Nelson, Paul K. Romano, Samuel Shaner, Benoit Forget, and
Kord Smith, "`Multigroup Cross-Section Generation with the OpenMC Monte Carlo
Particle Transport Code <https://doi.org/10.1080/00295450.2019.1571828>`_,"
*Nucl. Technol.* (2019).
- William Boyd, Benoit Forget, and Kord Smith, "`A single-step framework to
generate spatially self-shielded multi-group cross sections from Monte Carlo
transport simulations <https://doi.org/10.1016/j.anucene.2018.11.017>`_,"
*Ann. Nucl. Energy*, **125**, 261-271 (2019).
- Kun Zhuang, Xiaobin Tang, and Liangzhi Cao, "`Development and verification of
a model for generation of MSFR few-group homogenized cross-sections based on a
Monte Carlo code OpenMC <https://doi.org/10.1016/j.anucene.2018.09.037>`_,"
*Ann. Nucl. Energy*, **124**, 187-197 (2019).
- Changho Lee and Yeon Sang Jung, "Verification of the Cross Section Library
Generated Using OpenMC and MC\ :sup:`2`-3 for PROTEUS," *Proc. PHYSOR*, Cancun,
Mexico, Apr. 22-26 (2018).
@ -454,6 +469,11 @@ Parallelism
Depletion
---------
- Zhao-Qing Liu, Ze-Long Zhao, Yong-Wei Yang, Yu-Cui Gao, Hai-Yan Meng, and
Qing-Yu Gao, "`Development and validation of depletion code system IMPC-Burnup
for ADS <https://doi.org/10.1007/s41365-019-0560-z>`_," *Nucl. Sci. Tech.*,
**30**: 44 (2019).
- Colin Josey, Benoit Forget, and Kord Smith, "`High order methods for the
integration of the Bateman equations and other problems of the form of y' =
F(y,t)y <https://doi.org/10.1016/j.jcp.2017.08.025>`_," *J. Comput. Phys.*,

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

@ -11,6 +11,7 @@ Convenience Functions
:template: myfunction.rst
openmc.model.borated_water
openmc.model.cylinder_from_points
openmc.model.get_hexagonal_prism
openmc.model.get_rectangular_prism
openmc.model.subdivide

10
docs/source/usersguide/cross_sections.rst
View file

@ -192,11 +192,11 @@ Photon Cross Sections
Photon interaction data is needed to run OpenMC with photon transport enabled.
Some of this data, namely bremsstrahlung cross sections from `Seltzer and
Berger`_, stopping powers from the `NIST ESTAR database`_, and Compton profiles
calculated by `Biggs et al.`_ and available in the Geant4 G4EMLOW data file, is
distributed with OpenMC. The rest is available from the NNDC_, which provides
ENDF data from the photo-atomic and atomic relaxation sublibraries of the
ENDF/B-VII.1 library.
Berger`_, mean excitation energy from the `NIST ESTAR database`_, and Compton
profiles calculated by `Biggs et al.`_ and available in the Geant4 G4EMLOW data
file, is distributed with OpenMC. The rest is available from the NNDC_, which
provides ENDF data from the photo-atomic and atomic relaxation sublibraries of
the ENDF/B-VII.1 library.
Most of the pregenerated HDF5 libraries available at https://openmc.mcs.anl.gov
already have photon interaction data included. If you are building a data

12
docs/source/usersguide/geometry.rst
View file

@ -88,7 +88,7 @@ parameters for a sphere are the :math:`x,y,z` coordinates of the center of the
sphere and the radius of the sphere. All of these parameters can be set either
as optional keyword arguments to the class constructor or via attributes::
sphere = openmc.Sphere(R=10.0)
sphere = openmc.Sphere(r=10.0)
# This is equivalent
sphere = openmc.Sphere()
@ -98,7 +98,7 @@ Once a surface has been created, half-spaces can be obtained by applying the
unary ``-`` or ``+`` operators, corresponding to the negative and positive
half-spaces, respectively. For example::
>>> sphere = openmc.Sphere(R=10.0)
>>> sphere = openmc.Sphere(r=10.0)
>>> inside_sphere = -sphere
>>> outside_sphere = +sphere
>>> type(inside_sphere)
@ -140,10 +140,10 @@ may want to specify different behavior for particles passing through a
surface. To specify a vacuum boundary condition, simply change the
:attr:`Surface.boundary_type` attribute to 'vacuum'::
outer_surface = openmc.Sphere(R=100.0, boundary_type='vacuum')
outer_surface = openmc.Sphere(r=100.0, boundary_type='vacuum')
# This is equivalent
outer_surface = openmc.Sphere(R=100.0)
outer_surface = openmc.Sphere(r=100.0)
outer_surface.boundary_type = 'vacuum'
Reflective and periodic boundary conditions can be set with the strings
@ -154,8 +154,8 @@ can be determined automatically. For non-axis-aligned planes, it is necessary to
specify pairs explicitly using the :attr:`Surface.periodic_surface` attribute as
in the following example::
p1 = openmc.Plane(A=0.3, B=5.0, D=1.0, boundary_type='periodic')
p2 = openmc.Plane(A=0.3, B=5.0, D=-1.0, boundary_type='periodic')
p1 = openmc.Plane(a=0.3, b=5.0, d=1.0, boundary_type='periodic')
p2 = openmc.Plane(a=0.3, b=5.0, d=-1.0, boundary_type='periodic')
p1.periodic_surface = p2
Rotationally-periodic boundary conditions can be specified for a pair of

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

@ -243,9 +243,6 @@ coverage
Compile and link code instrumented for coverage analysis. This is typically
used in conjunction with gcov_.
maxcoord
Maximum number of nested coordinate levels in geometry. Defaults to 10.
To set any of these options (e.g. turning on debug mode), the following form
should be used:

18
docs/source/usersguide/tallies.rst
View file

@ -220,6 +220,16 @@ The following tables show all valid scores:
| |multiplicity from (n,2n), (n,3n), and (n,4n) |
| |reactions. |
+----------------------+---------------------------------------------------+
|H1-production |Total production of H1. |
+----------------------+---------------------------------------------------+
|H2-production |Total production of H2 (deuterium). |
+----------------------+---------------------------------------------------+
|H3-production |Total production of H3 (tritium). |
+----------------------+---------------------------------------------------+
|He3-production |Total production of He3. |
+----------------------+---------------------------------------------------+
|He4-production |Total production of He4 (alpha particles). |
+----------------------+---------------------------------------------------+
.. table:: **Miscellaneous scores: units are indicated for each.**
@ -246,6 +256,10 @@ The following tables show all valid scores:
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
+----------------------+---------------------------------------------------+
|heating |Total neutron heating in units of eV per source |
| |particle. This corresponds to MT=301 produced by |
| |NJOY's HEATR module. |
+----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |
| |fission product kinetic energy, prompt and delayed |
@ -281,3 +295,7 @@ The following tables show all valid scores:
|decay-rate |The delayed-nu-fission-weighted decay rate where |
| |the decay rate is in units of inverse seconds. |
+----------------------+---------------------------------------------------+
|damage-energy |Damage energy production in units of eV per source |
| |particle. This corresponds to MT=444 produced by |
| |NJOY's HEATR module. |
+----------------------+---------------------------------------------------+

2
docs/source/usersguide/volume.rst
View file

@ -30,7 +30,7 @@ arguments are not necessary. For example,
::
sphere = openmc.Sphere(R=10.0)
sphere = openmc.Sphere(r=10.0)
cell = openm.Cell(region=-sphere)
vol_calc = openmc.VolumeCalculation([cell], 1000000)

683
examples/jupyter/post-processing.ipynb
View file

File diff suppressed because one or more lines are too long

172
examples/jupyter/search.ipynb
View file

File diff suppressed because one or more lines are too long

6
examples/python/basic/build-xml.py
View file

@ -36,9 +36,9 @@ materials_file.export_to_xml()
###############################################################################
# Instantiate ZCylinder surfaces
surf1 = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=7, name='surf 1')
surf2 = openmc.ZCylinder(surface_id=2, x0=0, y0=0, R=9, name='surf 2')
surf3 = openmc.ZCylinder(surface_id=3, x0=0, y0=0, R=11, name='surf 3')
surf1 = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=7, name='surf 1')
surf2 = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=9, name='surf 2')
surf3 = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=11, name='surf 3')
surf3.boundary_type = 'vacuum'
# Instantiate Cells

2
examples/python/lattice/hexagonal/build-xml.py
View file

@ -43,7 +43,7 @@ left = openmc.XPlane(surface_id=1, x0=-3, name='left')
right = openmc.XPlane(surface_id=2, x0=3, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-4, name='bottom')
top = openmc.YPlane(surface_id=4, y0=4, name='top')
fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, r=0.4)
left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'

6
examples/python/lattice/nested/build-xml.py
View file

@ -39,9 +39,9 @@ left = openmc.XPlane(surface_id=1, x0=-2, name='left')
right = openmc.XPlane(surface_id=2, x0=2, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-2, name='bottom')
top = openmc.YPlane(surface_id=4, y0=2, name='top')
fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, R=0.3)
fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, R=0.2)
fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, r=0.4)
fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, r=0.3)
fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, r=0.2)
left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'

6
examples/python/lattice/simple/build-xml.py
View file

@ -39,9 +39,9 @@ left = openmc.XPlane(surface_id=1, x0=-2, name='left')
right = openmc.XPlane(surface_id=2, x0=2, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-2, name='bottom')
top = openmc.YPlane(surface_id=4, y0=2, name='top')
fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, R=0.3)
fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, R=0.2)
fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, r=0.4)
fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, r=0.3)
fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, r=0.2)
left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'

6
examples/python/pincell/build-xml.py
View file

@ -49,9 +49,9 @@ materials_file.export_to_xml()
###############################################################################
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, R=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, R=0.45720, name='Clad OR')
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=0.45720, name='Clad OR')
left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left')
right = openmc.XPlane(surface_id=5, x0=0.62992, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom')

9
examples/python/pincell_depletion/run_depletion.py
View file

@ -16,7 +16,6 @@ particles = 1000
time_step = 1*24*60*60 # s
final_time = 5*24*60*60 # s
time_steps = np.full(final_time // time_step, time_step)
chain_file = './chain_simple.xml'
power = 174 # W/cm, for 2D simulations only (use W for 3D)
@ -54,9 +53,9 @@ borated_water.add_s_alpha_beta('c_H_in_H2O')
###############################################################################
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, R=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, R=0.45720, name='Clad OR')
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=0.45720, name='Clad OR')
left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left')
right = openmc.XPlane(surface_id=5, x0=0.62992, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom')
@ -100,7 +99,7 @@ geometry = openmc.Geometry(root)
# Compute cell areas
area = {}
area[fuel] = np.pi * fuel_or.coefficients['R'] ** 2
area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2
# Set materials volume for depletion. Set to an area for 2D simulations
uo2.volume = area[fuel]

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

@ -98,7 +98,7 @@ materials_file.export_to_xml()
###############################################################################
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.54, name='Fuel OR')
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.54, name='Fuel OR')
left = openmc.XPlane(surface_id=4, x0=-0.63, name='left')
right = openmc.XPlane(surface_id=5, x0=0.63, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.63, name='bottom')

5
include/openmc/bank.h
View file

@ -20,15 +20,14 @@ namespace openmc {
namespace simulation {
extern "C" int64_t n_bank;
extern std::vector<Particle::Bank> source_bank;
extern std::vector<Particle::Bank> fission_bank;
extern std::vector<Particle::Bank> secondary_bank;
#ifdef _OPENMP
extern std::vector<Particle::Bank> master_fission_bank;
#endif
#pragma omp threadprivate(fission_bank, n_bank)
#pragma omp threadprivate(fission_bank, secondary_bank)
} // namespace simulation

2
include/openmc/capi.h
View file

@ -73,6 +73,8 @@ extern "C" {
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, const char** name);
int openmc_plot_geometry();
int openmc_id_map(const void* slice, int32_t* data_out);
int openmc_property_map(const void* slice, double* data_out);
int openmc_reset();
int openmc_run();
void openmc_set_seed(int64_t new_seed);

34
include/openmc/cell.h
View file

@ -43,6 +43,7 @@ constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
class Cell;
class Universe;
class UniversePartitioner;
namespace model {
extern std::vector<std::unique_ptr<Cell>> cells;
@ -65,6 +66,8 @@ public:
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
std::unique_ptr<UniversePartitioner> partitioner_;
};
//==============================================================================
@ -182,6 +185,7 @@ class DAGCell : public Cell
public:
moab::DagMC* dagmc_ptr_;
DAGCell();
int32_t dag_index_;
bool contains(Position r, Direction u, int32_t on_surface) const;
@ -192,6 +196,36 @@ public:
};
#endif
//==============================================================================
//! Speeds up geometry searches by grouping cells in a search tree.
//
//! Currently this object only works with universes that are divided up by a
//! bunch of z-planes. It could be generalized to other planes, cylinders,
//! and spheres.
//==============================================================================
class UniversePartitioner
{
public:
explicit UniversePartitioner(const Universe& univ);
//! Return the list of cells that could contain the given coordinates.
const std::vector<int32_t>& get_cells(Position r, Direction u) const;
private:
//! A sorted vector of indices to surfaces that partition the universe
std::vector<int32_t> surfs_;
//! Vectors listing the indices of the cells that lie within each partition
//
//! There are n+1 partitions with n surfaces. `partitions_.front()` gives the
//! cells that lie on the negative side of `surfs_.front()`.
//! `partitions_.back()` gives the cells that lie on the positive side of
//! `surfs_.back()`. Otherwise, `partitions_[i]` gives cells sandwiched
//! between `surfs_[i-1]` and `surfs_[i]`.
std::vector<std::vector<int32_t>> partitions_;
};
//==============================================================================
// Non-member functions
//==============================================================================

9
include/openmc/constants.h
View file

@ -225,6 +225,13 @@ constexpr int N_3P {197};
constexpr int N_N3P {198};
constexpr int N_3N2PA {199};
constexpr int N_5N2P {200};
constexpr int N_XP {203};
constexpr int N_XD {204};
constexpr int N_XT {205};
constexpr int N_X3HE {206};
constexpr int N_XA {207};
constexpr int HEATING {301};
constexpr int DAMAGE_ENERGY {444};
constexpr int COHERENT {502};
constexpr int INCOHERENT {504};
constexpr int PAIR_PROD_ELEC {515};
@ -396,8 +403,6 @@ constexpr int LEAKAGE {3};
// Miscellaneous
constexpr int C_NONE {-1};
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
constexpr int ERROR_INT {-2147483647}; // TODO: use <numeric_limits> when F90
// interop is gone
// Interpolation rules
enum class Interpolation {

10
include/openmc/error.h
View file

@ -7,6 +7,12 @@
#include "openmc/capi.h"
#ifdef __GNUC__
#define UNREACHABLE() __builtin_unreachable()
#else
#define UNREACHABLE() (void)0
#endif
namespace openmc {
inline void
@ -29,13 +35,13 @@ set_errmsg(const std::stringstream& message)
[[noreturn]] void fatal_error(const std::string& message, int err=-1);
inline
[[noreturn]] inline
void fatal_error(const std::stringstream& message)
{
fatal_error(message.str());
}
inline
[[noreturn]] inline
void fatal_error(const char* message)
{
fatal_error({message, std::strlen(message)});

37
include/openmc/geometry.h
View file

@ -1,6 +1,8 @@
#ifndef OPENMC_GEOMETRY_H
#define OPENMC_GEOMETRY_H
#include <array>
#include <cmath>
#include <cstdint>
#include <vector>
@ -15,18 +17,37 @@ namespace openmc {
namespace model {
extern "C" int root_universe;
extern int root_universe; //!< Index of root universe
extern int n_coord_levels; //!< Number of CSG coordinate levels
extern std::vector<int64_t> overlap_check_count;
} // namespace model
//==============================================================================
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface_index {0}; //!< if boundary is surface, index in surfaces vector
int coord_level; //!< coordinate level after crossing boundary
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
//==============================================================================
//! Check two distances by coincidence tolerance
//==============================================================================
inline bool coincident(double d1, double d2) {
return std::abs(d1 - d2) < FP_COINCIDENT;
}
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p);
bool check_cell_overlap(Particle* p);
//==============================================================================
//! Locate a particle in the geometry tree and set its geometry data fields.
@ -40,23 +61,19 @@ check_cell_overlap(Particle* p);
//! valid geometry coordinate stack.
//==============================================================================
extern "C" bool
find_cell(Particle* p, bool use_neighbor_lists);
bool find_cell(Particle* p, bool use_neighbor_lists);
//==============================================================================
//! Move a particle into a new lattice tile.
//==============================================================================
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3]);
void cross_lattice(Particle* p, const BoundaryInfo& boundary);
//==============================================================================
//! Find the next boundary a particle will intersect.
//==============================================================================
extern "C" void
distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
int lattice_translation[3], int* next_level);
BoundaryInfo distance_to_boundary(Particle* p);
} // namespace openmc

9
include/openmc/geometry_aux.h
View file

@ -57,7 +57,7 @@ void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
//! \return The index of the root universe.
//==============================================================================
extern "C" int32_t find_root_universe();
int32_t find_root_universe();
//==============================================================================
//! Populate all data structures needed for distribcells.
@ -74,7 +74,7 @@ void prepare_distribcell();
//! the root universe).
//==============================================================================
extern "C" void count_cell_instances(int32_t univ_indx);
void count_cell_instances(int32_t univ_indx);
//==============================================================================
//! Recursively search through universes and count universe instances.
@ -84,8 +84,7 @@ extern "C" void count_cell_instances(int32_t univ_indx);
//! search_univ.
//==============================================================================
extern "C" int
count_universe_instances(int32_t search_univ, int32_t target_univ_id);
int count_universe_instances(int32_t search_univ, int32_t target_univ_id);
//==============================================================================
//! Build a character array representing the path to a distribcell instance.
@ -107,7 +106,7 @@ distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset);
//! \return The number of coordinate levels.
//==============================================================================
extern "C" int maximum_levels(int32_t univ);
int maximum_levels(int32_t univ);
//==============================================================================
//! Deallocates global vectors and maps for cells, universes, and lattices.

3
include/openmc/lattice.h
View file

@ -152,8 +152,7 @@ public:
int indx_; //!< An index to a Lattice universes or offsets array.
LatticeIter(Lattice &lat, int indx)
: lat_(lat),
indx_(indx)
: indx_(indx), lat_(lat)
{}
bool operator==(const LatticeIter &rhs) {return (indx_ == rhs.indx_);}

19
include/openmc/material.h
View file

@ -47,7 +47,7 @@ public:
explicit Material(pugi::xml_node material_node);
// Methods
void calculate_xs(const Particle& p) const;
void calculate_xs(Particle& p) const;
//! Assign thermal scattering tables to specific nuclides within the material
//! so the code knows when to apply bound thermal scattering data
@ -95,20 +95,33 @@ public:
std::unique_ptr<Bremsstrahlung> ttb_;
private:
//! Calculate the collision stopping power
void collision_stopping_power(double* s_col, bool positron);
//! Initialize bremsstrahlung data
void init_bremsstrahlung();
//! Normalize density
void normalize_density();
void calculate_neutron_xs(const Particle& p) const;
void calculate_photon_xs(const Particle& p) const;
void calculate_neutron_xs(Particle& p) const;
void calculate_photon_xs(Particle& p) const;
};
//==============================================================================
// Non-member functions
//==============================================================================
//! Calculate Sternheimer adjustment factor
double sternheimer_adjustment(const std::vector<double>& f, const
std::vector<double>& e_b_sq, double e_p_sq, double n_conduction, double
log_I, double tol, int max_iter);
//! Calculate density effect correction
double density_effect(const std::vector<double>& f, const std::vector<double>&
e_b_sq, double e_p_sq, double n_conduction, double rho, double E, double tol,
int max_iter);
//! Read material data from materials.xml
void read_materials_xml();

42
include/openmc/mesh.h
View file

@ -43,7 +43,7 @@ public:
// Methods
//! Determine which bins were crossed by a particle
//!
//
//! \param[in] p Particle to check
//! \param[out] bins Bins that were crossed
//! \param[out] lengths Fraction of tracklength in each bin
@ -51,58 +51,56 @@ public:
std::vector<double>& lengths) const;
//! Determine which surface bins were crossed by a particle
//!
//
//! \param[in] p Particle to check
//! \param[out] bins Surface bins that were crossed
void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
//! Get bin at a given position in space
//!
//
//! \param[in] r Position to get bin for
//! \return Mesh bin
int get_bin(Position r) const;
//! Get bin given mesh indices
//!
//
//! \param[in] Array of mesh indices
//! \return Mesh bin
int get_bin_from_indices(const int* ijk) const;
//! Get mesh indices given a position
//!
//
//! \param[in] r Position to get indices for
//! \param[out] ijk Array of mesh indices
//! \param[out] in_mesh Whether position is in mesh
void get_indices(Position r, int* ijk, bool* in_mesh) const;
//! Get mesh indices corresponding to a mesh bin
//!
//
//! \param[in] bin Mesh bin
//! \param[out] ijk Mesh indices
void get_indices_from_bin(int bin, int* ijk) const;
//! Check if a line connected by two points intersects the mesh
//!
//! \param[in] r0 Starting position
//! Check where a line segment intersects the mesh and if it intersects at all
//
//! \param[in,out] r0 In: starting position, out: intersection point
//! \param[in] r1 Ending position
//! \return Whether line connecting r0 and r1 intersects mesh
bool intersects(Position r0, Position r1) const;
//! \param[out] ijk Indices of the mesh bin containing the intersection point
//! \return Whether the line segment connecting r0 and r1 intersects mesh
bool intersects(Position& r0, Position r1, int* ijk) const;
//! Write mesh data to an HDF5 group
//!
//
//! \param[in] group HDF5 group
void to_hdf5(hid_t group) const;
//! Count number of bank sites in each mesh bin / energy bin
//!
//! \param[in] n Number of bank sites
//
//! \param[in] bank Array of bank sites
//! \param[in] n_energy Number of energies
//! \param[in] energies Array of energies
//! \param[out] Whether any bank sites are outside the mesh
//! \return Array indicating number of sites in each mesh/energy bin
xt::xarray<double> count_sites(int64_t n, const Particle::Bank* bank,
int n_energy, const double* energies, bool* outside) const;
xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
bool* outside) const;
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
@ -113,9 +111,9 @@ public:
xt::xarray<double> width_; //!< Width of each mesh element
private:
bool intersects_1d(Position r0, Position r1) const;
bool intersects_2d(Position r0, Position r1) const;
bool intersects_3d(Position r0, Position r1) const;
bool intersects_1d(Position& r0, Position r1, int* ijk) const;
bool intersects_2d(Position& r0, Position r1, int* ijk) const;
bool intersects_3d(Position& r0, Position r1, int* ijk) const;
};
//==============================================================================
@ -123,10 +121,12 @@ private:
//==============================================================================
//! Read meshes from either settings/tallies
//
//! \param[in] root XML node
void read_meshes(pugi::xml_node root);
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
void meshes_to_hdf5(hid_t group);

82
include/openmc/nuclide.h
View file

@ -13,6 +13,7 @@
#include "openmc/constants.h"
#include "openmc/endf.h"
#include "openmc/particle.h"
#include "openmc/reaction.h"
#include "openmc/reaction_product.h"
#include "openmc/urr.h"
@ -20,69 +21,6 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
constexpr double CACHE_INVALID {-1.0};
//==============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy
//==============================================================================
struct NuclideMicroXS {
// Microscopic cross sections in barns
double total; //!< total cross section
double absorption; //!< absorption (disappearance)
double fission; //!< fission
double nu_fission; //!< neutron production from fission
double elastic; //!< If sab_frac is not 1 or 0, then this value is
//!< averaged over bound and non-bound nuclei
double thermal; //!< Bound thermal elastic & inelastic scattering
double thermal_elastic; //!< Bound thermal elastic scattering
double photon_prod; //!< microscopic photon production xs
// Cross sections for depletion reactions (note that these are not stored in
// macroscopic cache)
double reaction[DEPLETION_RX.size()];
// Indicies and factors needed to compute cross sections from the data tables
int index_grid; //!< Index on nuclide energy grid
int index_temp; //!< Temperature index for nuclide
double interp_factor; //!< Interpolation factor on nuc. energy grid
int index_sab {-1}; //!< Index in sab_tables
int index_temp_sab; //!< Temperature index for sab_tables
double sab_frac; //!< Fraction of atoms affected by S(a,b)
bool use_ptable; //!< In URR range with probability tables?
// Energy and temperature last used to evaluate these cross sections. If
// these values have changed, then the cross sections must be re-evaluated.
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
};
//==============================================================================
// MATERIALMACROXS contains cached macroscopic cross sections for the material a
// particle is traveling through
//==============================================================================
struct MaterialMacroXS {
double total; //!< macroscopic total xs
double absorption; //!< macroscopic absorption xs
double fission; //!< macroscopic fission xs
double nu_fission; //!< macroscopic production xs
double photon_prod; //!< macroscopic photon production xs
// Photon cross sections
double coherent; //!< macroscopic coherent xs
double incoherent; //!< macroscopic incoherent xs
double photoelectric; //!< macroscopic photoelectric xs
double pair_production; //!< macroscopic pair production xs
};
//==============================================================================
// Data for a nuclide
//==============================================================================
@ -102,14 +40,13 @@ public:
//! Initialize logarithmic grid for energy searches
void init_grid();
void calculate_xs(int i_sab, double E, int i_log_union,
double sqrtkT, double sab_frac);
void calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p);
void calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_frac);
void calculate_sab_xs(int i_sab, double sab_frac, Particle& p);
// Methods
double nu(double E, EmissionMode mode, int group=0) const;
void calculate_elastic_xs() const;
void calculate_elastic_xs(Particle& p) const;
//! Determines the microscopic 0K elastic cross section at a trial relative
//! energy used in resonance scattering
@ -117,7 +54,7 @@ public:
//! \brief Determines cross sections in the unresolved resonance range
//! from probability tables.
void calculate_urr_xs(int i_temp, double E) const;
void calculate_urr_xs(int i_temp, Particle& p) const;
// Data members
std::string name_; //!< Name of nuclide, e.g. "U235"
@ -194,15 +131,6 @@ extern std::unordered_map<std::string, int> nuclide_map;
} // namespace data
namespace simulation {
// Cross section caches
extern NuclideMicroXS* micro_xs;
extern MaterialMacroXS material_xs;
#pragma omp threadprivate(micro_xs, material_xs)
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================

214
include/openmc/particle.h
View file

@ -6,9 +6,11 @@
#include <array>
#include <cstdint>
#include <memory> // for unique_ptr
#include <sstream>
#include <string>
#include "openmc/constants.h"
#include "openmc/position.h"
namespace openmc {
@ -24,15 +26,14 @@ namespace openmc {
// use to store the bins for delayed group tallies.
constexpr int MAX_DELAYED_GROUPS {8};
// Maximum number of secondary particles created
constexpr int MAX_SECONDARY {1000};
// Maximum number of lost particles
constexpr int MAX_LOST_PARTICLES {10};
// Maximum number of lost particles, relative to the total number of particles
constexpr double REL_MAX_LOST_PARTICLES {1.0e-6};
constexpr double CACHE_INVALID {-1.0};
//==============================================================================
// Class declarations
//==============================================================================
@ -52,12 +53,88 @@ struct LocalCoord {
void reset();
};
//==============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy
//==============================================================================
struct NuclideMicroXS {
// Microscopic cross sections in barns
double total; //!< total cross section
double absorption; //!< absorption (disappearance)
double fission; //!< fission
double nu_fission; //!< neutron production from fission
double elastic; //!< If sab_frac is not 1 or 0, then this value is
//!< averaged over bound and non-bound nuclei
double thermal; //!< Bound thermal elastic & inelastic scattering
double thermal_elastic; //!< Bound thermal elastic scattering
double photon_prod; //!< microscopic photon production xs
// Cross sections for depletion reactions (note that these are not stored in
// macroscopic cache)
double reaction[DEPLETION_RX.size()];
// Indicies and factors needed to compute cross sections from the data tables
int index_grid; //!< Index on nuclide energy grid
int index_temp; //!< Temperature index for nuclide
double interp_factor; //!< Interpolation factor on nuc. energy grid
int index_sab {-1}; //!< Index in sab_tables
int index_temp_sab; //!< Temperature index for sab_tables
double sab_frac; //!< Fraction of atoms affected by S(a,b)
bool use_ptable; //!< In URR range with probability tables?
// Energy and temperature last used to evaluate these cross sections. If
// these values have changed, then the cross sections must be re-evaluated.
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
};
//==============================================================================
//! Cached microscopic photon cross sections for a particular element at the
//! current energy
//==============================================================================
struct ElementMicroXS {
int index_grid; //!< index on element energy grid
double last_E {0.0}; //!< last evaluated energy in [eV]
double interp_factor; //!< interpolation factor on energy grid
double total; //!< microscopic total photon xs
double coherent; //!< microscopic coherent xs
double incoherent; //!< microscopic incoherent xs
double photoelectric; //!< microscopic photoelectric xs
double pair_production; //!< microscopic pair production xs
};
//==============================================================================
// MACROXS contains cached macroscopic cross sections for the material a
// particle is traveling through
//==============================================================================
struct MacroXS {
double total; //!< macroscopic total xs
double absorption; //!< macroscopic absorption xs
double fission; //!< macroscopic fission xs
double nu_fission; //!< macroscopic production xs
double photon_prod; //!< macroscopic photon production xs
// Photon cross sections
double coherent; //!< macroscopic coherent xs
double incoherent; //!< macroscopic incoherent xs
double photoelectric; //!< macroscopic photoelectric xs
double pair_production; //!< macroscopic pair production xs
};
//============================================================================
//! State of a particle being transported through geometry
//============================================================================
class Particle {
public:
//==========================================================================
// Aliases and type definitions
//! Particle types
enum class Type {
neutron, photon, electron, positron
@ -73,19 +150,87 @@ public:
Type particle;
};
//==========================================================================
// Constructors
Particle();
//==========================================================================
// Methods and accessors
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
//! resets all coordinate levels for the particle
void clear();
//! create a secondary particle
//
//! stores the current phase space attributes of the particle in the
//! secondary bank and increments the number of sites in the secondary bank.
//! \param u Direction of the secondary particle
//! \param E Energy of the secondary particle in [eV]
//! \param type Particle type
void create_secondary(Direction u, double E, Type type) const;
//! initialize from a source site
//
//! initializes a particle from data stored in a source site. The source
//! site may have been produced from an external source, from fission, or
//! simply as a secondary particle.
//! \param src Source site data
void from_source(const Bank* src);
//! Transport a particle from birth to death
void transport();
//! Cross a surface and handle boundary conditions
void cross_surface();
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
void mark_as_lost(const char* message);
void mark_as_lost(const std::string& message)
{mark_as_lost(message.c_str());}
void mark_as_lost(const std::stringstream& message)
{mark_as_lost(message.str());}
//! create a particle restart HDF5 file
void write_restart() const;
//==========================================================================
// Data members
// Cross section caches
std::vector<NuclideMicroXS> neutron_xs_; //!< Microscopic neutron cross sections
std::vector<ElementMicroXS> photon_xs_; //!< Microscopic photon cross sections
MacroXS macro_xs_; //!< Macroscopic cross sections
int64_t id_; //!< Unique ID
Type type_ {Type::neutron}; //!< Particle type (n, p, e, etc.)
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels
std::vector<LocalCoord> coord_; //!< coordinates for all levels
// Particle coordinates before crossing a surface
int n_coord_last_ {1}; //!< number of current coordinates
int cell_last_[MAX_COORD]; //!< coordinates for all levels
std::vector<int> cell_last_; //!< coordinates for all levels
// Energy data
double E_; //!< post-collision energy in eV
@ -135,65 +280,6 @@ public:
// Track output
bool write_track_ {false};
// Secondary particles created
int64_t n_secondary_ {};
Bank secondary_bank_[MAX_SECONDARY];
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
//! resets all coordinate levels for the particle
void clear();
//! create a secondary particle
//
//! stores the current phase space attributes of the particle in the
//! secondary bank and increments the number of sites in the secondary bank.
//! \param u Direction of the secondary particle
//! \param E Energy of the secondary particle in [eV]
//! \param type Particle type
void create_secondary(Direction u, double E, Type type);
//! initialize from a source site
//
//! initializes a particle from data stored in a source site. The source
//! site may have been produced from an external source, from fission, or
//! simply as a secondary particle.
//! \param src Source site data
void from_source(const Bank* src);
//! Transport a particle from birth to death
void transport();
//! Cross a surface and handle boundary conditions
void cross_surface();
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
void mark_as_lost(const char* message);
void mark_as_lost(const std::string& message)
{mark_as_lost(message.c_str());}
void mark_as_lost(const std::stringstream& message)
{mark_as_lost(message.str());}
//! create a particle restart HDF5 file
void write_restart() const;
};
} // namespace openmc

26
include/openmc/photon.h
View file

@ -42,7 +42,7 @@ public:
PhotonInteraction(hid_t group, int i_element);
// Methods
void calculate_xs(double E) const;
void calculate_xs(Particle& p) const;
void compton_scatter(double alpha, bool doppler, double* alpha_out,
double* mu, int* i_shell) const;
@ -87,7 +87,8 @@ public:
// Stopping power data
double I_; // mean excitation energy
xt::xtensor<double, 1> stopping_power_collision_;
xt::xtensor<int, 1> n_electrons_;
xt::xtensor<double, 1> ionization_energy_;
xt::xtensor<double, 1> stopping_power_radiative_;
// Bremsstrahlung scaled DCS
@ -97,22 +98,6 @@ private:
void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const;
};
//==============================================================================
//! Cached microscopic photon cross sections for a particular element at the
//! current energy
//==============================================================================
struct ElementMicroXS {
int index_grid; //!< index on element energy grid
double last_E {0.0}; //!< last evaluated energy in [eV]
double interp_factor; //!< interpolation factor on energy grid
double total; //!< microscopic total photon xs
double coherent; //!< microscopic coherent xs
double incoherent; //!< microscopic incoherent xs
double photoelectric; //!< microscopic photoelectric xs
double pair_production; //!< microscopic pair production xs
};
//==============================================================================
// Non-member functions
//==============================================================================
@ -135,11 +120,6 @@ extern std::unordered_map<std::string, int> element_map;
} // namespace data
namespace simulation {
extern ElementMicroXS* micro_photon_xs;
#pragma omp threadprivate(micro_photon_xs)
} // namespace simulation
} // namespace openmc
#endif // OPENMC_PHOTON_H

15
include/openmc/physics.h
View file

@ -7,6 +7,8 @@
#include "openmc/position.h"
#include "openmc/reaction.h"
#include <vector>
namespace openmc {
//==============================================================================
@ -47,24 +49,23 @@ int sample_nuclide(const Particle* p);
//! Determine the average total, prompt, and delayed neutrons produced from
//! fission and creates appropriate bank sites.
void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Particle::Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
std::vector<Particle::Bank>& bank);
int sample_element(Particle* p);
Reaction* sample_fission(int i_nuclide, double E);
Reaction* sample_fission(int i_nuclide, const Particle* p);
void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product);
void absorption(Particle* p, int i_nuclide);
void scatter(Particle*, int i_nuclide);
//! Treats the elastic scattering of a neutron with a target.
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E,
Direction& u, double& mu_lab);
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
Particle* p);
void sab_scatter(int i_nuclide, int i_sab, double& E,
Direction& u, double& mu);
void sab_scatter(int i_nuclide, int i_sab, Particle* p);
//! samples the target velocity. The constant cross section free gas model is
//! the default method. Methods for correctly accounting for the energy

9
include/openmc/physics_mg.h
View file

@ -8,6 +8,8 @@
#include "openmc/particle.h"
#include "openmc/nuclide.h"
#include <vector>
namespace openmc {
//! \brief samples particle behavior after a collision event.
@ -31,12 +33,9 @@ scatter(Particle* p);
//! \brief Determines the average total, prompt and delayed neutrons produced
//! from fission and creates the appropriate bank sites.
//! \param p Particle to operate on
//! \param bank_array The particle bank to populate
//! \param size_bank Number of particles currently in the bank
//! \param bank_array_size Allocated size of the bank
//! \param bank The particle bank to populate
void
create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size);
create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank);
//! \brief Handles an absorption event
//! \param p Particle to operate on

128
include/openmc/plot.h
View file

@ -10,6 +10,8 @@
#include "hdf5.h"
#include "openmc/position.h"
#include "openmc/constants.h"
#include "openmc/cell.h"
#include "openmc/geometry.h"
#include "openmc/particle.h"
#include "openmc/xml_interface.h"
@ -53,6 +55,28 @@ struct RGBColor {
typedef xt::xtensor<RGBColor, 2> ImageData;
struct IdData {
// Constructor
IdData(size_t h_res, size_t v_res);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level);
// Members
xt::xtensor<int32_t, 3> data_; //!< 2D array of cell & material ids
};
struct PropertyData {
// Constructor
PropertyData(size_t h_res, size_t v_res);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level);
// Members
xt::xtensor<double, 3> data_; //!< 2D array of temperature & density data
};
enum class PlotType {
slice = 1,
voxel = 2
@ -65,16 +89,98 @@ enum class PlotBasis {
};
enum class PlotColorBy {
cells = 1,
mats = 2
cells = 0,
mats = 1
};
//===============================================================================
// Plot class
//===============================================================================
class PlotBase {
public:
template<class T> T get_map() const;
class Plot
{
// Members
public:
Position origin_; //!< Plot origin in geometry
Position width_; //!< Plot width in geometry
PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
std::array<size_t, 3> pixels_; //!< Plot size in pixels
int level_; //!< Plot universe level
};
template<class T>
T PlotBase::get_map() const {
size_t width = pixels_[0];
size_t height = pixels_[1];
// get pixel size
double in_pixel = (width_[0])/static_cast<double>(width);
double out_pixel = (width_[1])/static_cast<double>(height);
// size data array
T data(width, height);
// setup basis indices and initial position centered on pixel
int in_i, out_i;
Position xyz = origin_;
switch(basis_) {
case PlotBasis::xy :
in_i = 0;
out_i = 1;
break;
case PlotBasis::xz :
in_i = 0;
out_i = 2;
break;
case PlotBasis::yz :
in_i = 1;
out_i = 2;
break;
#ifdef __GNUC__
default:
__builtin_unreachable();
#endif
}
// set initial position
xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.;
xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.;
// arbitrary direction
Direction dir = {0.7071, 0.7071, 0.0};
#pragma omp parallel
{
Particle p;
p.r() = xyz;
p.u() = dir;
p.coord_[0].universe = model::root_universe;
int level = level_;
int j{};
#pragma omp for
for (int y = 0; y < height; y++) {
p.r()[out_i] = xyz[out_i] - out_pixel * y;
for (int x = 0; x < width; x++) {
p.r()[in_i] = xyz[in_i] + in_pixel * x;
p.n_coord_ = 1;
// local variables
bool found_cell = find_cell(&p, 0);
j = p.n_coord_ - 1;
if (level >=0) {j = level + 1;}
if (found_cell) {
data.set_value(y, x, p, j);
}
} // inner for
} // outer for
} // omp parallel
return data;
}
class Plot : public PlotBase {
public:
// Constructor
@ -95,17 +201,12 @@ private:
void set_meshlines(pugi::xml_node plot_node);
void set_mask(pugi::xml_node plot_node);
// Members
// Members
public:
int id_; //!< Plot ID
PlotType type_; //!< Plot type (Slice/Voxel)
PlotColorBy color_by_; //!< Plot coloring (cell/material)
Position origin_; //!< Plot origin in geometry
Position width_; //!< Plot width in geometry
PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
std::array<int, 3> pixels_; //!< Plot size in pixels
int meshlines_width_; //!< Width of lines added to the plot
int level_; //!< Plot universe level
int index_meshlines_mesh_; //!< Index of the mesh to draw on the plot
RGBColor meshlines_color_; //!< Color of meshlines on the plot
RGBColor not_found_; //!< Plot background color
@ -127,13 +228,6 @@ void draw_mesh_lines(Plot pl, ImageData& data);
//! \param[out] image data associated with the plot object
void output_ppm(Plot pl, const ImageData& data);
//! Get the rgb color for a given particle position in a plot
//! \param[in] particle with position for current pixel
//! \param[in] plot object
//! \param[out] rgb color
//! \param[out] cell or material id for particle position
void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id);
//! Initialize a voxel file
//! \param[in] id of an open hdf5 file
//! \param[in] dimensions of the voxel file (dx, dy, dz)

8
include/openmc/simulation.h
View file

@ -37,19 +37,15 @@ extern "C" int restart_batch; //!< batch at which a restart job resumed
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
extern "C" int total_gen; //!< total number of generations simulated
extern double total_weight; //!< Total source weight in a batch
extern "C" int64_t work; //!< number of particles per process
extern int64_t work_per_rank; //!< number of particles per MPI rank
extern std::vector<double> k_generation;
extern std::vector<int64_t> work_index;
// Threadprivate variables
extern "C" bool trace; //!< flag to show debug information
#ifdef _OPENMP
extern "C" int n_threads; //!< number of OpenMP threads
extern "C" int thread_id; //!< ID of a given thread
#endif
#pragma omp threadprivate(current_work, thread_id, trace)
#pragma omp threadprivate(current_work, trace)
} // namespace simulation

40
include/openmc/surface.h
View file

@ -131,6 +131,8 @@ class DAGSurface : public Surface
public:
moab::DagMC* dagmc_ptr_;
DAGSurface();
int32_t dag_index_;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
@ -178,7 +180,6 @@ public:
class SurfaceXPlane : public PeriodicSurface
{
double x0_;
public:
explicit SurfaceXPlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
@ -188,6 +189,8 @@ public:
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
double x0_;
};
//==============================================================================
@ -198,7 +201,6 @@ public:
class SurfaceYPlane : public PeriodicSurface
{
double y0_;
public:
explicit SurfaceYPlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
@ -208,6 +210,8 @@ public:
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
double y0_;
};
//==============================================================================
@ -218,7 +222,6 @@ public:
class SurfaceZPlane : public PeriodicSurface
{
double z0_;
public:
explicit SurfaceZPlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
@ -228,6 +231,8 @@ public:
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
double z0_;
};
//==============================================================================
@ -238,7 +243,6 @@ public:
class SurfacePlane : public PeriodicSurface
{
double A_, B_, C_, D_;
public:
explicit SurfacePlane(pugi::xml_node surf_node);
double evaluate(Position r) const;
@ -248,6 +252,8 @@ public:
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
double A_, B_, C_, D_;
};
//==============================================================================
@ -259,13 +265,14 @@ public:
class SurfaceXCylinder : public CSGSurface
{
double y0_, z0_, radius_;
public:
explicit SurfaceXCylinder(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double y0_, z0_, radius_;
};
//==============================================================================
@ -277,13 +284,14 @@ public:
class SurfaceYCylinder : public CSGSurface
{
double x0_, z0_, radius_;
public:
explicit SurfaceYCylinder(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double x0_, z0_, radius_;
};
//==============================================================================
@ -295,13 +303,14 @@ public:
class SurfaceZCylinder : public CSGSurface
{
double x0_, y0_, radius_;
public:
explicit SurfaceZCylinder(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double x0_, y0_, radius_;
};
//==============================================================================
@ -313,13 +322,14 @@ public:
class SurfaceSphere : public CSGSurface
{
double x0_, y0_, z0_, radius_;
public:
explicit SurfaceSphere(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double x0_, y0_, z0_, radius_;
};
//==============================================================================
@ -331,13 +341,14 @@ public:
class SurfaceXCone : public CSGSurface
{
double x0_, y0_, z0_, radius_sq_;
public:
explicit SurfaceXCone(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double x0_, y0_, z0_, radius_sq_;
};
//==============================================================================
@ -349,13 +360,14 @@ public:
class SurfaceYCone : public CSGSurface
{
double x0_, y0_, z0_, radius_sq_;
public:
explicit SurfaceYCone(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double x0_, y0_, z0_, radius_sq_;
};
//==============================================================================
@ -367,13 +379,14 @@ public:
class SurfaceZCone : public CSGSurface
{
double x0_, y0_, z0_, radius_sq_;
public:
explicit SurfaceZCone(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
double x0_, y0_, z0_, radius_sq_;
};
//==============================================================================
@ -384,14 +397,15 @@ public:
class SurfaceQuadric : public CSGSurface
{
// Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
public:
explicit SurfaceQuadric(pugi::xml_node surf_node);
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
// Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
};
//==============================================================================

6
include/openmc/tallies/tally_scoring.h
View file

@ -59,14 +59,14 @@ private:
//! since collisions do not occur in voids.
//
//! \param p The particle being tracked
void score_collision_tally(const Particle* p);
void score_collision_tally(Particle* p);
//! Score tallies based on a simple count of events (for continuous energy).
//
//! Analog tallies are triggered at every collision, not every event.
//
//! \param p The particle being tracked
void score_analog_tally_ce(const Particle* p);
void score_analog_tally_ce(Particle* p);
//! Score tallies based on a simple count of events (for multigroup).
//
@ -83,7 +83,7 @@ void score_analog_tally_mg(const Particle* p);
//
//! \param p The particle being tracked
//! \param distance The distance in [cm] traveled by the particle
void score_tracklength_tally(const Particle* p, double distance);
void score_tracklength_tally(Particle* p, double distance);
//! Score surface or mesh-surface tallies for particle currents.
//

2
include/openmc/thermal.h
View file

@ -10,7 +10,7 @@
#include "xtensor/xtensor.hpp"
#include "openmc/hdf5_interface.h"
#include "openmc/nuclide.h"
#include "openmc/particle.h"
namespace openmc {

2
openmc/capi/__init__.py
View file

@ -42,7 +42,6 @@ else:
def _dagmc_enabled():
return c_bool.in_dll(_dll, "dagmc_enabled").value
from .error import *
from .core import *
from .nuclide import *
@ -53,3 +52,4 @@ from .filter import *
from .tally import *
from .settings import settings
from .math import *
from .plot import *

3
openmc/capi/core.py
View file

@ -17,7 +17,8 @@ class _Bank(Structure):
('u', c_double*3),
('E', c_double),
('wgt', c_double),
('delayed_group', c_int)]
('delayed_group', c_int),
('particle', c_int)]
# Define input type for numpy arrays that will be passed into C++ functions

245
openmc/capi/plot.py Normal file
View file

@ -0,0 +1,245 @@
from ctypes import c_int, c_size_t, c_int32, c_double, Structure, POINTER
from . import _dll
from .error import _error_handler
import numpy as np
class _Position(Structure):
"""Definition of an xyz location in space with underlying c-types
C-type Attributes
-----------------
x : c_double
Position's x value (default: 0.0)
y : c_double
Position's y value (default: 0.0)
z : c_double
Position's z value (default: 0.0)
"""
_fields_ = [('x', c_double),
('y', c_double),
('z', c_double)]
def __getitem__(self, idx):
if idx == 0:
return self.x
elif idx == 1:
return self.y
elif idx == 2:
return self.z
else:
raise IndexError("{} index is invalid for _Position".format(idx))
def __setitem__(self, idx, val):
if idx == 0:
self.x = val
elif idx == 1:
self.y = val
elif idx == 2:
self.z = val
else:
raise IndexError("{} index is invalid for _Position".format(idx))
def __repr__(self):
return "({}, {}, {})".format(self.x, self.y, self.z)
class _PlotBase(Structure):
"""A structure defining a 2-D geometry slice with underlying c-types
C-Type Attributes
-----------------
origin : openmc.capi.plot._Position
A position defining the origin of the plot.
width_ : openmc.capi.plot._Position
The width of the plot along the x, y, and z axes, respectively
basis_ : c_int
The axes basis of the plot view.
pixels_ : c_size_t[3]
The resolution of the plot in the horizontal and vertical dimensions
level_ : c_int
The universe level for the plot view
Attributes
----------
origin : tuple or list of ndarray
Origin (center) of the plot
width : float
The horizontal dimension of the plot in geometry units (cm)
height : float
The vertical dimension of the plot in geometry units (cm)
basis : string
One of {'xy', 'xz', 'yz'} indicating the horizontal and vertical
axes of the plot.
h_res : int
The horizontal resolution of the plot in pixels
v_res : int
The vertical resolution of the plot in pixels
level : int
The universe level for the plot (default: -1 -> all universes shown)
"""
_fields_ = [('origin_', _Position),
('width_', _Position),
('basis_', c_int),
('pixels_', 3*c_size_t),
('level_', c_int)]
def __init__(self):
self.level_ = -1
@property
def origin(self):
return self.origin_
@property
def width(self):
return self.width_.x
@property
def height(self):
return self.width_.y
@property
def basis(self):
if self.basis_ == 1:
return 'xy'
elif self.basis_ == 2:
return 'xz'
elif self.basis_ == 3:
return 'yz'
raise ValueError("Plot basis {} is invalid".format(self.basis_))
@property
def h_res(self):
return self.pixels_[0]
@property
def v_res(self):
return self.pixels_[1]
@property
def level(self):
return int(self.level_)
@origin.setter
def origin(self, origin):
self.origin_.x = origin[0]
self.origin_.y = origin[1]
self.origin_.z = origin[2]
@width.setter
def width(self, width):
self.width_.x = width
@height.setter
def height(self, height):
self.width_.y = height
@basis.setter
def basis(self, basis):
if isinstance(basis, str):
valid_bases = ('xy', 'xz', 'yz')
basis = basis.lower()
if basis not in valid_bases:
raise ValueError("{} is not a valid plot basis.".format(basis))
if basis == 'xy':
self.basis_ = 1
elif basis == 'xz':
self.basis_ = 2
elif basis == 'yz':
self.basis_ = 3
return
if isinstance(basis, int):
valid_bases = (1, 2, 3)
if basis not in valid_bases:
raise ValueError("{} is not a valid plot basis.".format(basis))
self.basis_ = basis
return
raise ValueError("{} of type {} is an"
" invalid plot basis".format(basis, type(basis)))
@h_res.setter
def h_res(self, h_res):
self.pixels_[0] = h_res
@v_res.setter
def v_res(self, v_res):
self.pixels_[1] = v_res
@level.setter
def level(self, level):
self.level_ = level
def __repr__(self):
out_str = ["-----",
"Plot:",
"-----",
"Origin: {}".format(self.origin),
"Width: {}".format(self.width),
"Height: {}".format(self.height),
"Basis: {}".format(self.basis),
"HRes: {}".format(self.h_res),
"VRes: {}".format(self.v_res),
"Level: {}".format(self.level)]
return '\n'.join(out_str)
_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), POINTER(c_int32)]
_dll.openmc_id_map.restype = c_int
_dll.openmc_id_map.errcheck = _error_handler
def id_map(plot):
"""
Generate a 2-D map of cell and material IDs. Used for in-memory image
generation.
Parameters
----------
plot : openmc.capi.plot._PlotBase
Object describing the slice of the model to be generated
Returns
-------
id_map : numpy.ndarray
A NumPy array with shape (vertical pixels, horizontal pixels, 2) of
OpenMC property ids with dtype int32
"""
img_data = np.zeros((plot.v_res, plot.h_res, 2),
dtype=np.dtype('int32'))
_dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32)))
return img_data
_dll.openmc_property_map.argtypes = [POINTER(_PlotBase), POINTER(c_double)]
_dll.openmc_property_map.restype = c_int
_dll.openmc_property_map.errcheck = _error_handler
def property_map(plot):
"""
Generate a 2-D map of cell temperatures and material densities. Used for
in-memory image generation.
Parameters
----------
plot : openmc.capi.plot._PlotBase
Object describing the slice of the model to be generated
Returns
-------
property_map : numpy.ndarray
A NumPy array with shape (vertical pixels, horizontal pixels, 2) of
OpenMC property ids with dtype float
"""
prop_data = np.zeros((plot.v_res, plot.h_res, 2))
_dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double)))
return prop_data

5
openmc/cell.py
View file

@ -571,7 +571,10 @@ class Cell(IDManagerMixin):
# Check for other attributes
t = get_text(elem, 'temperature')
if t is not None:
c.temperature = float(t)
if ' ' in t:
c.temperature = [float(t_i) for t_i in t.split()]
else:
c.temperature = float(t)
for key in ('temperature', 'rotation', 'translation'):
value = get_text(elem, key)
if value is not None:

2
openmc/data/angle_distribution.py
View file

@ -115,7 +115,7 @@ class AngleDistribution(EqualityMixin):
Angular distribution
"""
energy = group['energy'].value
energy = group['energy'][()]
data = group['mu']
offsets = data.attrs['offsets']
interpolation = data.attrs['interpolation']

6
openmc/data/correlated.py
View file

@ -210,15 +210,15 @@ class CorrelatedAngleEnergy(AngleEnergy):
interp_data = group['energy'].attrs['interpolation']
energy_breakpoints = interp_data[0, :]
energy_interpolation = interp_data[1, :]
energy = group['energy'].value
energy = group['energy'][()]
offsets = group['energy_out'].attrs['offsets']
interpolation = group['energy_out'].attrs['interpolation']
n_discrete_lines = group['energy_out'].attrs['n_discrete_lines']
dset_eout = group['energy_out'].value
dset_eout = group['energy_out'][()]
energy_out = []
dset_mu = group['mu'].value
dset_mu = group['mu'][()]
mu = []
n_energy = len(energy)

BIN
openmc/data/density_effect.h5 Normal file
View file

Binary file not shown.

81
openmc/data/endf.py
View file

@ -66,7 +66,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):
@ -89,12 +89,15 @@ 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\3', s))
def _int_endf(s):
"""Convert string to int. Used for INTG records where blank entries
indicate a 0.
def int_endf(s):
"""Convert string of integer number in ENDF to int.
The ENDF-6 format technically allows integers to be represented by a field
of all blanks. This function acts like int(s) except when s is a string of
all whitespace, in which case zero is returned.
Parameters
----------
@ -106,8 +109,7 @@ def _int_endf(s):
integer
The number or 0
"""
s = s.strip()
return int(s) if s else 0
return 0 if s.isspace() else int(s)
def get_text_record(file_obj):
@ -127,35 +129,35 @@ def get_text_record(file_obj):
return file_obj.readline()[:66]
def get_cont_record(file_obj, skipC=False):
def get_cont_record(file_obj, skip_c=False):
"""Return data from a CONT record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
skipC : bool
skip_c : bool
Determine whether to skip the first two quantities (C1, C2) of the CONT
record.
Returns
-------
list
tuple
The six items within the CONT record
"""
line = file_obj.readline()
if skipC:
if skip_c:
C1 = None
C2 = None
else:
C1 = float_endf(line[:11])
C2 = float_endf(line[11:22])
L1 = int(line[22:33])
L2 = int(line[33:44])
N1 = int(line[44:55])
N2 = int(line[55:66])
return [C1, C2, L1, L2, N1, N2]
L1 = int_endf(line[22:33])
L2 = int_endf(line[33:44])
N1 = int_endf(line[44:55])
N2 = int_endf(line[55:66])
return (C1, C2, L1, L2, N1, N2)
def get_head_record(file_obj):
@ -168,18 +170,18 @@ def get_head_record(file_obj):
Returns
-------
list
tuple
The six items within the HEAD record
"""
line = file_obj.readline()
ZA = int(float_endf(line[:11]))
AWR = float_endf(line[11:22])
L1 = int(line[22:33])
L2 = int(line[33:44])
N1 = int(line[44:55])
N2 = int(line[55:66])
return [ZA, AWR, L1, L2, N1, N2]
L1 = int_endf(line[22:33])
L2 = int_endf(line[33:44])
N1 = int_endf(line[44:55])
N2 = int_endf(line[55:66])
return (ZA, AWR, L1, L2, N1, N2)
def get_list_record(file_obj):
@ -233,10 +235,10 @@ def get_tab1_record(file_obj):
line = file_obj.readline()
C1 = float_endf(line[:11])
C2 = float_endf(line[11:22])
L1 = int(line[22:33])
L2 = int(line[33:44])
n_regions = int(line[44:55])
n_pairs = int(line[55:66])
L1 = int_endf(line[22:33])
L2 = int_endf(line[33:44])
n_regions = int_endf(line[44:55])
n_pairs = int_endf(line[55:66])
params = [C1, C2, L1, L2]
# Read the interpolation region data, namely NBT and INT
@ -247,8 +249,8 @@ def get_tab1_record(file_obj):
line = file_obj.readline()
to_read = min(3, n_regions - m)
for j in range(to_read):
breakpoints[m] = int(line[0:11])
interpolation[m] = int(line[11:22])
breakpoints[m] = int_endf(line[0:11])
interpolation[m] = int_endf(line[11:22])
line = line[22:]
m += 1
@ -306,9 +308,9 @@ def get_intg_record(file_obj):
"""
# determine how many items are in list and NDIGIT
items = get_cont_record(file_obj)
ndigit = int(items[2])
npar = int(items[3]) # Number of parameters
nlines = int(items[4]) # Lines to read
ndigit = items[2]
npar = items[3] # Number of parameters
nlines = items[4] # Lines to read
NROW_RULES = {2: 18, 3: 12, 4: 11, 5: 9, 6: 8}
nrow = NROW_RULES[ndigit]
@ -316,13 +318,13 @@ def get_intg_record(file_obj):
corr = np.identity(npar)
for i in range(nlines):
line = file_obj.readline()
ii = _int_endf(line[:5]) - 1 # -1 to account for 0 indexing
jj = _int_endf(line[5:10]) - 1
ii = int_endf(line[:5]) - 1 # -1 to account for 0 indexing
jj = int_endf(line[5:10]) - 1
factor = 10**ndigit
for j in range(nrow):
if jj+j >= ii:
break
element = _int_endf(line[11+(ndigit+1)*j:11+(ndigit+1)*(j+1)])
element = int_endf(line[11+(ndigit+1)*j:11+(ndigit+1)*(j+1)])
if element > 0:
corr[ii, jj] = (element+0.5)/factor
elif element < 0:
@ -507,16 +509,7 @@ class Evaluation(object):
# File numbers, reaction designations, and number of records
for i in range(NXC):
line = file_obj.readline()
mf = int(line[22:33])
mt = int(line[33:44])
nc = int(line[44:55])
try:
mod = int(line[55:66])
except ValueError:
# In JEFF 3.2, a few isotopes of U have MOD values that are
# missing. This prevents failure on these isotopes.
mod = 0
_, _, mf, mt, nc, mod = get_cont_record(file_obj, skip_c=True)
self.reaction_list.append((mf, mt, nc, mod))
@property

2
openmc/data/energy_distribution.py
View file

@ -1144,7 +1144,7 @@ class ContinuousTabular(EnergyDistribution):
interp_data = group['energy'].attrs['interpolation']
energy_breakpoints = interp_data[0, :]
energy_interpolation = interp_data[1, :]
energy = group['energy'].value
energy = group['energy'][()]
data = group['distribution']
offsets = data.attrs['offsets']

6
openmc/data/function.py
View file

@ -349,8 +349,8 @@ class Tabulated1D(Function1D):
raise ValueError("Expected an HDF5 attribute 'type' equal to '"
+ cls.__name__ + "'")
x = dataset.value[0, :]
y = dataset.value[1, :]
x = dataset[0, :]
y = dataset[1, :]
breakpoints = dataset.attrs['breakpoints']
interpolation = dataset.attrs['interpolation']
return cls(x, y, breakpoints, interpolation)
@ -434,7 +434,7 @@ class Polynomial(np.polynomial.Polynomial, Function1D):
if dataset.attrs['type'].decode() != cls.__name__:
raise ValueError("Expected an HDF5 attribute 'type' equal to '"
+ cls.__name__ + "'")
return cls(dataset.value)
return cls(dataset[()])
class Combination(EqualityMixin):

2
openmc/data/kalbach_mann.py
View file

@ -202,7 +202,7 @@ class KalbachMann(AngleEnergy):
interp_data = group['energy'].attrs['interpolation']
energy_breakpoints = interp_data[0, :]
energy_interpolation = interp_data[1, :]
energy = group['energy'].value
energy = group['energy'][()]
data = group['distribution']
offsets = data.attrs['offsets']

16
openmc/data/multipole.py
View file

@ -356,24 +356,24 @@ class WindowedMultipole(EqualityMixin):
# Read scalars.
out.spacing = group['spacing'].value
out.sqrtAWR = group['sqrtAWR'].value
out.E_min = group['E_min'].value
out.E_max = group['E_max'].value
out.spacing = group['spacing'][()]
out.sqrtAWR = group['sqrtAWR'][()]
out.E_min = group['E_min'][()]
out.E_max = group['E_max'][()]
# Read arrays.
err = "WMP '{}' array shape is not consistent with the '{}' array shape"
out.data = group['data'].value
out.data = group['data'][()]
out.windows = group['windows'].value
out.windows = group['windows'][()]
out.broaden_poly = group['broaden_poly'].value.astype(np.bool)
out.broaden_poly = group['broaden_poly'][...].astype(np.bool)
if out.broaden_poly.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('broaden_poly', 'windows'))
out.curvefit = group['curvefit'].value
out.curvefit = group['curvefit'][()]
if out.curvefit.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('curvefit', 'windows'))

43
openmc/data/neutron.py
View file

@ -448,11 +448,13 @@ class IncidentNeutron(EqualityMixin):
for rx in self.reactions.values():
# Skip writing redundant reaction if it doesn't have photon
# production or is a summed transmutation reaction. MT=4 is also
# sometimes needed for probability tables.
# sometimes needed for probability tables. Also write gas
# production, heating, and damage energy production.
if rx.redundant:
photon_rx = any(p.particle == 'photon' for p in rx.products)
transmutation_rx = (rx.mt in (16, 103, 104, 105, 106, 107))
if not (photon_rx or transmutation_rx or rx.mt == 4):
keep_mts = (4, 16, 103, 104, 105, 106, 107,
203, 204, 205, 206, 207, 301, 444)
if not (photon_rx or rx.mt in keep_mts):
continue
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
@ -519,7 +521,7 @@ class IncidentNeutron(EqualityMixin):
kTg = group['kTs']
kTs = []
for temp in kTg:
kTs.append(kTg[temp].value)
kTs.append(kTg[temp][()])
data = cls(name, atomic_number, mass_number, metastable,
atomic_weight_ratio, kTs)
@ -527,7 +529,7 @@ class IncidentNeutron(EqualityMixin):
# Read energy grid
e_group = group['energy']
for temperature, dset in e_group.items():
data.energy[temperature] = dset.value
data.energy[temperature] = dset[()]
# Read reaction data
rxs_group = group['reactions']
@ -615,13 +617,14 @@ class IncidentNeutron(EqualityMixin):
# Read energy grid
n_energy = ace.nxs[3]
energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy]*EV_PER_MEV
i = ace.jxs[1]
energy = ace.xss[i : i + n_energy]*EV_PER_MEV
data.energy[strT] = energy
total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2 * n_energy]
absorption_xs = ace.xss[ace.jxs[1] + 2 * n_energy:ace.jxs[1] +
3 * n_energy]
total_xs = ace.xss[i + n_energy : i + 2*n_energy]
absorption_xs = ace.xss[i + 2*n_energy : i + 3*n_energy]
heating_number = ace.xss[i + 4*n_energy : i + 5*n_energy]*EV_PER_MEV
# Create redundant reactions (total and absorption)
# Create redundant reactions (total, absorption, and heating)
total = Reaction(1)
total.xs[strT] = Tabulated1D(energy, total_xs)
total.redundant = True
@ -633,13 +636,15 @@ class IncidentNeutron(EqualityMixin):
absorption.redundant = True
data.reactions[101] = absorption
heating = Reaction(301)
heating.xs[strT] = Tabulated1D(energy, heating_number*total_xs)
heating.redundant = True
data.reactions[301] = heating
# Read each reaction
n_reaction = ace.nxs[4] + 1
for i in range(n_reaction):
rx = Reaction.from_ace(ace, i)
# Don't include gas production / damage cross sections
if 200 < rx.mt < 219 or rx.mt == 444:
continue
data.reactions[rx.mt] = rx
# Some photon production reactions may be assigned to MTs that don't
@ -690,6 +695,8 @@ class IncidentNeutron(EqualityMixin):
mts = data.get_reaction_components(rx.mt)
if mts != [rx.mt]:
rx.redundant = True
if rx.mt in (203, 204, 205, 206, 207, 444):
rx.redundant = True
# Read unresolved resonance probability tables
urr = ProbabilityTables.from_ace(ace)
@ -784,16 +791,19 @@ class IncidentNeutron(EqualityMixin):
return data
@classmethod
def from_njoy(cls, filename, temperatures=None, **kwargs):
def from_njoy(cls, filename, temperatures=None, evaluation=None, **kwargs):
"""Generate incident neutron data by running NJOY.
Parameters
----------
filename : str
Path to ENDF evaluation
Path to ENDF file
temperatures : iterable of float
Temperatures in Kelvin to produce data at. If omitted, data is
produced at room temperature (293.6 K)
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this
argument indicates which evaluation to use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace`
@ -808,6 +818,7 @@ class IncidentNeutron(EqualityMixin):
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
pendf_file = os.path.join(tmpdir, 'pendf')
kwargs['evaluation'] = evaluation
make_ace(filename, temperatures, ace_file, xsdir_file,
pendf_file, **kwargs)
@ -818,7 +829,7 @@ class IncidentNeutron(EqualityMixin):
data.add_temperature_from_ace(table)
# Add fission energy release data
ev = Evaluation(filename)
ev = evaluation if evaluation is not None else Evaluation(filename)
if (1, 458) in ev.section:
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)

43
openmc/data/njoy.py
View file

@ -72,8 +72,13 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr} /
{mat} 3 /
302 318 402 /
{mat} 4 /
302 318 402 444 /
"""
_TEMPLATE_GASPR = """
gaspr / %%%%%%%%%%%%%%%%%%%%%%%%% Add gas production %%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {ngaspr_in} {ngaspr} /
"""
_TEMPLATE_PURR = """
@ -186,7 +191,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
"""Generate ACE file from an ENDF file
"""Generate pointwise ENDF file from an ENDF file
Parameters
----------
@ -211,8 +216,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
error=0.001, broadr=True, heatr=True, purr=True, acer=True,
**kwargs):
error=0.001, broadr=True, heatr=True, gaspr=True, purr=True,
acer=True, evaluation=None, **kwargs):
"""Generate incident neutron ACE file from an ENDF file
Parameters
@ -234,10 +239,15 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
Indicating whether to Doppler broaden XS when running NJOY
heatr : bool, optional
Indicating whether to add heating kerma when running NJOY
gaspr : bool, optional
Indicating whether to add gas production data when running NJOY
purr : bool, optional
Indicating whether to add probability table when running NJOY
acer : bool, optional
Indicating whether to generate ACE file when running NJOY
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this argument
indicates which evaluation should be used.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -247,7 +257,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
If the NJOY process returns with a non-zero status
"""
ev = endf.Evaluation(filename)
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
@ -285,6 +295,13 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
commands += _TEMPLATE_HEATR
nlast = nheatr
# gaspr
if gaspr:
ngaspr_in = nlast
ngaspr = ngaspr_in + 1
commands += _TEMPLATE_GASPR
nlast = ngaspr
# purr
if purr:
npurr_in = nlast
@ -338,7 +355,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir='xsdir', error=0.001, **kwargs):
ace='ace', xsdir='xsdir', error=0.001, evaluation=None,
evaluation_thermal=None, **kwargs):
"""Generate thermal scattering ACE file from ENDF files
Parameters
@ -356,6 +374,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
Path of xsdir file to write
error : float, optional
Fractional error tolerance for NJOY processing
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF neutron sublibrary file contains multiple material
evaluations, this argument indicates which evaluation to use.
evaluation_thermal : openmc.data.endf.Evaluation, optional
If the ENDF thermal scattering sublibrary file contains multiple
material evaluations, this argument indicates which evaluation to use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -365,11 +389,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
If the NJOY process returns with a non-zero status
"""
ev = endf.Evaluation(filename)
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
ev_thermal = endf.Evaluation(filename_thermal)
ev_thermal = (evaluation_thermal if evaluation_thermal is not None
else endf.Evaluation(filename_thermal))
mat_thermal = ev_thermal.material
zsymam_thermal = ev_thermal.target['zsymam']

678
openmc/data/photon.py
View file

@ -19,71 +19,62 @@ from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
from .function import Tabulated1D
_SUBSHELLS = ['K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2',
'O3', 'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2',
'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11',
'Q1', 'Q2', 'Q3']
# Helper function to map designator to subshell string or None
def _subshell(i):
if i == 0:
return None
else:
return _SUBSHELLS[i - 1]
# Electron subshell labels
_SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3',
'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2', 'P3', 'P4',
'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11','Q1', 'Q2', 'Q3']
_REACTION_NAME = {
501: 'Total photon interaction',
502: 'Photon coherent scattering',
504: 'Photon incoherent scattering',
515: 'Pair production, electron field',
516: 'Total pair production',
517: 'Pair production, nuclear field',
522: 'Photoelectric absorption',
526: 'Electro-atomic scattering',
527: 'Electro-atomic bremsstrahlung',
528: 'Electro-atomic excitation',
534: 'K (1s1/2) subshell photoelectric',
535: 'L1 (2s1/2) subshell photoelectric',
536: 'L2 (2p1/2) subshell photoelectric',
537: 'L3 (2p3/2) subshell photoelectric',
538: 'M1 (3s1/2) subshell photoelectric',
539: 'M2 (3p1/2) subshell photoelectric',
540: 'M3 (3p3/2) subshell photoelectric',
541: 'M4 (3d3/2) subshell photoelectric',
542: 'M5 (3d5/2) subshell photoelectric',
543: 'N1 (4s1/2) subshell photoelectric',
544: 'N2 (4p1/2) subshell photoelectric',
545: 'N3 (4p3/2) subshell photoelectric',
546: 'N4 (4d3/2) subshell photoelectric',
547: 'N5 (4d5/2) subshell photoelectric',
548: 'N6 (4f5/2) subshell photoelectric',
549: 'N7 (4f7/2) subshell photoelectric',
550: 'O1 (5s1/2) subshell photoelectric',
551: 'O2 (5p1/2) subshell photoelectric',
552: 'O3 (5p3/2) subshell photoelectric',
553: 'O4 (5d3/2) subshell photoelectric',
554: 'O5 (5d5/2) subshell photoelectric',
555: 'O6 (5f5/2) subshell photoelectric',
556: 'O7 (5f7/2) subshell photoelectric',
557: 'O8 (5g7/2) subshell photoelectric',
558: 'O9 (5g9/2) subshell photoelectric',
559: 'P1 (6s1/2) subshell photoelectric',
560: 'P2 (6p1/2) subshell photoelectric',
561: 'P3 (6p3/2) subshell photoelectric',
562: 'P4 (6d3/2) subshell photoelectric',
563: 'P5 (6d5/2) subshell photoelectric',
564: 'P6 (6f5/2) subshell photoelectric',
565: 'P7 (6f7/2) subshell photoelectric',
566: 'P8 (6g7/2) subshell photoelectric',
567: 'P9 (6g9/2) subshell photoelectric',
568: 'P10 (6h9/2) subshell photoelectric',
569: 'P11 (6h11/2) subshell photoelectric',
570: 'Q1 (7s1/2) subshell photoelectric',
571: 'Q2 (7p1/2) subshell photoelectric',
572: 'Q3 (7p3/2) subshell photoelectric'
501: ('Total photon interaction', 'total'),
502: ('Photon coherent scattering', 'coherent'),
504: ('Photon incoherent scattering', 'incoherent'),
515: ('Pair production, electron field', 'pair_production_electron'),
516: ('Total pair production', 'pair_production_total'),
517: ('Pair production, nuclear field', 'pair_production_nuclear'),
522: ('Photoelectric absorption', 'photoelectric'),
526: ('Electro-atomic scattering', 'electro_atomic_scat'),
527: ('Electro-atomic bremsstrahlung', 'electro_atomic_brem'),
528: ('Electro-atomic excitation', 'electro_atomic_excit'),
534: ('K (1s1/2) subshell photoelectric', 'K'),
535: ('L1 (2s1/2) subshell photoelectric', 'L1'),
536: ('L2 (2p1/2) subshell photoelectric', 'L2'),
537: ('L3 (2p3/2) subshell photoelectric', 'L3'),
538: ('M1 (3s1/2) subshell photoelectric', 'M1'),
539: ('M2 (3p1/2) subshell photoelectric', 'M2'),
540: ('M3 (3p3/2) subshell photoelectric', 'M3'),
541: ('M4 (3d3/2) subshell photoelectric', 'M4'),
542: ('M5 (3d5/2) subshell photoelectric', 'M5'),
543: ('N1 (4s1/2) subshell photoelectric', 'N1'),
544: ('N2 (4p1/2) subshell photoelectric', 'N2'),
545: ('N3 (4p3/2) subshell photoelectric', 'N3'),
546: ('N4 (4d3/2) subshell photoelectric', 'N4'),
547: ('N5 (4d5/2) subshell photoelectric', 'N5'),
548: ('N6 (4f5/2) subshell photoelectric', 'N6'),
549: ('N7 (4f7/2) subshell photoelectric', 'N7'),
550: ('O1 (5s1/2) subshell photoelectric', 'O1'),
551: ('O2 (5p1/2) subshell photoelectric', 'O2'),
552: ('O3 (5p3/2) subshell photoelectric', 'O3'),
553: ('O4 (5d3/2) subshell photoelectric', 'O4'),
554: ('O5 (5d5/2) subshell photoelectric', 'O5'),
555: ('O6 (5f5/2) subshell photoelectric', 'O6'),
556: ('O7 (5f7/2) subshell photoelectric', 'O7'),
557: ('O8 (5g7/2) subshell photoelectric', 'O8'),
558: ('O9 (5g9/2) subshell photoelectric', 'O9'),
559: ('P1 (6s1/2) subshell photoelectric', 'P1'),
560: ('P2 (6p1/2) subshell photoelectric', 'P2'),
561: ('P3 (6p3/2) subshell photoelectric', 'P3'),
562: ('P4 (6d3/2) subshell photoelectric', 'P4'),
563: ('P5 (6d5/2) subshell photoelectric', 'P5'),
564: ('P6 (6f5/2) subshell photoelectric', 'P6'),
565: ('P7 (6f7/2) subshell photoelectric', 'P7'),
566: ('P8 (6g7/2) subshell photoelectric', 'P8'),
567: ('P9 (6g9/2) subshell photoelectric', 'P9'),
568: ('P10 (6h9/2) subshell photoelectric', 'P10'),
569: ('P11 (6h11/2) subshell photoelectric', 'P11'),
570: ('Q1 (7s1/2) subshell photoelectric', 'Q1'),
571: ('Q2 (7p1/2) subshell photoelectric', 'Q2'),
572: ('Q3 (7p3/2) subshell photoelectric', 'Q3')
}
# Compton profiles are read from a pre-generated HDF5 file when they are first
@ -92,18 +83,15 @@ _REACTION_NAME = {
# is a 2D array with shape (n_shells, n_momentum_values) stored on the key Z
_COMPTON_PROFILES = {}
# Stopping powers are read from a pre-generated HDF5 file when they are first
# needed. The dictionary stores an array of energy values at which the other
# quantities are tabulated with the key 'energy' and for each element has the
# mean excitation energy and arrays containing the collision stopping powers
# and radiative stopping powers stored on the key 'Z'.
_STOPPING_POWERS = {}
# Scaled bremsstrahlung DCSs are read from a data file provided by Selzter and
# Berger when they are first needed. The dictionary stores an array of n
# incident electron kinetic energies with key 'electron_energies', an array of
# k reduced photon energies with key 'photon_energies', and the cross sections
# for each element are in a 2D array with shape (n, k) stored on the key 'Z'.
# It also stores data used for calculating the density effect correction and
# stopping power, namely, the mean excitation energy with the key 'I', number
# of electrons per subshell with the key 'num_electrons', and binding energies
# with the key 'ionization_energy'.
_BREMSSTRAHLUNG = {}
@ -228,7 +216,7 @@ class AtomicRelaxation(EqualityMixin):
# Get shell designators
n = ace.nxs[7]
idx = ace.jxs[11]
shells = [_subshell(int(i)) for i in ace.xss[idx : idx+n]]
shells = [_SUBSHELLS[int(i)] for i in ace.xss[idx : idx+n]]
# Get number of electrons for each shell
idx = ace.jxs[12]
@ -248,8 +236,8 @@ class AtomicRelaxation(EqualityMixin):
if n_transitions > 0:
records = []
for j in range(n_transitions):
subj = _subshell(int(ace.xss[idx]))
subk = _subshell(int(ace.xss[idx + 1]))
subj = _SUBSHELLS[int(ace.xss[idx])]
subk = _SUBSHELLS[int(ace.xss[idx + 1])]
etr = ace.xss[idx + 2]*EV_PER_MEV
if j == 0:
ftr = ace.xss[idx + 3]
@ -304,7 +292,7 @@ class AtomicRelaxation(EqualityMixin):
# Read data for each subshell
for i in range(n_subshells):
params, list_items = get_list_record(file_obj)
subi = _subshell(int(params[0]))
subi = _SUBSHELLS[int(params[0])]
n_transitions = int(params[5])
binding_energy[subi] = list_items[0]
num_electrons[subi] = list_items[1]
@ -313,8 +301,8 @@ class AtomicRelaxation(EqualityMixin):
# Read transition data
records = []
for j in range(n_transitions):
subj = _subshell(int(list_items[6*(j+1)]))
subk = _subshell(int(list_items[6*(j+1) + 1]))
subj = _SUBSHELLS[int(list_items[6*(j+1)])]
subk = _SUBSHELLS[int(list_items[6*(j+1) + 1])]
etr = list_items[6*(j+1) + 2]
ftr = list_items[6*(j+1) + 3]
records.append((subj, subk, etr, ftr))
@ -326,8 +314,70 @@ class AtomicRelaxation(EqualityMixin):
# Return instance of class
return cls(binding_energy, num_electrons, transitions)
def to_hdf5(self, group):
raise NotImplementedError
@classmethod
def from_hdf5(cls, group):
"""Generate atomic relaxation data from an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to read from
Returns
-------
openmc.data.AtomicRelaxation
Atomic relaxation data
"""
# Create data dictionaries
binding_energy = {}
num_electrons = {}
transitions = {}
designators = [s.decode() for s in group.attrs['designators']]
columns = ['secondary', 'tertiary', 'energy (eV)', 'probability']
for shell in designators:
# Shell group
sub_group = group[shell]
# Read subshell binding energy and number of electrons
if 'binding_energy' in sub_group.attrs:
binding_energy[shell] = sub_group.attrs['binding_energy']
if 'num_electrons' in sub_group.attrs:
num_electrons[shell] = sub_group.attrs['num_electrons']
# Read transition data
if 'transitions' in sub_group:
df = pd.DataFrame(sub_group['transitions'][()],
columns=columns)
# Replace float indexes back to subshell strings
df[columns[:2]] = df[columns[:2]].replace(
np.arange(float(len(_SUBSHELLS))), _SUBSHELLS)
transitions[shell] = df
return cls(binding_energy, num_electrons, transitions)
def to_hdf5(self, group, shell):
"""Write atomic relaxation data to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
shell : str
The subshell to write data for
"""
# Write subshell binding energy and number of electrons
group.attrs['binding_energy'] = self.binding_energy[shell]
group.attrs['num_electrons'] = self.num_electrons[shell]
# Write transition data with replacements
if shell in self.transitions:
df = self.transitions[shell].replace(
_SUBSHELLS, range(len(_SUBSHELLS)))
group.create_dataset('transitions', data=df.values.astype(float))
class IncidentPhoton(EqualityMixin):
@ -352,12 +402,16 @@ class IncidentPhoton(EqualityMixin):
atomic_relaxation : openmc.data.AtomicRelaxation or None
Atomic relaxation data
bremsstrahlung : dict
Dictionary of bremsstrahlung DCS data with keys 'electron_energy'
(incident electron kinetic energy values in [eV]), 'photon_energy'
(ratio of the energy of the emitted photon to the incident electron
kinetic energy), and 'dcs' (cross section values in [b]). The cross
sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`,
where :math:`E_k` is the energy of the emitted photon.
Dictionary of bremsstrahlung data with keys 'I' (mean excitation energy
in [eV]), 'num_electrons' (number of electrons in each subshell),
'ionization_energy' (ionization potential of each subshell),
'electron_energy' (incident electron kinetic energy values in [eV]),
'photon_energy' (ratio of the energy of the emitted photon to the
incident electron kinetic energy), and 'dcs' (cross section values in
[b]). The cross sections are in scaled form: :math:`(\beta^2/Z^2) E_k
(d\sigma/dE_k)`, where :math:`E_k` is the energy of the emitted photon.
A negative number of electrons in a subshell indicates conduction
electrons.
compton_profiles : dict
Dictionary of Compton profile data with keys 'num_electrons' (number of
electrons in each subshell), 'binding_energy' (ionization potential of
@ -368,11 +422,6 @@ class IncidentPhoton(EqualityMixin):
reactions : collections.OrderedDict
Contains the cross sections for each photon reaction. The keys are MT
values and the values are instances of :class:`PhotonReaction`.
stopping_powers : dict
Dictionary of stopping power data with keys 'energy' (in [eV]), 'I' (mean
excitation energy), 's_collision' (collision stopping power in
[eV cm\ :sup:`2`/g]), and 's_radiative' (radiative stopping power in
[eV cm\ :sup:`2`/g])
"""
@ -381,7 +430,6 @@ class IncidentPhoton(EqualityMixin):
self._atomic_relaxation = None
self.reactions = OrderedDict()
self.compton_profiles = {}
self.stopping_powers = {}
self.bremsstrahlung = {}
def __contains__(self, mt):
@ -514,10 +562,13 @@ class IncidentPhoton(EqualityMixin):
idx += n_energy
# Copy binding energy
shell = _subshell(d)
shell = _SUBSHELLS[d]
e = data.atomic_relaxation.binding_energy[shell]
rx.subshell_binding_energy = e
# Add bremsstrahlung DCS data
data._add_bremsstrahlung()
return data
@classmethod
@ -560,12 +611,12 @@ class IncidentPhoton(EqualityMixin):
if not _COMPTON_PROFILES:
filename = os.path.join(os.path.dirname(__file__), 'compton_profiles.h5')
with h5py.File(filename, 'r') as f:
_COMPTON_PROFILES['pz'] = f['pz'].value
_COMPTON_PROFILES['pz'] = f['pz'][()]
for i in range(1, 101):
group = f['{:03}'.format(i)]
num_electrons = group['num_electrons'].value
binding_energy = group['binding_energy'].value*EV_PER_MEV
J = group['J'].value
num_electrons = group['num_electrons'][()]
binding_energy = group['binding_energy'][()]*EV_PER_MEV
J = group['J'][()]
_COMPTON_PROFILES[i] = {'num_electrons': num_electrons,
'binding_energy': binding_energy,
'J': J}
@ -577,29 +628,189 @@ class IncidentPhoton(EqualityMixin):
data.compton_profiles['binding_energy'] = profile['binding_energy']
data.compton_profiles['J'] = [Tabulated1D(pz, J_k) for J_k in profile['J']]
# Load stopping power data if it has not yet been loaded
if not _STOPPING_POWERS:
filename = os.path.join(os.path.dirname(__file__), 'stopping_powers.h5')
with h5py.File(filename, 'r') as f:
# Units are in MeV; convert to eV
_STOPPING_POWERS['energy'] = f['energy'].value*EV_PER_MEV
for i in range(1, 99):
group = f['{:03}'.format(i)]
# Add bremsstrahlung DCS data
data._add_bremsstrahlung()
# Units are in MeV cm^2/g; convert to eV cm^2/g
_STOPPING_POWERS[i] = {
'I': group.attrs['I'],
's_collision': group['s_collision'].value*EV_PER_MEV,
's_radiative': group['s_radiative'].value*EV_PER_MEV
}
return data
# Add stopping power data
if Z < 99:
data.stopping_powers['energy'] = _STOPPING_POWERS['energy']
data.stopping_powers.update(_STOPPING_POWERS[Z])
@classmethod
def from_hdf5(cls, group_or_filename):
"""Generate photon reaction from an HDF5 group
Parameters
----------
group_or_filename : h5py.Group or str
HDF5 group containing interaction data. If given as a string, it is
assumed to be the filename for the HDF5 file, and the first group is
used to read from.
Returns
-------
openmc.data.IncidentPhoton
Photon interaction data
"""
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
h5file = h5py.File(str(group_or_filename), 'r')
# Make sure version matches
if 'version' in h5file.attrs:
major, minor = h5file.attrs['version']
# For now all versions of HDF5 data can be read
else:
raise IOError(
'HDF5 data does not indicate a version. Your installation '
'of the OpenMC Python API expects version {}.x data.'
.format(HDF5_VERSION_MAJOR))
group = list(h5file.values())[0]
Z = group.attrs['Z']
data = cls(Z)
# Read energy grid
energy = group['energy'][()]
# Read cross section data
for mt, (name, key) in _REACTION_NAME.items():
if key in group:
rgroup = group[key]
elif key in group['subshells']:
rgroup = group['subshells'][key]
else:
continue
data.reactions[mt] = PhotonReaction.from_hdf5(rgroup, mt, energy)
# Check for necessary reactions
for mt in (502, 504, 522):
assert mt in data, "Reaction {} not found".format(mt)
# Read atomic relaxation
data.atomic_relaxation = AtomicRelaxation.from_hdf5(group['subshells'])
# Read Compton profiles
if 'compton_profiles' in group:
rgroup = group['compton_profiles']
profile = data.compton_profiles
profile['num_electrons'] = rgroup['num_electrons'][()]
profile['binding_energy'] = rgroup['binding_energy'][()]
# Get electron momentum values
pz = rgroup['pz'][()]
J = rgroup['J'][()]
if pz.size != J.shape[1]:
raise ValueError("'J' array shape is not consistent with the "
"'pz' array shape")
profile['J'] = [Tabulated1D(pz, Jk) for Jk in J]
# Read bremsstrahlung
if 'bremsstrahlung' in group:
rgroup = group['bremsstrahlung']
data.bremsstrahlung['I'] = rgroup.attrs['I']
for key in ('dcs', 'electron_energy', 'ionization_energy',
'num_electrons', 'photon_energy'):
data.bremsstrahlung[key] = rgroup[key][()]
return data
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export incident photon data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
mode : {'r', r+', 'w', 'x', 'a'}
Mode that is used to open the HDF5 file. This is the second argument
to the :class:`h5py.File` constructor.
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
f = h5py.File(str(path), mode, libver=libver)
f.attrs['filetype'] = np.string_('data_photon')
if 'version' not in f.attrs:
f.attrs['version'] = np.array(HDF5_VERSION)
group = f.create_group(self.name)
group.attrs['Z'] = Z = self.atomic_number
# Determine union energy grid
union_grid = np.array([])
for rx in self:
union_grid = np.union1d(union_grid, rx.xs.x)
group.create_dataset('energy', data=union_grid)
# Write cross sections
shell_group = group.create_group('subshells')
designators = []
for mt, rx in self.reactions.items():
name, key = _REACTION_NAME[mt]
if mt in [502, 504, 515, 517, 522]:
sub_group = group.create_group(key)
elif mt >= 534 and mt <= 572:
# Subshell
designators.append(key)
sub_group = shell_group.create_group(key)
# Write atomic relaxation
if key in self.atomic_relaxation.subshells:
self.atomic_relaxation.to_hdf5(sub_group, key)
else:
continue
rx.to_hdf5(sub_group, union_grid, Z)
shell_group.attrs['designators'] = np.array(designators, dtype='S')
# Write Compton profiles
if self.compton_profiles:
compton_group = group.create_group('compton_profiles')
profile = self.compton_profiles
compton_group.create_dataset('num_electrons',
data=profile['num_electrons'])
compton_group.create_dataset('binding_energy',
data=profile['binding_energy'])
# Get electron momentum values
compton_group.create_dataset('pz', data=profile['J'][0].x)
# Create/write 2D array of profiles
J = np.array([Jk.y for Jk in profile['J']])
compton_group.create_dataset('J', data=J)
# Write bremsstrahlung
if self.bremsstrahlung:
brem_group = group.create_group('bremsstrahlung')
for key, value in self.bremsstrahlung.items():
if key == 'I':
brem_group.attrs[key] = value
else:
brem_group.create_dataset(key, data=value)
def _add_bremsstrahlung(self):
"""Add the data used in the thick-target bremsstrahlung approximation
"""
# Load bremsstrahlung data if it has not yet been loaded
if not _BREMSSTRAHLUNG:
# Add data used for density effect correction
filename = os.path.join(os.path.dirname(__file__), 'density_effect.h5')
with h5py.File(filename, 'r') as f:
for i in range(1, 101):
group = f['{:03}'.format(i)]
_BREMSSTRAHLUNG[i] = {
'I': group.attrs['I'],
'num_electrons': group['num_electrons'][()],
'ionization_energy': group['ionization_energy'][()]
}
filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT')
brem = open(filename, 'r').read().split()
@ -640,153 +851,12 @@ class IncidentPhoton(EqualityMixin):
# Get scaled DCS values (millibarns) on new energy grid
dcs[:,j] = cs(log_energy)
_BREMSSTRAHLUNG[i] = {'dcs': dcs}
_BREMSSTRAHLUNG[i]['dcs'] = dcs
# Add bremsstrahlung DCS data
data.bremsstrahlung['electron_energy'] = _BREMSSTRAHLUNG['electron_energy']
data.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy']
data.bremsstrahlung['dcs'] = _BREMSSTRAHLUNG[Z]['dcs']
return data
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export incident photon data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
mode : {'r', r+', 'w', 'x', 'a'}
Mode that is used to open the HDF5 file. This is the second argument
to the :class:`h5py.File` constructor.
"""
# Open file and write version
f = h5py.File(str(path), mode, libver=libver)
f.attrs['filetype'] = np.string_('data_photon')
if 'version' not in f.attrs:
f.attrs['version'] = np.array(HDF5_VERSION)
group = f.create_group(self.name)
group.attrs['Z'] = Z = self.atomic_number
# Determine union energy grid
union_grid = np.array([])
for rx in self:
union_grid = np.union1d(union_grid, rx.xs.x)
group.create_dataset('energy', data=union_grid)
# Write coherent scattering cross section
rx = self.reactions[502]
coh_group = group.create_group('coherent')
coh_group.create_dataset('xs', data=rx.xs(union_grid))
if rx.scattering_factor is not None:
# Create integrated form factor
ff = deepcopy(rx.scattering_factor)
ff.x *= ff.x
ff.y *= ff.y/Z**2
int_ff = Tabulated1D(ff.x, ff.integral())
int_ff.to_hdf5(coh_group, 'integrated_scattering_factor')
if rx.anomalous_real is not None:
rx.anomalous_real.to_hdf5(coh_group, 'anomalous_real')
if rx.anomalous_imag is not None:
rx.anomalous_imag.to_hdf5(coh_group, 'anomalous_imag')
# Write incoherent scattering cross section
rx = self[504]
incoh_group = group.create_group('incoherent')
incoh_group.create_dataset('xs', data=rx.xs(union_grid))
if rx.scattering_factor is not None:
rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor')
# Write electron-field pair production cross section
if 515 in self:
pair_group = group.create_group('pair_production_electron')
pair_group.create_dataset('xs', data=self[515].xs(union_grid))
# Write nuclear-field pair production cross section
if 517 in self:
pair_group = group.create_group('pair_production_nuclear')
pair_group.create_dataset('xs', data=self[517].xs(union_grid))
# Write photoelectric cross section
photoelec_group = group.create_group('photoelectric')
photoelec_group.create_dataset('xs', data=self[522].xs(union_grid))
# Write photoionization cross sections
shell_group = group.create_group('subshells')
designators = []
for mt, rx in self.reactions.items():
if mt >= 534 and mt <= 572:
# Get name of subshell
shell = _SUBSHELLS[mt - 534]
designators.append(shell)
sub_group = shell_group.create_group(shell)
if self.atomic_relaxation is not None:
relax = self.atomic_relaxation
# Write subshell binding energy and number of electrons
sub_group.attrs['binding_energy'] = relax.binding_energy[shell]
sub_group.attrs['num_electrons'] = relax.num_electrons[shell]
# Write transition data with replacements
if shell in relax.transitions:
shell_values = _SUBSHELLS.copy()
shell_values.insert(0, None)
df = relax.transitions[shell].replace(
shell_values, range(len(shell_values)))
sub_group.create_dataset(
'transitions', data=df.values.astype(float))
# Determine threshold
threshold = rx.xs.x[0]
idx = np.searchsorted(union_grid, threshold, side='right') - 1
# Interpolate cross section onto union grid and write
photoionization = rx.xs(union_grid[idx:])
sub_group.create_dataset('xs', data=photoionization)
assert len(union_grid) == len(photoionization) + idx
sub_group['xs'].attrs['threshold_idx'] = idx
shell_group.attrs['designators'] = np.array(designators, dtype='S')
# Write Compton profiles
if self.compton_profiles:
compton_group = group.create_group('compton_profiles')
profile = self.compton_profiles
compton_group.create_dataset('num_electrons',
data=profile['num_electrons'])
compton_group.create_dataset('binding_energy',
data=profile['binding_energy'])
# Get electron momentum values
compton_group.create_dataset('pz', data=profile['J'][0].x)
# Create/write 2D array of profiles
J = np.array([Jk.y for Jk in profile['J']])
compton_group.create_dataset('J', data=J)
# Write stopping powers
if self.stopping_powers:
s_group = group.create_group('stopping_powers')
for key, value in self.stopping_powers.items():
if key == 'I':
s_group.attrs[key] = value
else:
s_group.create_dataset(key, data=value)
# Write bremsstrahlung
if self.bremsstrahlung:
brem_group = group.create_group('bremsstrahlung')
brem = self.bremsstrahlung
brem_group.create_dataset('electron_energy',
data=brem['electron_energy'])
brem_group.create_dataset('photon_energy',
data=brem['photon_energy'])
brem_group.create_dataset('dcs', data=brem['dcs'])
self.bremsstrahlung['electron_energy'] = _BREMSSTRAHLUNG['electron_energy']
self.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy']
self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number])
class PhotonReaction(EqualityMixin):
@ -822,7 +892,7 @@ class PhotonReaction(EqualityMixin):
def __repr__(self):
if self.mt in _REACTION_NAME:
return "<Photon Reaction: MT={} {}>".format(
self.mt, _REACTION_NAME[self.mt])
self.mt, _REACTION_NAME[self.mt][0])
else:
return "<Photon Reaction: MT={}>".format(self.mt)
@ -999,3 +1069,93 @@ class PhotonReaction(EqualityMixin):
params, rx.anomalous_imag = get_tab1_record(file_obj)
return rx
@classmethod
def from_hdf5(cls, group, mt, energy):
"""Generate photon reaction from an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to read from
mt : int
The MT value of the reaction to get data for
energy : Iterable of float
arrays of energies at which cross sections are tabulated at
Returns
-------
openmc.data.PhotonReaction
Photon reaction data
"""
# Create instance
rx = cls(mt)
# Cross sections
xs = group['xs'][()]
# Replace zero elements to small non-zero to enable log-log
xs[xs == 0.0] = np.exp(-500.0)
# Threshold
threshold_idx = 0
if 'threshold_idx' in group['xs'].attrs:
threshold_idx = group['xs'].attrs['threshold_idx']
# Store cross section
rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5])
# Check for anomalous scattering factor
if 'anomalous_real' in group:
rx.anomalous_real = Tabulated1D.from_hdf5(group['anomalous_real'])
if 'anomalous_imag' in group:
rx.anomalous_imag = Tabulated1D.from_hdf5(group['anomalous_imag'])
# Check for factors / scattering functions
if 'scattering_factor' in group:
rx.scattering_factor = Tabulated1D.from_hdf5(group['scattering_factor'])
return rx
def to_hdf5(self, group, energy, Z):
"""Write photon reaction to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
energy : Iterable of float
arrays of energies at which cross sections are tabulated at
Z : int
atomic number
"""
# Write cross sections
if self.mt >= 534 and self.mt <= 572:
# Determine threshold
threshold = self.xs.x[0]
idx = np.searchsorted(energy, threshold, side='right') - 1
# Interpolate cross section onto union grid and write
photoionization = self.xs(energy[idx:])
group.create_dataset('xs', data=photoionization)
assert len(energy) == len(photoionization) + idx
group['xs'].attrs['threshold_idx'] = idx
else:
group.create_dataset('xs', data=self.xs(energy))
# Write scattering factor
if self.scattering_factor is not None:
if self.mt == 502:
# Create integrated form factor
ff = deepcopy(self.scattering_factor)
ff.x *= ff.x
ff.y *= ff.y/Z**2
int_ff = Tabulated1D(ff.x, ff.integral())
int_ff.to_hdf5(group, 'integrated_scattering_factor')
self.scattering_factor.to_hdf5(group, 'scattering_factor')
if self.anomalous_real is not None:
self.anomalous_real.to_hdf5(group, 'anomalous_real')
if self.anomalous_imag is not None:
self.anomalous_imag.to_hdf5(group, 'anomalous_imag')

10
openmc/data/reaction.py
View file

@ -51,7 +51,9 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
189: '(n,nta)', 190: '(n,2n2p)', 191: '(n,p3He)',
192: '(n,d3He)', 193: '(n,3Hea)', 194: '(n,4n2p)',
195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)',
198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)',
198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 203: '(n,Xp)',
204: '(n,Xd)', 205: '(n,Xt)', 206: '(n,X3He)', 207: '(n,Xa)',
301: 'heating', 444: 'damage-energy',
649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
849: '(n,ac)', 891: '(n,2nc)'}
REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)})
@ -936,7 +938,7 @@ class Reaction(EqualityMixin):
'Could not create reaction cross section for MT={} '
'at T={} because no corresponding energy grid '
'exists.'.format(mt, T))
xs = Tgroup['xs'].value
xs = Tgroup['xs'][()]
threshold_idx = Tgroup['xs'].attrs['threshold_idx'] - 1
tabulated_xs = Tabulated1D(energy[T][threshold_idx:], xs)
tabulated_xs._threshold_idx = threshold_idx
@ -988,6 +990,10 @@ class Reaction(EqualityMixin):
# Read reaction cross section
xs = ace.xss[ace.jxs[7] + loc + 1:ace.jxs[7] + loc + 1 + n_energy]
# For damage energy production, convert to eV
if mt == 444:
xs *= EV_PER_MEV
# Fix negatives -- known issue for Y89 in JEFF 3.2
if np.any(xs < 0.0):
warn("Negative cross sections found for MT={} in {}. Setting "

BIN
openmc/data/stopping_powers.h5
View file

Binary file not shown.

24
openmc/data/thermal.py
View file

@ -200,8 +200,8 @@ class CoherentElastic(EqualityMixin):
Coherent elastic scattering cross section
"""
bragg_edges = dataset.value[0, :]
factors = dataset.value[1, :]
bragg_edges = dataset[0, :]
factors = dataset[1, :]
return cls(bragg_edges, factors)
@ -414,7 +414,7 @@ class ThermalScattering(EqualityMixin):
kTg = group['kTs']
kTs = []
for temp in kTg:
kTs.append(kTg[temp].value)
kTs.append(kTg[temp][()])
temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" for kT in kTs]
table = cls(name, atomic_weight_ratio, kTs)
@ -438,7 +438,7 @@ class ThermalScattering(EqualityMixin):
# Angular distribution
if 'mu_out' in elastic_group:
table.elastic_mu_out[T] = elastic_group['mu_out'].value
table.elastic_mu_out[T] = elastic_group['mu_out'][()]
# Read thermal inelastic scattering
if 'inelastic' in Tgroup:
@ -446,8 +446,8 @@ class ThermalScattering(EqualityMixin):
table.inelastic_xs[T] = Tabulated1D.from_hdf5(
inelastic_group['xs'])
if table.secondary_mode in ('equal', 'skewed'):
table.inelastic_e_out[T] = inelastic_group['energy_out'].value
table.inelastic_mu_out[T] = inelastic_group['mu_out'].value
table.inelastic_e_out[T] = inelastic_group['energy_out'][()]
table.inelastic_mu_out[T] = inelastic_group['mu_out'][()]
elif table.secondary_mode == 'continuous':
table.inelastic_dist[T] = AngleEnergy.from_hdf5(
inelastic_group)
@ -610,7 +610,8 @@ class ThermalScattering(EqualityMixin):
return table
@classmethod
def from_njoy(cls, filename, filename_thermal, temperatures=None, **kwargs):
def from_njoy(cls, filename, filename_thermal, temperatures=None,
evaluation=None, evaluation_thermal=None, **kwargs):
"""Generate incident neutron data by running NJOY.
Parameters
@ -623,6 +624,13 @@ class ThermalScattering(EqualityMixin):
Temperatures in Kelvin to produce data at. If omitted, data is
produced at all temperatures in the ENDF thermal scattering
sublibrary.
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF neutron sublibrary file contains multiple material
evaluations, this argument indicates which evaluation to use.
evaluation_thermal : openmc.data.endf.Evaluation, optional
If the ENDF thermal scattering sublibrary file contains multiple
material evaluations, this argument indicates which evaluation to
use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace_thermal`
@ -636,6 +644,8 @@ class ThermalScattering(EqualityMixin):
# Run NJOY to create an ACE library
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
kwargs['evaluation'] = evaluation
kwargs['evaluation_thermal'] = evaluation_thermal
make_ace_thermal(filename, filename_thermal, temperatures,
ace_file, xsdir_file, **kwargs)

4
openmc/data/urr.py
View file

@ -166,8 +166,8 @@ class ProbabilityTables(EqualityMixin):
absorption_flag = group.attrs['absorption']
multiply_smooth = bool(group.attrs['multiply_smooth'])
energy = group['energy'].value
table = group['table'].value
energy = group['energy'][()]
table = group['table'][()]
return cls(energy, table, interpolation, inelastic_flag,
absorption_flag, multiply_smooth)

2
openmc/deplete/results_list.py
View file

@ -20,7 +20,7 @@ class ResultsList(list):
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
# Get number of results stored
n = fh["number"].value.shape[0]
n = fh["number"][...].shape[0]
for i in range(n):
self.append(Results.from_hdf5(fh, i))

24
openmc/examples.py
View file

@ -51,8 +51,8 @@ def pwr_pin_cell():
# Instantiate ZCylinder surfaces
pitch = 1.26
fuel_or = openmc.ZCylinder(x0=0, y0=0, R=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(x0=0, y0=0, R=0.45720, name='Clad OR')
fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
left = openmc.XPlane(x0=-pitch/2, name='left', boundary_type='reflective')
right = openmc.XPlane(x0=pitch/2, name='right', boundary_type='reflective')
bottom = openmc.YPlane(y0=-pitch/2, name='bottom',
@ -256,14 +256,14 @@ def pwr_core():
bot_nozzle, top_nozzle, top_fa, bot_fa)
# Define surfaces.
s1 = openmc.ZCylinder(R=0.41, surface_id=1)
s2 = openmc.ZCylinder(R=0.475, surface_id=2)
s3 = openmc.ZCylinder(R=0.56, surface_id=3)
s4 = openmc.ZCylinder(R=0.62, surface_id=4)
s5 = openmc.ZCylinder(R=187.6, surface_id=5)
s6 = openmc.ZCylinder(R=209.0, surface_id=6)
s7 = openmc.ZCylinder(R=229.0, surface_id=7)
s8 = openmc.ZCylinder(R=249.0, surface_id=8, boundary_type='vacuum')
s1 = openmc.ZCylinder(r=0.41, surface_id=1)
s2 = openmc.ZCylinder(r=0.475, surface_id=2)
s3 = openmc.ZCylinder(r=0.56, surface_id=3)
s4 = openmc.ZCylinder(r=0.62, surface_id=4)
s5 = openmc.ZCylinder(r=187.6, surface_id=5)
s6 = openmc.ZCylinder(r=209.0, surface_id=6)
s7 = openmc.ZCylinder(r=229.0, surface_id=7)
s8 = openmc.ZCylinder(r=249.0, surface_id=8, boundary_type='vacuum')
s31 = openmc.ZPlane(z0=-229.0, surface_id=31, boundary_type='vacuum')
s32 = openmc.ZPlane(z0=-199.0, surface_id=32)
@ -473,8 +473,8 @@ def pwr_assembly():
model.materials = (fuel, clad, hot_water)
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(x0=0, y0=0, R=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(x0=0, y0=0, R=0.45720, name='Clad OR')
fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
# Create boundary planes to surround the geometry
pitch = 21.42

32
openmc/filter.py
View file

@ -170,19 +170,19 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
# If the HDF5 'type' variable matches this class's short_name, then
# there is no overriden from_hdf5 method. Pass the bins to __init__.
if group['type'].value.decode() == cls.short_name.lower():
out = cls(group['bins'].value, filter_id=filter_id)
out._num_bins = group['n_bins'].value
if group['type'][()].decode() == cls.short_name.lower():
out = cls(group['bins'][()], filter_id=filter_id)
out._num_bins = group['n_bins'][()]
return out
# Search through all subclasses and find the one matching the HDF5
# 'type'. Call that class's from_hdf5 method.
for subclass in cls._recursive_subclasses():
if group['type'].value.decode() == subclass.short_name.lower():
if group['type'][()].decode() == subclass.short_name.lower():
return subclass.from_hdf5(group, **kwargs)
raise ValueError("Unrecognized Filter class: '"
+ group['type'].value.decode() + "'")
+ group['type'][()].decode() + "'")
@property
def bins(self):
@ -618,16 +618,16 @@ class MeshFilter(Filter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
if 'meshes' not in kwargs:
raise ValueError(cls.__name__ + " requires a 'meshes' keyword "
"argument.")
mesh_id = group['bins'].value
mesh_id = group['bins'][()]
mesh_obj = kwargs['meshes'][mesh_id]
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
@ -1191,15 +1191,15 @@ class DistribcellFilter(Filter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(group['bins'].value, filter_id=filter_id)
out._num_bins = group['n_bins'].value
out = cls(group['bins'][()], filter_id=filter_id)
out._num_bins = group['n_bins'][()]
return out
@ -1638,13 +1638,13 @@ class EnergyFunctionFilter(Filter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
energy = group['energy'].value
y = group['y'].value
energy = group['energy'][()]
y = group['y'][()]
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
return cls(energy, y, filter_id=filter_id)

32
openmc/filter_expansion.py
View file

@ -92,14 +92,14 @@ class LegendreFilter(ExpansionFilter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(group['order'].value, filter_id)
out = cls(group['order'][()], filter_id)
return out
@ -198,15 +198,15 @@ class SpatialLegendreFilter(ExpansionFilter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
order = group['order'].value
axis = group['axis'].value.decode()
min_, max_ = group['min'].value, group['max'].value
order = group['order'][()]
axis = group['axis'][()].decode()
min_, max_ = group['min'][()], group['max'][()]
return cls(order, axis, min_, max_, filter_id)
@ -294,15 +294,15 @@ class SphericalHarmonicsFilter(ExpansionFilter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(group['order'].value, filter_id)
out.cosine = group['cosine'].value.decode()
out = cls(group['order'][()], filter_id)
out.cosine = group['cosine'][()].decode()
return out
@ -437,14 +437,14 @@ class ZernikeFilter(ExpansionFilter):
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
+ group['type'][()].decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
order = group['order'].value
x, y, r = group['x'].value, group['y'].value, group['r'].value
order = group['order'][()]
x, y, r = group['x'][()], group['y'][()], group['r'][()]
return cls(order, x, y, r, filter_id)

7
openmc/geometry.py
View file

@ -97,9 +97,14 @@ class Geometry(object):
# Clean the indentation in the file to be user-readable
xml.clean_indentation(root_element)
# Check if path is a directory
p = Path(path)
if p.is_dir():
p /= 'geometry.xml'
# Write the XML Tree to the geometry.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True, encoding='utf-8')
tree.write(str(p), xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml(cls, path='geometry.xml', materials=None):

18
openmc/lattice.py
View file

@ -100,14 +100,14 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
"""
lattice_id = int(group.name.split('/')[-1].lstrip('lattice '))
name = group['name'].value.decode() if 'name' in group else ''
lattice_type = group['type'].value.decode()
name = group['name'][()].decode() if 'name' in group else ''
lattice_type = group['type'][()].decode()
if lattice_type == 'rectangular':
dimension = group['dimension'][...]
lower_left = group['lower_left'][...]
pitch = group['pitch'][...]
outer = group['outer'].value
outer = group['outer'][()]
universe_ids = group['universes'][...]
# Create the Lattice
@ -136,13 +136,13 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
lattice.universes = uarray
elif lattice_type == 'hexagonal':
n_rings = group['n_rings'].value
n_axial = group['n_axial'].value
center = group['center'][...]
pitch = group['pitch'][...]
outer = group['outer'].value
n_rings = group['n_rings'][()]
n_axial = group['n_axial'][()]
center = group['center'][()]
pitch = group['pitch'][()]
outer = group['outer'][()]
universe_ids = group['universes'][...]
universe_ids = group['universes'][()]
# Create the Lattice
lattice = openmc.HexLattice(lattice_id, name)

20
openmc/material.py
View file

@ -1,6 +1,7 @@
from collections import OrderedDict
from copy import deepcopy
from numbers import Real, Integral
from pathlib import Path
import warnings
from xml.etree import ElementTree as ET
@ -276,10 +277,10 @@ class Material(IDManagerMixin):
"""
mat_id = int(group.name.split('/')[-1].lstrip('material '))
name = group['name'].value.decode() if 'name' in group else ''
density = group['atom_density'].value
name = group['name'][()].decode() if 'name' in group else ''
density = group['atom_density'][()]
if 'nuclide_densities' in group:
nuc_densities = group['nuclide_densities'][...]
nuc_densities = group['nuclide_densities'][()]
# Create the Material
material = cls(mat_id, name)
@ -289,7 +290,7 @@ class Material(IDManagerMixin):
# Read the names of the S(a,b) tables for this Material and add them
if 'sab_names' in group:
sab_tables = group['sab_names'].value
sab_tables = group['sab_names'][()]
for sab_table in sab_tables:
name = sab_table.decode()
material.add_s_alpha_beta(name)
@ -298,13 +299,13 @@ class Material(IDManagerMixin):
material.set_density(density=density, units='atom/b-cm')
if 'nuclides' in group:
nuclides = group['nuclides'].value
nuclides = group['nuclides'][()]
# Add all nuclides to the Material
for fullname, density in zip(nuclides, nuc_densities):
name = fullname.decode().strip()
material.add_nuclide(name, percent=density, percent_type='ao')
if 'macroscopics' in group:
macroscopics = group['macroscopics'].value
macroscopics = group['macroscopics'][()]
# Add all macroscopics to the Material
for fullname in macroscopics:
name = fullname.decode().strip()
@ -1065,9 +1066,14 @@ class Materials(cv.CheckedList):
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)
# Check if path is a directory
p = Path(path)
if p.is_dir():
p /= 'materials.xml'
# Write the XML Tree to the materials.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True, encoding='utf-8')
tree.write(str(p), xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml(cls, path='materials.xml'):

10
openmc/mesh.py
View file

@ -174,11 +174,11 @@ class Mesh(IDManagerMixin):
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.type = group['type'].value.decode()
mesh.dimension = group['dimension'].value
mesh.lower_left = group['lower_left'].value
mesh.upper_right = group['upper_right'].value
mesh.width = group['width'].value
mesh.type = group['type'][()].decode()
mesh.dimension = group['dimension'][()]
mesh.lower_left = group['lower_left'][()]
mesh.upper_right = group['upper_right'][()]
mesh.width = group['width'][()]
return mesh

152
openmc/mgxs/__init__.py
View file

@ -6,26 +6,50 @@ from openmc.mgxs.mgxs import *
from openmc.mgxs.mdgxs import *
GROUP_STRUCTURES = {}
"""Dictionary of commonly used energy group structures, including "CASMO-X" (where X
is 2, 4, 8, 16, 25, 40 or 70) from the CASMO_ lattice physics code and other commonly
used activation_ energy group structures "VITAMIN-J-175", "TRIPOLI-315", "CCFE-709_"
and "UKAEA-1102_"
"""Dictionary of commonly used energy group structures:
- "CASMO-X" (where X is 2, 4, 8, 16, 25, 40 or 70) from the CASMO_ lattice
physics code
- "XMAS-172_" designed for LWR analysis ([SAR1990]_, [SAN2004]_)
- "SHEM-361_" designed for LWR analysis to eliminate self-shielding calculations
of thermal resonances ([HFA2005]_, [SAN2007]_, [HEB2008]_)
- activation_ energy group structures "VITAMIN-J-175", "TRIPOLI-315",
"CCFE-709_" and "UKAEA-1102_"
.. _CASMO: https://www.studsvik.com/SharepointFiles/CASMO-5%20Development%20and%20Applications.pdf
.. _XMAS-172: https://www-nds.iaea.org/wimsd/energy.htm
.. _SHEM-361: https://www.polymtl.ca/merlin/libraries.htm
.. _activation: https://fispact.ukaea.uk/wiki/Keyword:GETXS
.. _CCFE-709: https://fispact.ukaea.uk/wiki/CCFE-709_group_structure
.. _UKAEA-1102: https://fispact.ukaea.uk/wiki/UKAEA-1102_group_structure
.. [SAR1990] Sartori, E., OECD/NEA Data Bank: Standard Energy Group Structures
of Cross Section Libraries for Reactor Shielding, Reactor Cell and Fusion
Neutronics Applications: VITAMIN-J, ECCO-33, ECCO-2000 and XMAS JEF/DOC-315
Revision 3 - DRAFT (December 11, 1990).
.. [SAN2004] Santamarina, A., Collignon, C., & Garat, C. (2004). French
calculation schemes for light water reactor analysis. United States:
American Nuclear Society - ANS.
.. [HFA2005] Hfaiedh, N. & Santamarina, A., "Determination of the Optimized
SHEM Mesh for Neutron Transport Calculations," Proc. Top. Mtg. in
Mathematics & Computations, Supercomputing, Reactor Physics and Nuclear and
Biological Applications, September 12-15, Avignon, France, 2005.
.. [SAN2007] Santamarina, A. & Hfaiedh, N. (2007). The SHEM energy mesh for
accurate fuel depletion and BUC calculations. Proceedings of the International
Conference on Safety Criticality ICNC 2007, St Peterburg (Russia), Vol. I pp.
446-452.
.. [HEB2008] Hébert, Alain & Santamarina, Alain. (2008). Refinement of the
Santamarina-Hfaiedh energy mesh between 22.5 eV and 11.4 keV. International
Conference on the Physics of Reactors 2008, PHYSOR 08. 2. 929-938.
"""
GROUP_STRUCTURES['CASMO-2'] = np.array([
0., 6.25e-1, 2.e7])
0., 6.25e-1, 2.e7])
GROUP_STRUCTURES['CASMO-4'] = np.array([
0., 6.25e-1, 5.53e3, 8.21e5, 2.e7])
0., 6.25e-1, 5.53e3, 8.21e5, 2.e7])
GROUP_STRUCTURES['CASMO-8'] = np.array([
0., 5.8e-2, 1.4e-1, 2.8e-1, 6.25e-1, 4., 5.53e3, 8.21e5, 2.e7])
0., 5.8e-2, 1.4e-1, 2.8e-1, 6.25e-1, 4., 5.53e3, 8.21e5, 2.e7])
GROUP_STRUCTURES['CASMO-16'] = np.array([
0., 3.e-2, 5.8e-2, 1.4e-1, 2.8e-1, 3.5e-1, 6.25e-1, 8.5e-1,
9.72e-1, 1.02, 1.097, 1.15, 1.3, 4., 5.53e3, 8.21e5, 2.e7])
0., 3.e-2, 5.8e-2, 1.4e-1, 2.8e-1, 3.5e-1, 6.25e-1, 8.5e-1,
9.72e-1, 1.02, 1.097, 1.15, 1.3, 4., 5.53e3, 8.21e5, 2.e7])
GROUP_STRUCTURES['CASMO-25'] = np.array([
0., 3.e-2, 5.8e-2, 1.4e-1, 2.8e-1, 3.5e-1, 6.25e-1, 9.72e-1, 1.02, 1.097,
1.15, 1.855, 4., 9.877, 1.5968e1, 1.4873e2, 5.53e3, 9.118e3, 1.11e5, 5.e5,
@ -46,6 +70,42 @@ GROUP_STRUCTURES['CASMO-70'] = np.array([
3.6726e2, 9.069e2, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3,
9.118e3, 1.503e4, 2.478e4, 4.085e4, 6.734e4, 1.11e5, 1.83e5,
3.025e5, 5.e5, 8.21e5, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2.e7])
GROUP_STRUCTURES['XMAS-172'] = np.array([
1.00001e-05, 3.00000e-03, 5.00000e-03, 6.90000e-03, 1.00000e-02,
1.50000e-02, 2.00000e-02, 2.50000e-02, 3.00000e-02, 3.50000e-02,
4.20000e-02, 5.00000e-02, 5.80000e-02, 6.70000e-02, 7.70000e-02,
8.00000e-02, 9.50000e-02, 1.00001e-01, 1.15000e-01, 1.34000e-01,
1.40000e-01, 1.60000e-01, 1.80000e-01, 1.89000e-01, 2.20000e-01,
2.48000e-01, 2.80000e-01, 3.00000e-01, 3.14500e-01, 3.20000e-01,
3.50000e-01, 3.91000e-01, 4.00000e-01, 4.33000e-01, 4.85000e-01,
5.00000e-01, 5.40000e-01, 6.25000e-01, 7.05000e-01, 7.80000e-01,
7.90000e-01, 8.50000e-01, 8.60000e-01, 9.10000e-01, 9.30000e-01,
9.50000e-01, 9.72000e-01, 9.86000e-01, 9.96000e-01, 1.02000e+00,
1.03500e+00, 1.04500e+00, 1.07100e+00, 1.09700e+00, 1.11000e+00,
1.12535e+00, 1.15000e+00, 1.17000e+00, 1.23500e+00, 1.30000e+00,
1.33750e+00, 1.37000e+00, 1.44498e+00, 1.47500e+00, 1.50000e+00,
1.59000e+00, 1.67000e+00, 1.75500e+00, 1.84000e+00, 1.93000e+00,
2.02000e+00, 2.10000e+00, 2.13000e+00, 2.36000e+00, 2.55000e+00,
2.60000e+00, 2.72000e+00, 2.76792e+00, 3.30000e+00, 3.38075e+00,
4.00000e+00, 4.12925e+00, 5.04348e+00, 5.34643e+00, 6.16012e+00,
7.52398e+00, 8.31529e+00, 9.18981e+00, 9.90555e+00, 1.12245e+01,
1.37096e+01, 1.59283e+01, 1.94548e+01, 2.26033e+01, 2.49805e+01,
2.76077e+01, 3.05113e+01, 3.37201e+01, 3.72665e+01, 4.01690e+01,
4.55174e+01, 4.82516e+01, 5.15780e+01, 5.55951e+01, 6.79041e+01,
7.56736e+01, 9.16609e+01, 1.36742e+02, 1.48625e+02, 2.03995e+02,
3.04325e+02, 3.71703e+02, 4.53999e+02, 6.77287e+02, 7.48518e+02,
9.14242e+02, 1.01039e+03, 1.23410e+03, 1.43382e+03, 1.50733e+03,
2.03468e+03, 2.24867e+03, 3.35463e+03, 3.52662e+03, 5.00451e+03,
5.53084e+03, 7.46586e+03, 9.11882e+03, 1.11378e+04, 1.50344e+04,
1.66156e+04, 2.47875e+04, 2.73944e+04, 2.92830e+04, 3.69786e+04,
4.08677e+04, 5.51656e+04, 6.73795e+04, 8.22975e+04, 1.11090e+05,
1.22773e+05, 1.83156e+05, 2.47235e+05, 2.73237e+05, 3.01974e+05,
4.07622e+05, 4.50492e+05, 4.97871e+05, 5.50232e+05, 6.08101e+05,
8.20850e+05, 9.07180e+05, 1.00259e+06, 1.10803e+06, 1.22456e+06,
1.35335e+06, 1.65299e+06, 2.01897e+06, 2.23130e+06, 2.46597e+06,
3.01194e+06, 3.67879e+06, 4.49329e+06, 5.48812e+06, 6.06531e+06,
6.70320e+06, 8.18731e+06, 1.00000e+07, 1.16183e+07, 1.38403e+07,
1.49182e+07, 1.73325e+07, 1.96403e+07])
GROUP_STRUCTURES['VITAMIN-J-175'] = np.array([
1.0000e-5, 1.0000e-1, 4.1399e-1, 5.3158e-1, 6.8256e-1,
8.7643e-1, 1.1253, 1.4450, 1.8554, 2.3824, 3.0590,
@ -123,6 +183,80 @@ GROUP_STRUCTURES['TRIPOLI-315,'] = np.array([
7.788e6, 8.187e6, 8.607e6, 9.048e6, 9.512e6, 1.000e7, 1.051e7,
1.105e7, 1.162e7, 1.221e7, 1.284e7, 1.350e7, 1.384e7, 1.419e7,
1.455e7, 1.492e7, 1.568e7, 1.649e7, 1.691e7, 1.733e7, 1.964e7])
GROUP_STRUCTURES['SHEM-361'] = np.array([
0.00000e+00, 2.49990e-03, 4.55602e-03, 7.14526e-03, 1.04505e-02,
1.48300e-02, 2.00104e-02, 2.49394e-02, 2.92989e-02, 3.43998e-02,
4.02999e-02, 4.73019e-02, 5.54982e-02, 6.51999e-02, 7.64969e-02,
8.97968e-02, 1.04298e-01, 1.19995e-01, 1.37999e-01, 1.61895e-01,
1.90005e-01, 2.09610e-01, 2.31192e-01, 2.54997e-01, 2.79989e-01,
3.05012e-01, 3.25008e-01, 3.52994e-01, 3.90001e-01, 4.31579e-01,
4.75017e-01, 5.20011e-01, 5.54990e-01, 5.94993e-01, 6.24999e-01,
7.19999e-01, 8.00371e-01, 8.80024e-01, 9.19978e-01, 9.44022e-01,
9.63960e-01, 9.81959e-01, 9.96501e-01, 1.00904e+00, 1.02101e+00,
1.03499e+00, 1.07799e+00, 1.09198e+00, 1.10395e+00, 1.11605e+00,
1.12997e+00, 1.14797e+00, 1.16999e+00, 1.21397e+00, 1.25094e+00,
1.29304e+00, 1.33095e+00, 1.38098e+00, 1.41001e+00, 1.44397e+00,
1.51998e+00, 1.58803e+00, 1.66895e+00, 1.77997e+00, 1.90008e+00,
1.98992e+00, 2.07010e+00, 2.15695e+00, 2.21709e+00, 2.27299e+00,
2.33006e+00, 2.46994e+00, 2.55000e+00, 2.59009e+00, 2.62005e+00,
2.64004e+00, 2.70012e+00, 2.71990e+00, 2.74092e+00, 2.77512e+00,
2.88405e+00, 3.14211e+00, 3.54307e+00, 3.71209e+00, 3.88217e+00,
4.00000e+00, 4.21983e+00, 4.30981e+00, 4.41980e+00, 4.76785e+00,
4.93323e+00, 5.10997e+00, 5.21008e+00, 5.32011e+00, 5.38003e+00,
5.41025e+00, 5.48817e+00, 5.53004e+00, 5.61979e+00, 5.72015e+00,
5.80021e+00, 5.96014e+00, 6.05991e+00, 6.16011e+00, 6.28016e+00,
6.35978e+00, 6.43206e+00, 6.48178e+00, 6.51492e+00, 6.53907e+00,
6.55609e+00, 6.57184e+00, 6.58829e+00, 6.60611e+00, 6.63126e+00,
6.71668e+00, 6.74225e+00, 6.75981e+00, 6.77605e+00, 6.79165e+00,
6.81070e+00, 6.83526e+00, 6.87021e+00, 6.91778e+00, 6.99429e+00,
7.13987e+00, 7.38015e+00, 7.60035e+00, 7.73994e+00, 7.83965e+00,
7.97008e+00, 8.13027e+00, 8.30032e+00, 8.52407e+00, 8.67369e+00,
8.80038e+00, 8.97995e+00, 9.14031e+00, 9.50002e+00, 1.05793e+01,
1.08038e+01, 1.10529e+01, 1.12694e+01, 1.15894e+01, 1.17094e+01,
1.18153e+01, 1.19795e+01, 1.21302e+01, 1.23086e+01, 1.24721e+01,
1.26000e+01, 1.33297e+01, 1.35460e+01, 1.40496e+01, 1.42505e+01,
1.44702e+01, 1.45952e+01, 1.47301e+01, 1.48662e+01, 1.57792e+01,
1.60498e+01, 1.65501e+01, 1.68305e+01, 1.74457e+01, 1.75648e+01,
1.77590e+01, 1.79591e+01, 1.90848e+01, 1.91997e+01, 1.93927e+01,
1.95974e+01, 2.00734e+01, 2.02751e+01, 2.04175e+01, 2.05199e+01,
2.06021e+01, 2.06847e+01, 2.07676e+01, 2.09763e+01, 2.10604e+01,
2.11448e+01, 2.12296e+01, 2.13360e+01, 2.14859e+01, 2.17018e+01,
2.20011e+01, 2.21557e+01, 2.23788e+01, 2.25356e+01, 2.46578e+01,
2.78852e+01, 3.16930e+01, 3.30855e+01, 3.45392e+01, 3.56980e+01,
3.60568e+01, 3.64191e+01, 3.68588e+01, 3.73038e+01, 3.77919e+01,
3.87874e+01, 3.97295e+01, 4.12270e+01, 4.21441e+01, 4.31246e+01,
4.41721e+01, 4.52904e+01, 4.62053e+01, 4.75173e+01, 4.92591e+01,
5.17847e+01, 5.29895e+01, 5.40600e+01, 5.70595e+01, 5.99250e+01,
6.23083e+01, 6.36306e+01, 6.45923e+01, 6.50460e+01, 6.55029e+01,
6.58312e+01, 6.61612e+01, 6.64929e+01, 6.68261e+01, 6.90682e+01,
7.18869e+01, 7.35595e+01, 7.63322e+01, 7.93679e+01, 8.39393e+01,
8.87741e+01, 9.33256e+01, 9.73287e+01, 1.00594e+02, 1.01098e+02,
1.01605e+02, 1.02115e+02, 1.03038e+02, 1.05646e+02, 1.10288e+02,
1.12854e+02, 1.15480e+02, 1.16524e+02, 1.17577e+02, 1.20554e+02,
1.26229e+02, 1.32701e+02, 1.39504e+02, 1.46657e+02, 1.54176e+02,
1.63056e+02, 1.67519e+02, 1.75229e+02, 1.83295e+02, 1.84952e+02,
1.86251e+02, 1.87559e+02, 1.88877e+02, 1.90204e+02, 1.93078e+02,
1.95996e+02, 2.00958e+02, 2.12108e+02, 2.24325e+02, 2.35590e+02,
2.41796e+02, 2.56748e+02, 2.68297e+02, 2.76468e+02, 2.84888e+02,
2.88327e+02, 2.95922e+02, 3.19928e+02, 3.35323e+02, 3.53575e+02,
3.71703e+02, 3.90760e+02, 4.19094e+02, 4.53999e+02, 5.01746e+02,
5.39204e+02, 5.77146e+02, 5.92941e+02, 6.00099e+02, 6.12834e+02,
6.46837e+02, 6.77287e+02, 7.48517e+02, 8.32218e+02, 9.09681e+02,
9.82494e+02, 1.06432e+03, 1.13467e+03, 1.34358e+03, 1.58620e+03,
1.81183e+03, 2.08410e+03, 2.39729e+03, 2.70024e+03, 2.99618e+03,
3.48107e+03, 4.09735e+03, 5.00451e+03, 6.11252e+03, 7.46585e+03,
9.11881e+03, 1.11377e+04, 1.36037e+04, 1.48997e+04, 1.62005e+04,
1.85847e+04, 2.26994e+04, 2.49991e+04, 2.61001e+04, 2.73944e+04,
2.92810e+04, 3.34596e+04, 3.69786e+04, 4.08677e+04, 4.99159e+04,
5.51656e+04, 6.73794e+04, 8.22974e+04, 9.46645e+04, 1.15624e+05,
1.22773e+05, 1.40000e+05, 1.64999e+05, 1.95008e+05, 2.30014e+05,
2.67826e+05, 3.20646e+05, 3.83884e+05, 4.12501e+05, 4.56021e+05,
4.94002e+05, 5.78443e+05, 7.06511e+05, 8.60006e+05, 9.51119e+05,
1.05115e+06, 1.16205e+06, 1.28696e+06, 1.33694e+06, 1.40577e+06,
1.63654e+06, 1.90139e+06, 2.23130e+06, 2.72531e+06, 3.32871e+06,
4.06569e+06, 4.96585e+06, 6.06530e+06, 6.70319e+06, 7.40817e+06,
8.18730e+06, 9.04836e+06, 9.99999e+06, 1.16183e+07, 1.38403e+07,
1.49182e+07, 1.96403e+07])
GROUP_STRUCTURES['CCFE-709'] = np.array([
1.e-5, 1.0471e-5, 1.0965e-5, 1.1482e-5, 1.2023e-5,
1.2589e-5, 1.3183e-5, 1.3804e-5, 1.4454e-5, 1.5136e-5,

8
openmc/mgxs_library.py
View file

@ -2183,7 +2183,7 @@ class XSdata(object):
kTs_group = group['kTs']
float_temperatures = []
for temperature in temperatures:
kT = kTs_group[temperature].value
kT = kTs_group[temperature][()]
float_temperatures.append(kT / openmc.data.K_BOLTZMANN)
attrs = group.attrs.keys()
@ -2219,7 +2219,7 @@ class XSdata(object):
for xs_type in xs_types:
set_func = 'set_' + xs_type.replace(' ', '_').replace('-', '_')
if xs_type in temperature_group:
getattr(data, set_func)(temperature_group[xs_type].value,
getattr(data, set_func)(temperature_group[xs_type][()],
float_temp)
scatt_group = temperature_group['scatter_data']
@ -2227,7 +2227,7 @@ class XSdata(object):
# Get scatter matrix and 'un-flatten' it
g_max = scatt_group['g_max']
g_min = scatt_group['g_min']
flat_scatter = scatt_group['scatter_matrix'].value
flat_scatter = scatt_group['scatter_matrix'][()]
scatter_matrix = np.zeros(data.xs_shapes["[G][G'][Order]"])
G = data.energy_groups.num_groups
if data.representation == 'isotropic':
@ -2259,7 +2259,7 @@ class XSdata(object):
# Repeat for multiplicity
if 'multiplicity_matrix' in scatt_group:
flat_mult = scatt_group['multiplicity_matrix'].value
flat_mult = scatt_group['multiplicity_matrix'][()]
mult_matrix = np.zeros(data.xs_shapes["[G][G']"])
flat_index = 0
for p in range(Np):

71
openmc/model/funcs.py
View file

@ -2,8 +2,9 @@ from collections import OrderedDict
from collections.abc import Iterable
from math import sqrt
from numbers import Real
from functools import partial
from openmc import XPlane, YPlane, Plane, ZCylinder
from openmc import XPlane, YPlane, Plane, ZCylinder, Quadric
from openmc.checkvalue import check_type, check_value
import openmc.data
@ -237,16 +238,16 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
c = sqrt(3.)/3.
# y = -x/sqrt(3) + a
upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type)
upper_right = Plane(a=c, b=1., d=l+x*c+y, boundary_type=boundary_type)
# y = x/sqrt(3) + a
upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type)
upper_left = Plane(a=-c, b=1., d=l-x*c+y, boundary_type=boundary_type)
# y = x/sqrt(3) - a
lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type)
lower_right = Plane(a=-c, b=1., d=-l-x*c+y, boundary_type=boundary_type)
# y = -x/sqrt(3) - a
lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type)
lower_left = Plane(a=c, b=1., d=-l+x*c+y, boundary_type=boundary_type)
prism = -right & +left & -upper_right & -upper_left & \
+lower_right & +lower_left
@ -262,17 +263,17 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
c = sqrt(3.)
# y = -sqrt(3)*(x - a)
upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type)
upper_right = Plane(a=c, b=1., d=c*l+x*c+y, boundary_type=boundary_type)
# y = sqrt(3)*(x + a)
lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y,
lower_right = Plane(a=-c, b=1., d=-c*l-x*c+y,
boundary_type=boundary_type)
# y = -sqrt(3)*(x + a)
lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type)
lower_left = Plane(a=c, b=1., d=-c*l+x*c+y, boundary_type=boundary_type)
# y = sqrt(3)*(x + a)
upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type)
upper_left = Plane(a=-c, b=1., d=c*l-x*c+y, boundary_type=boundary_type)
prism = -top & +bottom & -upper_right & +lower_right & \
+lower_left & -upper_left
@ -292,8 +293,8 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
t = l - corner_radius/c
# Cylinder with corner radius and boundary type pre-applied
cyl1 = partial(ZCylinder, R=corner_radius, boundary_type=boundary_type)
cyl2 = partial(ZCylinder, R=corner_radius/(2*c),
cyl1 = partial(ZCylinder, r=corner_radius, boundary_type=boundary_type)
cyl2 = partial(ZCylinder, r=corner_radius/(2*c),
boundary_type=boundary_type)
if orientation == 'x':
@ -345,6 +346,54 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
return prism
def cylinder_from_points(p1, p2, r, **kwargs):
"""Return cylinder defined by two points passing through its center.
Parameters
----------
p1, p2 : 3-tuples
Coordinates of two points that pass through the center of the cylinder
r : float
Radius of the cylinder
kwargs : dict
Keyword arguments passed to the :class:`openmc.Quadric` constructor
Returns
-------
openmc.Quadric
Quadric surface representing the cylinder.
"""
# Get x, y, z coordinates of two points
x1, y1, z1 = p1
x2, y2, z2 = p2
# Define intermediate terms
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
cx = y1*z2 + y2*z1
cy = -(x1*z2 + x2*z1)
cz = x1*y2 + x2*y1
# Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the
# cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|.
# Expanding out all terms and grouping according to what Quadric expects
# gives the following coefficients.
kwargs['a'] = dy*dy + dz*dz
kwargs['b'] = dx*dx + dz*dz
kwargs['c'] = dx*dx + dy*dy
kwargs['d'] = -2*dx*dy
kwargs['e'] = -2*dy*dz
kwargs['f'] = -2*dx*dz
kwargs['g'] = cy*dz - cz*dy
kwargs['h'] = cz*dx - cx*dz
kwargs['j'] = cx*dy - cy*dx
kwargs['k'] = -(dx*dx + dy*dy + dz*dz)*r*r
return openmc.Quadric(**kwargs)
def subdivide(surfaces):
"""Create regions separated by a series of surfaces.

29
openmc/model/model.py
View file

@ -1,4 +1,5 @@
from collections.abc import Iterable
from pathlib import Path
import openmc
from openmc.checkvalue import check_type, check_value
@ -154,27 +155,39 @@ class Model(object):
'si_celi', 'si_leqi', 'celi', 'leqi'))
getattr(dep.integrator, method)(op, timesteps, **kwargs)
def export_to_xml(self):
"""Export model to XML files."""
def export_to_xml(self, directory='.'):
"""Export model to XML files.
self.settings.export_to_xml()
Parameters
----------
directory : str
Directory to write XML files to. If it doesn't exist already, it
will be created.
"""
# Create directory if
d = Path(directory)
if not d.is_dir():
d.mkdir(parents=True)
self.settings.export_to_xml(d)
if not self.settings.dagmc:
self.geometry.export_to_xml()
self.geometry.export_to_xml(d)
# If a materials collection was specified, export it. Otherwise, look
# for all materials in the geometry and use that to automatically build
# a collection.
if self.materials:
self.materials.export_to_xml()
self.materials.export_to_xml(d)
else:
materials = openmc.Materials(self.geometry.get_all_materials()
.values())
materials.export_to_xml()
materials.export_to_xml(d)
if self.tallies:
self.tallies.export_to_xml()
self.tallies.export_to_xml(d)
if self.plots:
self.plots.export_to_xml()
self.plots.export_to_xml(d)
def run(self, **kwargs):
"""Creates the XML files, runs OpenMC, and returns k-effective

2
openmc/model/triso.py
View file

@ -49,7 +49,7 @@ class TRISO(openmc.Cell):
"""
def __init__(self, outer_radius, fill, center=(0., 0., 0.)):
self._surface = openmc.Sphere(R=outer_radius)
self._surface = openmc.Sphere(r=outer_radius)
super().__init__(fill=fill, region=-self._surface)
self.center = np.asarray(center)

22
openmc/particle_restart.py
View file

@ -49,44 +49,44 @@ class Particle(object):
@property
def current_batch(self):
return self._f['current_batch'].value
return self._f['current_batch'][()]
@property
def current_generation(self):
return self._f['current_generation'].value
return self._f['current_generation'][()]
@property
def energy(self):
return self._f['energy'].value
return self._f['energy'][()]
@property
def generations_per_batch(self):
return self._f['generations_per_batch'].value
return self._f['generations_per_batch'][()]
@property
def id(self):
return self._f['id'].value
return self._f['id'][()]
@property
def type(self):
return self._f['type'].value
return self._f['type'][()]
@property
def n_particles(self):
return self._f['n_particles'].value
return self._f['n_particles'][()]
@property
def run_mode(self):
return self._f['run_mode'].value.decode()
return self._f['run_mode'][()].decode()
@property
def uvw(self):
return self._f['uvw'].value
return self._f['uvw'][()]
@property
def weight(self):
return self._f['weight'].value
return self._f['weight'][()]
@property
def xyz(self):
return self._f['xyz'].value
return self._f['xyz'][()]

10
openmc/plots.py
View file

@ -1,9 +1,10 @@
from collections.abc import Iterable, Mapping
from numbers import Real, Integral
from xml.etree import ElementTree as ET
from pathlib import Path
import subprocess
import sys
import warnings
from xml.etree import ElementTree as ET
import numpy as np
@ -818,6 +819,11 @@ class Plots(cv.CheckedList):
# Clean the indentation in the file to be user-readable
clean_indentation(self._plots_file)
# Check if path is a directory
p = Path(path)
if p.is_dir():
p /= 'plots.xml'
# Write the XML Tree to the plots.xml file
tree = ET.ElementTree(self._plots_file)
tree.write(path, xml_declaration=True, encoding='utf-8', method="xml")
tree.write(str(p), xml_declaration=True, encoding='utf-8')

93
openmc/region.py
View file

@ -229,14 +229,28 @@ class Region(metaclass=ABCMeta):
clone : openmc.Region
The clone of this region
Raises
------
NotImplementedError
This method is not implemented for the abstract region class.
"""
pass
@abstractmethod
def translate(self, vector, memo=None):
"""Translate region in given direction
Parameters
----------
vector : iterable of float
Direction in which region should be translated
memo : dict or None
Dictionary used for memoization. This parameter is used internally
and should not be specified by the user.
Returns
-------
openmc.Region
Translated region
"""
raise NotImplementedError('The clone method is not implemented for '
'the abstract region class.')
pass
class Intersection(Region, MutableSequence):
@ -247,7 +261,7 @@ class Intersection(Region, MutableSequence):
following example:
>>> equator = openmc.ZPlane(z0=0.0)
>>> earth = openmc.Sphere(R=637.1e6)
>>> earth = openmc.Sphere(r=637.1e6)
>>> northern_hemisphere = -earth & +equator
>>> southern_hemisphere = -earth & -equator
>>> type(northern_hemisphere)
@ -352,6 +366,27 @@ class Intersection(Region, MutableSequence):
clone[:] = [n.clone(memo) for n in self]
return clone
def translate(self, vector, memo=None):
"""Translate region in given direction
Parameters
----------
vector : iterable of float
Direction in which region should be translated
memo : dict or None
Dictionary used for memoization. This parameter is used internally
and should not be specified by the user.
Returns
-------
openmc.Intersection
Translated region
"""
if memo is None:
memo = {}
return type(self)(n.translate(vector, memo) for n in self)
class Union(Region, MutableSequence):
r"""Union of two or more regions.
@ -361,7 +396,7 @@ class Union(Region, MutableSequence):
example:
>>> s1 = openmc.ZPlane(z0=0.0)
>>> s2 = openmc.Sphere(R=637.1e6)
>>> s2 = openmc.Sphere(r=637.1e6)
>>> type(-s2 | +s1)
<class 'openmc.region.Union'>
@ -464,6 +499,27 @@ class Union(Region, MutableSequence):
clone[:] = [n.clone(memo) for n in self]
return clone
def translate(self, vector, memo=None):
"""Translate region in given direction
Parameters
----------
vector : iterable of float
Direction in which region should be translated
memo : dict or None
Dictionary used for memoization. This parameter is used internally
and should not be specified by the user.
Returns
-------
openmc.Union
Translated region
"""
if memo is None:
memo = {}
return type(self)(n.translate(vector, memo) for n in self)
class Complement(Region):
"""Complement of a region.
@ -584,3 +640,24 @@ class Complement(Region):
clone = deepcopy(self)
clone.node = self.node.clone(memo)
return clone
def translate(self, vector, memo=None):
"""Translate region in given direction
Parameters
----------
vector : iterable of float
Direction in which region should be translated
memo : dict or None
Dictionary used for memoization. This parameter is used internally
and should not be specified by the user.
Returns
-------
openmc.Complement
Translated region
"""
if memo is None:
memo = {}
return type(self)(self.node.translate(vector, memo))

27
openmc/settings.py
View file

@ -1,4 +1,5 @@
from collections.abc import Iterable, MutableSequence, Mapping
from pathlib import Path
from numbers import Real, Integral
import warnings
from xml.etree import ElementTree as ET
@ -129,8 +130,6 @@ class Settings(object):
range. 'multipole' is a boolean indicating whether or not the windowed
multipole method should be used to evaluate resolved resonance cross
sections.
threads : int
Number of OpenMP threads
trace : tuple or list
Show detailed information about a single particle, indicated by three
integers: the batch number, generation number, and particle number
@ -196,7 +195,6 @@ class Settings(object):
self._statepoint = {}
self._sourcepoint = {}
self._threads = None
self._no_reduce = None
self._verbosity = None
@ -311,10 +309,6 @@ class Settings(object):
def statepoint(self):
return self._statepoint
@property
def threads(self):
return self._threads
@property
def no_reduce(self):
return self._no_reduce
@ -622,12 +616,6 @@ class Settings(object):
self._temperature = temperature
@threads.setter
def threads(self, threads):
cv.check_type('number of threads', threads, Integral)
cv.check_greater_than('number of threads', threads, 0)
self._threads = threads
@trace.setter
def trace(self, trace):
cv.check_type('trace', trace, Iterable, Integral)
@ -884,11 +872,6 @@ class Settings(object):
else:
element.text = str(value)
def _create_threads_subelement(self, root):
if self._threads is not None:
element = ET.SubElement(root, "threads")
element.text = str(self._threads)
def _create_trace_subelement(self, root):
if self._trace is not None:
element = ET.SubElement(root, "trace")
@ -979,7 +962,6 @@ class Settings(object):
self._create_entropy_mesh_subelement(root_element)
self._create_trigger_subelement(root_element)
self._create_no_reduce_subelement(root_element)
self._create_threads_subelement(root_element)
self._create_verbosity_subelement(root_element)
self._create_tabular_legendre_subelements(root_element)
self._create_temperature_subelements(root_element)
@ -995,6 +977,11 @@ class Settings(object):
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)
# Check if path is a directory
p = Path(path)
if p.is_dir():
p /= 'settings.xml'
# Write the XML Tree to the settings.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True, encoding='utf-8', method="xml")
tree.write(str(p), xml_declaration=True, encoding='utf-8')

76
openmc/statepoint.py
View file

@ -161,35 +161,35 @@ class StatePoint(object):
@property
def cmfd_balance(self):
return self._f['cmfd/cmfd_balance'].value if self.cmfd_on else None
return self._f['cmfd/cmfd_balance'][()] if self.cmfd_on else None
@property
def cmfd_dominance(self):
return self._f['cmfd/cmfd_dominance'].value if self.cmfd_on else None
return self._f['cmfd/cmfd_dominance'][()] if self.cmfd_on else None
@property
def cmfd_entropy(self):
return self._f['cmfd/cmfd_entropy'].value if self.cmfd_on else None
return self._f['cmfd/cmfd_entropy'][()] if self.cmfd_on else None
@property
def cmfd_indices(self):
return self._f['cmfd/indices'].value if self.cmfd_on else None
return self._f['cmfd/indices'][()] if self.cmfd_on else None
@property
def cmfd_src(self):
if self.cmfd_on:
data = self._f['cmfd/cmfd_src'].value
data = self._f['cmfd/cmfd_src'][()]
return np.reshape(data, tuple(self.cmfd_indices), order='F')
else:
return None
@property
def cmfd_srccmp(self):
return self._f['cmfd/cmfd_srccmp'].value if self.cmfd_on else None
return self._f['cmfd/cmfd_srccmp'][()] if self.cmfd_on else None
@property
def current_batch(self):
return self._f['current_batch'].value
return self._f['current_batch'][()]
@property
def date_and_time(self):
@ -199,7 +199,7 @@ class StatePoint(object):
@property
def entropy(self):
if self.run_mode == 'eigenvalue':
return self._f['entropy'].value
return self._f['entropy'][()]
else:
return None
@ -220,14 +220,14 @@ class StatePoint(object):
@property
def generations_per_batch(self):
if self.run_mode == 'eigenvalue':
return self._f['generations_per_batch'].value
return self._f['generations_per_batch'][()]
else:
return None
@property
def global_tallies(self):
if self._global_tallies is None:
data = self._f['global_tallies'].value
data = self._f['global_tallies'][()]
gt = np.zeros(data.shape[0], dtype=[
('name', 'a14'), ('sum', 'f8'), ('sum_sq', 'f8'),
('mean', 'f8'), ('std_dev', 'f8')])
@ -248,42 +248,42 @@ class StatePoint(object):
@property
def k_cmfd(self):
if self.cmfd_on:
return self._f['cmfd/k_cmfd'].value
return self._f['cmfd/k_cmfd'][()]
else:
return None
@property
def k_generation(self):
if self.run_mode == 'eigenvalue':
return self._f['k_generation'].value
return self._f['k_generation'][()]
else:
return None
@property
def k_combined(self):
if self.run_mode == 'eigenvalue':
return ufloat(*self._f['k_combined'].value)
return ufloat(*self._f['k_combined'][()])
else:
return None
@property
def k_col_abs(self):
if self.run_mode == 'eigenvalue':
return self._f['k_col_abs'].value
return self._f['k_col_abs'][()]
else:
return None
@property
def k_col_tra(self):
if self.run_mode == 'eigenvalue':
return self._f['k_col_tra'].value
return self._f['k_col_tra'][()]
else:
return None
@property
def k_abs_tra(self):
if self.run_mode == 'eigenvalue':
return self._f['k_abs_tra'].value
return self._f['k_abs_tra'][()]
else:
return None
@ -303,22 +303,22 @@ class StatePoint(object):
@property
def n_batches(self):
return self._f['n_batches'].value
return self._f['n_batches'][()]
@property
def n_inactive(self):
if self.run_mode == 'eigenvalue':
return self._f['n_inactive'].value
return self._f['n_inactive'][()]
else:
return None
@property
def n_particles(self):
return self._f['n_particles'].value
return self._f['n_particles'][()]
@property
def n_realizations(self):
return self._f['n_realizations'].value
return self._f['n_realizations'][()]
@property
def path(self):
@ -330,20 +330,20 @@ class StatePoint(object):
@property
def run_mode(self):
return self._f['run_mode'].value.decode()
return self._f['run_mode'][()].decode()
@property
def runtime(self):
return {name: dataset.value
return {name: dataset[()]
for name, dataset in self._f['runtime'].items()}
@property
def seed(self):
return self._f['seed'].value
return self._f['seed'][()]
@property
def source(self):
return self._f['source_bank'].value if self.source_present else None
return self._f['source_bank'][()] if self.source_present else None
@property
def source_present(self):
@ -376,24 +376,24 @@ class StatePoint(object):
group = tallies_group['tally {}'.format(tally_id)]
# Read the number of realizations
n_realizations = group['n_realizations'].value
n_realizations = group['n_realizations'][()]
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id)
tally._sp_filename = self._f.filename
tally.name = group['name'].value.decode() if 'name' in group else ''
tally.estimator = group['estimator'].value.decode()
tally.name = group['name'][()].decode() if 'name' in group else ''
tally.estimator = group['estimator'][()].decode()
tally.num_realizations = n_realizations
# Read derivative information.
if 'derivative' in group:
deriv_id = group['derivative'].value
deriv_id = group['derivative'][()]
tally.derivative = self.tally_derivatives[deriv_id]
# Read all filters
n_filters = group['n_filters'].value
n_filters = group['n_filters'][()]
if n_filters > 0:
filter_ids = group['filters'].value
filter_ids = group['filters'][()]
filters_group = self._f['tallies/filters']
for filter_id in filter_ids:
filter_group = filters_group['filter {}'.format(
@ -403,15 +403,15 @@ class StatePoint(object):
tally.filters.append(new_filter)
# Read nuclide bins
nuclide_names = group['nuclides'].value
nuclide_names = group['nuclides'][()]
# Add all nuclides to the Tally
for name in nuclide_names:
nuclide = openmc.Nuclide(name.decode().strip())
tally.nuclides.append(nuclide)
scores = group['score_bins'].value
n_score_bins = group['n_score_bins'].value
scores = group['score_bins'][()]
n_score_bins = group['n_score_bins'][()]
# Add the scores to the Tally
for j, score in enumerate(scores):
@ -445,14 +445,14 @@ class StatePoint(object):
group = self._f['tallies/derivatives/derivative {}'
.format(d_id)]
deriv = openmc.TallyDerivative(derivative_id=d_id)
deriv.variable = group['independent variable'].value.decode()
deriv.variable = group['independent variable'][()].decode()
if deriv.variable == 'density':
deriv.material = group['material'].value
deriv.material = group['material'][()]
elif deriv.variable == 'nuclide_density':
deriv.material = group['material'].value
deriv.nuclide = group['nuclide'].value.decode()
deriv.material = group['material'][()]
deriv.nuclide = group['nuclide'][()].decode()
elif deriv.variable == 'temperature':
deriv.material = group['material'].value
deriv.material = group['material'][()]
self._derivs[d_id] = deriv
self._derivs_read = True

24
openmc/summary.py
View file

@ -85,14 +85,14 @@ class Summary(object):
def _read_nuclides(self):
if 'nuclides/names' in self._f:
names = self._f['nuclides/names'].value
awrs = self._f['nuclides/awrs'].value
names = self._f['nuclides/names'][()]
awrs = self._f['nuclides/awrs'][()]
for name, awr in zip(names, awrs):
self._nuclides[name.decode()] = awr
def _read_macroscopics(self):
if 'macroscopics/names' in self._f:
names = self._f['macroscopics/names'].value
names = self._f['macroscopics/names'][()]
for name in names:
self._macroscopics = name.decode()
@ -130,34 +130,34 @@ class Summary(object):
for key, group in self._f['geometry/cells'].items():
cell_id = int(key.lstrip('cell '))
name = group['name'].value.decode() if 'name' in group else ''
fill_type = group['fill_type'].value.decode()
name = group['name'][()].decode() if 'name' in group else ''
fill_type = group['fill_type'][()].decode()
if fill_type == 'material':
fill = group['material'].value
fill = group['material'][()]
elif fill_type == 'universe':
fill = group['fill'].value
fill = group['fill'][()]
else:
fill = group['lattice'].value
fill = group['lattice'][()]
region = group['region'].value.decode() if 'region' in group else ''
region = group['region'][()].decode() if 'region' in group else ''
# Create this Cell
cell = openmc.Cell(cell_id=cell_id, name=name)
if fill_type == 'universe':
if 'translation' in group:
translation = group['translation'][...]
translation = group['translation'][()]
translation = np.asarray(translation, dtype=np.float64)
cell.translation = translation
if 'rotation' in group:
rotation = group['rotation'][...]
rotation = group['rotation'][()]
rotation = np.asarray(rotation, dtype=np.int)
cell._rotation = rotation
elif fill_type == 'material':
cell.temperature = group['temperature'][...]
cell.temperature = group['temperature'][()]
# Store Cell fill information for after Universe/Lattice creation
cell_fills[cell.id] = (fill_type, fill)

516
openmc/surface.py
View file

@ -1,9 +1,9 @@
from abc import ABCMeta
from abc import ABCMeta, abstractmethod
from collections import OrderedDict
from copy import deepcopy
from functools import partial
from numbers import Real, Integral
from xml.etree import ElementTree as ET
from warnings import warn
import numpy as np
@ -14,8 +14,13 @@ from openmc.mixin import IDManagerMixin
_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic']
_WARNING_UPPER = """\
"{}(...) accepts an argument named '{}', not '{}'. Future versions of OpenMC \
will not accept the capitalized version.\
"""
class Surface(IDManagerMixin):
class Surface(IDManagerMixin, metaclass=ABCMeta):
"""An implicit surface with an associated boundary condition.
An implicit surface is defined as the set of zeros of a function of the
@ -142,7 +147,6 @@ class Surface(IDManagerMixin):
desired half-space
"""
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
@ -175,6 +179,14 @@ class Surface(IDManagerMixin):
return memo[self]
@abstractmethod
def evaluate(self, point):
pass
@abstractmethod
def translate(self, vector):
pass
def to_xml_element(self):
"""Return XML representation of the surface
@ -257,9 +269,9 @@ class Surface(IDManagerMixin):
"""
surface_id = int(group.name.split('/')[-1].lstrip('surface '))
name = group['name'].value.decode() if 'name' in group else ''
surf_type = group['type'].value.decode()
bc = group['boundary_type'].value.decode()
name = group['name'][()].decode() if 'name' in group else ''
surf_type = group['type'][()].decode()
bc = group['boundary_type'][()].decode()
coeffs = group['coefficients'][...]
# Create the Surface based on its type
@ -280,29 +292,29 @@ class Surface(IDManagerMixin):
surface = Plane(surface_id, bc, A, B, C, D, name)
elif surf_type == 'x-cylinder':
y0, z0, R = coeffs
surface = XCylinder(surface_id, bc, y0, z0, R, name)
y0, z0, r = coeffs
surface = XCylinder(surface_id, bc, y0, z0, r, name)
elif surf_type == 'y-cylinder':
x0, z0, R = coeffs
surface = YCylinder(surface_id, bc, x0, z0, R, name)
x0, z0, r = coeffs
surface = YCylinder(surface_id, bc, x0, z0, r, name)
elif surf_type == 'z-cylinder':
x0, y0, R = coeffs
surface = ZCylinder(surface_id, bc, x0, y0, R, name)
x0, y0, r = coeffs
surface = ZCylinder(surface_id, bc, x0, y0, r, name)
elif surf_type == 'sphere':
x0, y0, z0, R = coeffs
surface = Sphere(surface_id, bc, x0, y0, z0, R, name)
x0, y0, z0, r = coeffs
surface = Sphere(surface_id, bc, x0, y0, z0, r, name)
elif surf_type in ['x-cone', 'y-cone', 'z-cone']:
x0, y0, z0, R2 = coeffs
x0, y0, z0, r2 = coeffs
if surf_type == 'x-cone':
surface = XCone(surface_id, bc, x0, y0, z0, R2, name)
surface = XCone(surface_id, bc, x0, y0, z0, r2, name)
elif surf_type == 'y-cone':
surface = YCone(surface_id, bc, x0, y0, z0, R2, name)
surface = YCone(surface_id, bc, x0, y0, z0, r2, name)
elif surf_type == 'z-cone':
surface = ZCone(surface_id, bc, x0, y0, z0, R2, name)
surface = ZCone(surface_id, bc, x0, y0, z0, r2, name)
elif surf_type == 'quadric':
a, b, c, d, e, f, g, h, j, k = coeffs
@ -324,13 +336,13 @@ class Plane(Surface):
Boundary condition that defines the behavior for particles hitting the
surface. Defaults to transmissive boundary condition where particles
freely pass through the surface.
A : float, optional
a : float, optional
The 'A' parameter for the plane. Defaults to 1.
B : float, optional
b : float, optional
The 'B' parameter for the plane. Defaults to 0.
C : float, optional
c : float, optional
The 'C' parameter for the plane. Defaults to 0.
D : float, optional
d : float, optional
The 'D' parameter for the plane. Defaults to 0.
name : str, optional
Name of the plane. If not specified, the name will be the empty string.
@ -363,56 +375,61 @@ class Plane(Surface):
"""
_type = 'plane'
_coeff_keys = ('A', 'B', 'C', 'D')
_coeff_keys = ('a', 'b', 'c', 'd')
def __init__(self, surface_id=None, boundary_type='transmission',
A=1., B=0., C=0., D=0., name=''):
a=1., b=0., c=0., d=0., name='', **kwargs):
super().__init__(surface_id, boundary_type, name=name)
self._periodic_surface = None
self.a = A
self.b = B
self.c = C
self.d = D
self.a = a
self.b = b
self.c = c
self.d = d
for k, v in kwargs.items():
if k in 'ABCD':
warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k),
FutureWarning)
setattr(self, k.lower(), v)
@property
def a(self):
return self.coefficients['A']
return self.coefficients['a']
@property
def b(self):
return self.coefficients['B']
return self.coefficients['b']
@property
def c(self):
return self.coefficients['C']
return self.coefficients['c']
@property
def d(self):
return self.coefficients['D']
return self.coefficients['d']
@property
def periodic_surface(self):
return self._periodic_surface
@a.setter
def a(self, A):
check_type('A coefficient', A, Real)
self._coefficients['A'] = A
def a(self, a):
check_type('A coefficient', a, Real)
self._coefficients['a'] = a
@b.setter
def b(self, B):
check_type('B coefficient', B, Real)
self._coefficients['B'] = B
def b(self, b):
check_type('B coefficient', b, Real)
self._coefficients['b'] = b
@c.setter
def c(self, C):
check_type('C coefficient', C, Real)
self._coefficients['C'] = C
def c(self, c):
check_type('C coefficient', c, Real)
self._coefficients['c'] = c
@d.setter
def d(self, D):
check_type('D coefficient', D, Real)
self._coefficients['D'] = D
def d(self, d):
check_type('D coefficient', d, Real)
self._coefficients['d'] = d
@periodic_surface.setter
def periodic_surface(self, periodic_surface):
@ -432,13 +449,34 @@ class Plane(Surface):
Returns
-------
float
:math:`Ax' + By' + Cz' - d`
:math:`Ax' + By' + Cz' - D`
"""
x, y, z = point
return self.a*x + self.b*y + self.c*z - self.d
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.Plane
Translated surface
"""
vx, vy, vz = vector
d = self.d + self.a*vx + self.b*vy + self.c*vz
if d == self.d:
return self
else:
return type(self)(a=self.a, b=self.b, c=self.c, d=d)
def to_xml_element(self):
"""Return XML representation of the surface
@ -456,6 +494,38 @@ class Plane(Surface):
element.set("periodic_surface_id", str(self.periodic_surface.id))
return element
@classmethod
def from_points(cls, p1, p2, p3, **kwargs):
"""Return a plane given three points that pass through it.
Parameters
----------
p1, p2, p3 : 3-tuples
Points that pass through the plane
kwargs : dict
Keyword arguments passed to the :class:`Plane` constructor
Returns
-------
Plane
Plane that passes through the three points
"""
# Convert to numpy arrays
p1 = np.asarray(p1)
p2 = np.asarray(p2)
p3 = np.asarray(p3)
# Find normal vector to plane by taking cross product of two vectors
# connecting p1->p2 and p1->p3
n = np.cross(p2 - p1, p3 - p1)
# The equation of the plane will by n·(<x,y,z> - p1) = 0. Determine
# coefficients a, b, c, and d based on that
a, b, c = n
d = np.dot(n, p1)
return cls(a=a, b=b, c=c, d=d, **kwargs)
class XPlane(Plane):
"""A plane perpendicular to the x axis of the form :math:`x - x_0 = 0`
@ -561,6 +631,26 @@ class XPlane(Plane):
"""
return point[0] - self.x0
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.XPlane
Translated surface
"""
vx = vector[0]
if vx == 0:
return self
else:
return type(self)(x0=self.x0 + vx)
class YPlane(Plane):
"""A plane perpendicular to the y axis of the form :math:`y - y_0 = 0`
@ -667,6 +757,26 @@ class YPlane(Plane):
"""
return point[1] - self.y0
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.YPlane
Translated surface
"""
vy = vector[1]
if vy == 0.0:
return self
else:
return type(self)(y0=self.y0 + vy)
class ZPlane(Plane):
"""A plane perpendicular to the z axis of the form :math:`z - z_0 = 0`
@ -773,8 +883,28 @@ class ZPlane(Plane):
"""
return point[2] - self.z0
def translate(self, vector):
"""Translate surface in given direction
class Cylinder(Surface, metaclass=ABCMeta):
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.ZPlane
Translated surface
"""
vz = vector[2]
if vz == 0.0:
return self
else:
return type(self)(z0=self.z0 + vz)
class Cylinder(Surface):
"""A cylinder whose length is parallel to the x-, y-, or z-axis.
Parameters
@ -786,7 +916,7 @@ class Cylinder(Surface, metaclass=ABCMeta):
Boundary condition that defines the behavior for particles hitting the
surface. Defaults to transmissive boundary condition where particles
freely pass through the surface.
R : float, optional
r : float, optional
Radius of the cylinder. Defaults to 1.
name : str, optional
Name of the cylinder. If not specified, the name will be the empty
@ -810,23 +940,23 @@ class Cylinder(Surface, metaclass=ABCMeta):
"""
def __init__(self, surface_id=None, boundary_type='transmission',
R=1., name=''):
r=1., name=''):
super().__init__(surface_id, boundary_type, name=name)
self.r = R
self.r = r
@property
def r(self):
return self.coefficients['R']
return self.coefficients['r']
@r.setter
def r(self, R):
check_type('R coefficient', R, Real)
self._coefficients['R'] = R
def r(self, r):
check_type('r coefficient', r, Real)
self._coefficients['r'] = r
class XCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the x-axis of the form
:math:`(y - y_0)^2 + (z - z_0)^2 = R^2`.
:math:`(y - y_0)^2 + (z - z_0)^2 = r^2`.
Parameters
----------
@ -841,7 +971,7 @@ class XCylinder(Cylinder):
y-coordinate of the center of the cylinder. Defaults to 0.
z0 : float, optional
z-coordinate of the center of the cylinder. Defaults to 0.
R : float, optional
r : float, optional
Radius of the cylinder. Defaults to 0.
name : str, optional
Name of the cylinder. If not specified, the name will be the empty
@ -868,11 +998,14 @@ class XCylinder(Cylinder):
"""
_type = 'x-cylinder'
_coeff_keys = ('y0', 'z0', 'R')
_coeff_keys = ('y0', 'z0', 'r')
def __init__(self, surface_id=None, boundary_type='transmission',
y0=0., z0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, R, name=name)
y0=0., z0=0., r=1., name='', *, R=None):
if R is not None:
warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning)
r = R
super().__init__(surface_id, boundary_type, r, name=name)
self.y0 = y0
self.z0 = z0
@ -938,17 +1071,39 @@ class XCylinder(Cylinder):
Returns
-------
float
:math:`(y' - y_0)^2 + (z' - z_0)^2 - R^2`
:math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2`
"""
y = point[1] - self.y0
z = point[2] - self.z0
return y**2 + z**2 - self.r**2
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.XCylinder
Translated surface
"""
vx, vy, vz = vector
if vy == 0.0 and vz == 0.0:
return self
else:
y0 = self.y0 + vy
z0 = self.z0 + vz
return type(self)(y0=y0, z0=z0, r=self.r)
class YCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the y-axis of the form
:math:`(x - x_0)^2 + (z - z_0)^2 = R^2`.
:math:`(x - x_0)^2 + (z - z_0)^2 = r^2`.
Parameters
----------
@ -963,7 +1118,7 @@ class YCylinder(Cylinder):
x-coordinate of the center of the cylinder. Defaults to 0.
z0 : float, optional
z-coordinate of the center of the cylinder. Defaults to 0.
R : float, optional
r : float, optional
Radius of the cylinder. Defaults to 1.
name : str, optional
Name of the cylinder. If not specified, the name will be the empty
@ -990,11 +1145,14 @@ class YCylinder(Cylinder):
"""
_type = 'y-cylinder'
_coeff_keys = ('x0', 'z0', 'R')
_coeff_keys = ('x0', 'z0', 'r')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., z0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, R, name=name)
x0=0., z0=0., r=1., name='', *, R=None):
if R is not None:
warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning)
r = R
super().__init__(surface_id, boundary_type, r, name=name)
self.x0 = x0
self.z0 = z0
@ -1060,17 +1218,39 @@ class YCylinder(Cylinder):
Returns
-------
float
:math:`(x' - x_0)^2 + (z' - z_0)^2 - R^2`
:math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2`
"""
x = point[0] - self.x0
z = point[2] - self.z0
return x**2 + z**2 - self.r**2
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.YCylinder
Translated surface
"""
vx, vy, vz = vector
if vx == 0.0 and vz == 0.0:
return self
else:
x0 = self.x0 + vx
z0 = self.z0 + vz
return type(self)(x0=x0, z0=z0, r=self.r)
class ZCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the z-axis of the form
:math:`(x - x_0)^2 + (y - y_0)^2 = R^2`.
:math:`(x - x_0)^2 + (y - y_0)^2 = r^2`.
Parameters
----------
@ -1085,7 +1265,7 @@ class ZCylinder(Cylinder):
x-coordinate of the center of the cylinder. Defaults to 0.
y0 : float, optional
y-coordinate of the center of the cylinder. Defaults to 0.
R : float, optional
r : float, optional
Radius of the cylinder. Defaults to 1.
name : str, optional
Name of the cylinder. If not specified, the name will be the empty
@ -1112,11 +1292,14 @@ class ZCylinder(Cylinder):
"""
_type = 'z-cylinder'
_coeff_keys = ('x0', 'y0', 'R')
_coeff_keys = ('x0', 'y0', 'r')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, R, name=name)
x0=0., y0=0., r=1., name='', *, R=None):
if R is not None:
warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning)
r = R
super().__init__(surface_id, boundary_type, r, name=name)
self.x0 = x0
self.y0 = y0
@ -1182,16 +1365,38 @@ class ZCylinder(Cylinder):
Returns
-------
float
:math:`(x' - x_0)^2 + (y' - y_0)^2 - R^2`
:math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2`
"""
x = point[0] - self.x0
y = point[1] - self.y0
return x**2 + y**2 - self.r**2
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.ZCylinder
Translated surface
"""
vx, vy, vz = vector
if vx == 0.0 and vy == 0.0:
return self
else:
x0 = self.x0 + vx
y0 = self.y0 + vy
return type(self)(x0=x0, y0=y0, r=self.r)
class Sphere(Surface):
"""A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`.
"""A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`.
Parameters
----------
@ -1208,7 +1413,7 @@ class Sphere(Surface):
y-coordinate of the center of the sphere. Defaults to 0.
z0 : float, optional
z-coordinate of the center of the sphere. Defaults to 0.
R : float, optional
r : float, optional
Radius of the sphere. Defaults to 1.
name : str, optional
Name of the sphere. If not specified, the name will be the empty string.
@ -1238,15 +1443,18 @@ class Sphere(Surface):
"""
_type = 'sphere'
_coeff_keys = ('x0', 'y0', 'z0', 'R')
_coeff_keys = ('x0', 'y0', 'z0', 'r')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R=1., name=''):
x0=0., y0=0., z0=0., r=1., name='', *, R=None):
if R is not None:
warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning)
r = R
super().__init__(surface_id, boundary_type, name=name)
self.x0 = x0
self.y0 = y0
self.z0 = z0
self.r = R
self.r = r
@property
def x0(self):
@ -1262,7 +1470,7 @@ class Sphere(Surface):
@property
def r(self):
return self.coefficients['R']
return self.coefficients['r']
@x0.setter
def x0(self, x0):
@ -1280,9 +1488,9 @@ class Sphere(Surface):
self._coefficients['z0'] = z0
@r.setter
def r(self, R):
check_type('R coefficient', R, Real)
self._coefficients['R'] = R
def r(self, r):
check_type('r coefficient', r, Real)
self._coefficients['r'] = r
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
@ -1329,7 +1537,7 @@ class Sphere(Surface):
Returns
-------
float
:math:`(x' - x_0)^2 + (y' - y_0)^2 + (z' - z_0)^2 - R^2`
:math:`(x' - x_0)^2 + (y' - y_0)^2 + (z' - z_0)^2 - r^2`
"""
x = point[0] - self.x0
@ -1337,8 +1545,31 @@ class Sphere(Surface):
z = point[2] - self.z0
return x**2 + y**2 + z**2 - self.r**2
def translate(self, vector):
"""Translate surface in given direction
class Cone(Surface, metaclass=ABCMeta):
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.Sphere
Translated surface
"""
vx, vy, vz = vector
if vx == 0.0 and vy == 0.0 and vz == 0.0:
return self
else:
x0 = self.x0 + vx
y0 = self.y0 + vy
z0 = self.z0 + vz
return type(self)(x0=x0, y0=y0, z0=z0, r=self.r)
class Cone(Surface):
"""A conical surface parallel to the x-, y-, or z-axis.
Parameters
@ -1356,7 +1587,7 @@ class Cone(Surface, metaclass=ABCMeta):
y-coordinate of the apex. Defaults to 0.
z0 : float
z-coordinate of the apex. Defaults to 0.
R2 : float
r2 : float
Parameter related to the aperature. Defaults to 1.
name : str
Name of the cone. If not specified, the name will be the empty string.
@ -1385,15 +1616,18 @@ class Cone(Surface, metaclass=ABCMeta):
"""
_coeff_keys = ('x0', 'y0', 'z0', 'R2')
_coeff_keys = ('x0', 'y0', 'z0', 'r2')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
x0=0., y0=0., z0=0., r2=1., name='', *, R2=None):
if R2 is not None:
warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning)
r2 = R2
super().__init__(surface_id, boundary_type, name=name)
self.x0 = x0
self.y0 = y0
self.z0 = z0
self.r2 = R2
self.r2 = r2
@property
def x0(self):
@ -1409,7 +1643,7 @@ class Cone(Surface, metaclass=ABCMeta):
@property
def r2(self):
return self.coefficients['R2']
return self.coefficients['r2']
@x0.setter
def x0(self, x0):
@ -1427,14 +1661,37 @@ class Cone(Surface, metaclass=ABCMeta):
self._coefficients['z0'] = z0
@r2.setter
def r2(self, R2):
check_type('R^2 coefficient', R2, Real)
self._coefficients['R2'] = R2
def r2(self, r2):
check_type('r^2 coefficient', r2, Real)
self._coefficients['r2'] = r2
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.Cone
Translated surface
"""
vx, vy, vz = vector
if vx == 0.0 and vy == 0.0 and vz == 0.0:
return self
else:
x0 = self.x0 + vx
y0 = self.y0 + vy
z0 = self.z0 + vz
return type(self)(x0=x0, y0=y0, z0=z0, r2=self.r2)
class XCone(Cone):
"""A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
R^2 (x - x_0)^2`.
r^2 (x - x_0)^2`.
Parameters
----------
@ -1451,7 +1708,7 @@ class XCone(Cone):
y-coordinate of the apex. Defaults to 0.
z0 : float, optional
z-coordinate of the apex. Defaults to 0.
R2 : float, optional
r2 : float, optional
Parameter related to the aperature. Defaults to 1.
name : str, optional
Name of the cone. If not specified, the name will be the empty string.
@ -1464,7 +1721,7 @@ class XCone(Cone):
y-coordinate of the apex
z0 : float
z-coordinate of the apex
R2 : float
r2 : float
Parameter related to the aperature
boundary_type : {'transmission, 'vacuum', 'reflective'}
Boundary condition that defines the behavior for particles hitting the
@ -1494,7 +1751,7 @@ class XCone(Cone):
Returns
-------
float
:math:`(y' - y_0)^2 + (z' - z_0)^2 - R^2(x' - x_0)^2`
:math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2(x' - x_0)^2`
"""
x = point[0] - self.x0
@ -1505,7 +1762,7 @@ class XCone(Cone):
class YCone(Cone):
"""A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 =
R^2 (y - y_0)^2`.
r^2 (y - y_0)^2`.
Parameters
----------
@ -1522,7 +1779,7 @@ class YCone(Cone):
y-coordinate of the apex. Defaults to 0.
z0 : float, optional
z-coordinate of the apex. Defaults to 0.
R2 : float, optional
r2 : float, optional
Parameter related to the aperature. Defaults to 1.
name : str, optional
Name of the cone. If not specified, the name will be the empty string.
@ -1535,7 +1792,7 @@ class YCone(Cone):
y-coordinate of the apex
z0 : float
z-coordinate of the apex
R2 : float
r2 : float
Parameter related to the aperature
boundary_type : {'transmission, 'vacuum', 'reflective'}
Boundary condition that defines the behavior for particles hitting the
@ -1565,7 +1822,7 @@ class YCone(Cone):
Returns
-------
float
:math:`(x' - x_0)^2 + (z' - z_0)^2 - R^2(y' - y_0)^2`
:math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2(y' - y_0)^2`
"""
x = point[0] - self.x0
@ -1576,7 +1833,7 @@ class YCone(Cone):
class ZCone(Cone):
"""A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 =
R^2 (z - z_0)^2`.
r^2 (z - z_0)^2`.
Parameters
----------
@ -1593,7 +1850,7 @@ class ZCone(Cone):
y-coordinate of the apex. Defaults to 0.
z0 : float, optional
z-coordinate of the apex. Defaults to 0.
R2 : float, optional
r2 : float, optional
Parameter related to the aperature. Defaults to 1.
name : str, optional
Name of the cone. If not specified, the name will be the empty string.
@ -1606,7 +1863,7 @@ class ZCone(Cone):
y-coordinate of the apex
z0 : float
z-coordinate of the apex
R2 : float
r2 : float
Parameter related to the aperature
boundary_type : {'transmission, 'vacuum', 'reflective'}
Boundary condition that defines the behavior for particles hitting the
@ -1636,7 +1893,7 @@ class ZCone(Cone):
Returns
-------
float
:math:`(x' - x_0)^2 + (y' - y_0)^2 - R^2(z' - z_0)^2`
:math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2(z' - z_0)^2`
"""
x = point[0] - self.x0
@ -1810,6 +2067,30 @@ class Quadric(Surface):
y*(self.b*y + self.e*z + self.h) + \
z*(self.c*z + self.f*x + self.j) + self.k
def translate(self, vector):
"""Translate surface in given direction
Parameters
----------
vector : iterable of float
Direction in which surface should be translated
Returns
-------
openmc.Quadric
Translated surface
"""
vx, vy, vz = vector
a, b, c, d, e, f, g, h, j, k = (getattr(self, key) for key in
self._coeff_keys)
k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz
- g*vx - h*vy - j*vz)
g = g - 2*a*vx - d*vy - f*vz
h = h - 2*b*vy - d*vx - e*vz
j = j - 2*c*vz - e*vy - f*vx
return type(self)(a=a, b=b, c=c, d=d, e=e, f=f, g=g, h=h, j=j, k=k)
class Halfspace(Region):
"""A positive or negative half-space region.
@ -1824,7 +2105,7 @@ class Halfspace(Region):
can be created from an existing Surface through the __neg__ and __pos__
operators, as the following example demonstrates:
>>> sphere = openmc.Sphere(surface_id=1, R=10.0)
>>> sphere = openmc.Sphere(surface_id=1, r=10.0)
>>> inside_sphere = -sphere
>>> outside_sphere = +sphere
>>> type(inside_sphere)
@ -1955,3 +2236,30 @@ class Halfspace(Region):
clone = deepcopy(self)
clone.surface = self.surface.clone(memo)
return clone
def translate(self, vector, memo=None):
"""Translate half-space in given direction
Parameters
----------
vector : iterable of float
Direction in which region should be translated
memo : dict or None
Dictionary used for memoization
Returns
-------
openmc.Halfspace
Translated half-space
"""
if memo is None:
memo = {}
# If translated surface not in memo, add it
key = (self.surface, tuple(vector))
if key not in memo:
memo[key] = self.surface.translate(vector)
# Return translated surface
return type(self)(memo[key], self.side)

11
openmc/tallies.py
View file

@ -5,6 +5,7 @@ from functools import partial, reduce
from itertools import product
from numbers import Integral, Real
import operator
from pathlib import Path
import warnings
from xml.etree import ElementTree as ET
@ -216,7 +217,7 @@ class Tally(IDManagerMixin):
f = h5py.File(self._sp_filename, 'r')
# Extract Tally data from the file
data = f['tallies/tally {0}/results'.format(self.id)].value
data = f['tallies/tally {0}/results'.format(self.id)]
sum = data[:, :, 0]
sum_sq = data[:, :, 1]
@ -3189,7 +3190,11 @@ class Tallies(cv.CheckedList):
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)
# Check if path is a directory
p = Path(path)
if p.is_dir():
p /= 'tallies.xml'
# Write the XML Tree to the tallies.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True,
encoding='utf-8', method="xml")
tree.write(str(p), xml_declaration=True, encoding='utf-8')

2
openmc/universe.py
View file

@ -124,7 +124,7 @@ class Universe(IDManagerMixin):
"""
universe_id = int(group.name.split('/')[-1].lstrip('universe '))
cell_ids = group['cells'].value
cell_ids = group['cells'][()]
# Create this Universe
universe = cls(universe_id)

6
openmc/volume.py
View file

@ -211,9 +211,9 @@ class VolumeCalculation(object):
domain_id = int(obj_name[7:])
ids.append(domain_id)
group = f[obj_name]
volume = ufloat(*group['volume'].value)
nucnames = group['nuclides'].value
atoms_ = group['atoms'].value
volume = ufloat(*group['volume'][()])
nucnames = group['nuclides'][()]
atoms_ = group['atoms'][()]
atom_dict = OrderedDict()
for name_i, atoms_i in zip(nucnames, atoms_):

34
scripts/openmc-plot-mesh-tally
View file

@ -216,24 +216,40 @@ class MeshPlotter(tk.Frame):
index = self.filterBoxes[f.short_name].current()
spec_list.append((type(f), (index,)))
dims = (self.nx, self.ny, self.nz)
text = self.basisBox.get()
if text == 'xy':
dims = (self.nx, self.ny)
h_ind = 0
v_ind = 1
elif text == 'yz':
dims = (self.ny, self.nz)
h_ind = 1
v_ind = 2
else:
dims = (self.nx, self.nz)
h_ind = 0
v_ind = 2
axial_ind = 3 - (h_ind + v_ind)
dims = (dims[h_ind], dims[v_ind])
mesh_dim = len(self.mesh.dimension)
if mesh_dim == 3:
mesh_indices = [0,0,0]
else:
mesh_indices = [0,0]
matrix = np.zeros(dims)
for i in range(dims[0]):
for j in range(dims[1]):
if text == 'xy':
meshtuple = (i + 1, j + 1, axial_level)
elif text == 'yz':
meshtuple = (axial_level, i + 1, j + 1)
if mesh_dim == 3:
mesh_indices[h_ind] = i + 1
mesh_indices[v_ind] = j + 1
mesh_indices[axial_ind] = axial_level
else:
meshtuple = (i + 1, axial_level, j + 1)
mesh_indices[0] = i + 1
mesh_indices[1] = j + 1
filters, filter_bins = zip(*spec_list + [
(type(mesh_filter), (meshtuple,))])
(type(mesh_filter), (tuple(mesh_indices),))])
mean = selectedTally.get_values(
[self.scoreBox.get()], filters, filter_bins)
stdev = selectedTally.get_values(

5
src/bank.cpp
View file

@ -16,10 +16,9 @@ namespace openmc {
namespace simulation {
int64_t n_bank;
std::vector<Particle::Bank> source_bank;
std::vector<Particle::Bank> fission_bank;
std::vector<Particle::Bank> secondary_bank;
#ifdef _OPENMP
std::vector<Particle::Bank> master_fission_bank;
#endif
@ -33,7 +32,7 @@ std::vector<Particle::Bank> master_fission_bank;
void free_memory_bank()
{
simulation::source_bank.clear();
#pragma omp parallel
#pragma omp parallel
{
simulation::fission_bank.clear();
}

2
src/bremsstrahlung.cpp
View file

@ -45,7 +45,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
auto n_e = data::ttb_e_grid.size();
// Find the lower bounding index of the incident electron energy
int j = lower_bound_index(data::ttb_e_grid.cbegin(),
size_t j = lower_bound_index(data::ttb_e_grid.cbegin(),
data::ttb_e_grid.cend(), e);
if (j == n_e - 1) --j;

185
src/cell.cpp
View file

@ -2,6 +2,7 @@
#include <cmath>
#include <sstream>
#include <set>
#include <string>
#include "openmc/capi.h"
@ -325,7 +326,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// Convert the infix region spec to RPN.
rpn_ = generate_rpn(id_, region_);
rpn_.shrink_to_fit();
// Check if this is a simple cell.
simple_ = true;
@ -336,6 +336,21 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
}
}
// If this cell is simple, remove all the superfluous operator tokens.
if (simple_) {
size_t i0 = 0;
size_t i1 = 0;
while (i1 < rpn_.size()) {
if (rpn_[i1] < OP_UNION) {
rpn_[i0] = rpn_[i1];
++i0;
}
++i1;
}
rpn_.resize(i0);
}
rpn_.shrink_to_fit();
// Read the translation vector.
if (check_for_node(cell_node, "translation")) {
if (fill_ == C_NONE) {
@ -522,19 +537,17 @@ bool
CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
{
for (int32_t token : rpn_) {
if (token < OP_UNION) {
// If the token is not an operator, evaluate the sense of particle with
// respect to the surface and see if the token matches the sense. If the
// particle's surface attribute is set and matches the token, that
// overrides the determination based on sense().
if (token == on_surface) {
} else if (-token == on_surface) {
return false;
} else {
// Note the off-by-one indexing
bool sense = model::surfaces[abs(token)-1]->sense(r, u);
if (sense != (token > 0)) {return false;}
}
// Assume that no tokens are operators. Evaluate the sense of particle with
// respect to the surface and see if the token matches the sense. If the
// particle's surface attribute is set and matches the token, that
// overrides the determination based on sense().
if (token == on_surface) {
} else if (-token == on_surface) {
return false;
} else {
// Note the off-by-one indexing
bool sense = model::surfaces[abs(token)-1]->sense(r, u);
if (sense != (token > 0)) {return false;}
}
}
return true;
@ -601,7 +614,7 @@ std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_id(3, id_);
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
@ -621,7 +634,7 @@ DAGCell::distance(Position r, Direction u, int32_t on_surface) const
bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_id(3, id_);
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
int result = 0;
double pnt[3] = {r.x, r.y, r.z};
@ -635,6 +648,142 @@ void DAGCell::to_hdf5(hid_t group_id) const { return; }
#endif
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
UniversePartitioner::UniversePartitioner(const Universe& univ)
{
// Define an ordered set of surface indices that point to z-planes. Use a
// functor to to order the set by the z0_ values of the corresponding planes.
struct compare_surfs {
bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
{
const auto* surf = model::surfaces[i_surf].get();
const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zi = zplane->z0_;
surf = model::surfaces[j_surf].get();
zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zj = zplane->z0_;
return zi < zj;
}
};
std::set<int32_t, compare_surfs> surf_set;
// Find all of the z-planes in this universe. A set is used here for the
// O(log(n)) insertions that will ensure entries are not repeated.
for (auto i_cell : univ.cells_) {
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
auto i_surf = std::abs(token) - 1;
const auto* surf = model::surfaces[i_surf].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
surf_set.insert(i_surf);
}
}
}
// Populate the surfs_ vector from the ordered set.
surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());
// Populate the partition lists.
partitions_.resize(surfs_.size() + 1);
for (auto i_cell : univ.cells_) {
// Find the tokens for bounding z-planes.
int32_t lower_token = 0, upper_token = 0;
double min_z, max_z;
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
const auto* surf = model::surfaces[std::abs(token) - 1].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {
if (lower_token == 0 || zplane->z0_ < min_z) {
lower_token = token;
min_z = zplane->z0_;
}
if (upper_token == 0 || zplane->z0_ > max_z) {
upper_token = token;
max_z = zplane->z0_;
}
}
}
}
// If there are no bounding z-planes, add this cell to all partitions.
if (lower_token == 0) {
for (auto& p : partitions_) p.push_back(i_cell);
continue;
}
// Find the first partition this cell lies in. If the lower_token indicates
// a negative halfspace, then the cell is unbounded in the lower direction
// and it lies in the first partition onward. Otherwise, it is bounded by
// the positive halfspace given by the lower_token.
int first_partition = 0;
if (lower_token > 0) {
for (int i = 0; i < surfs_.size(); ++i) {
if (lower_token == surfs_[i] + 1) {
first_partition = i + 1;
break;
}
}
}
// Find the last partition this cell lies in. The logic is analogous to the
// logic for first_partition.
int last_partition = surfs_.size();
if (upper_token < 0) {
for (int i = first_partition; i < surfs_.size(); ++i) {
if (upper_token == -(surfs_[i] + 1)) {
last_partition = i;
break;
}
}
}
// Add the cell to all relevant partitions.
for (int i = first_partition; i <= last_partition; ++i) {
partitions_[i].push_back(i_cell);
}
}
}
const std::vector<int32_t>&
UniversePartitioner::get_cells(Position r, Direction u) const
{
// Perform a binary search for the partition containing the given coordinates.
int left = 0;
int middle = (surfs_.size() - 1) / 2;
int right = surfs_.size() - 1;
while (true) {
// Check the sense of the coordinates for the current surface.
const auto& surf = *model::surfaces[surfs_[middle]];
if (surf.sense(r, u)) {
// The coordinates lie in the positive halfspace. Recurse if there are
// more surfaces to check. Otherwise, return the cells on the positive
// side of this surface.
int right_leaf = right - (right - middle) / 2;
if (right_leaf != middle) {
left = middle + 1;
middle = right_leaf;
} else {
return partitions_[middle+1];
}
} else {
// The coordinates lie in the negative halfspace. Recurse if there are
// more surfaces to check. Otherwise, return the cells on the negative
// side of this surface.
int left_leaf = left + (middle - left) / 2;
if (left_leaf != middle) {
right = middle-1;
middle = left_leaf;
} else {
return partitions_[middle];
}
}
}
}
//==============================================================================
// Non-method functions
//==============================================================================
@ -830,9 +979,9 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed)
{
moab::EntityHandle surf =
surf_xed->dagmc_ptr_->entity_by_id(2, surf_xed->id_);
surf_xed->dagmc_ptr_->entity_by_index(2, surf_xed->dag_index_);
moab::EntityHandle vol =
cur_cell->dagmc_ptr_->entity_by_id(3, cur_cell->id_);
cur_cell->dagmc_ptr_->entity_by_index(3, cur_cell->dag_index_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);

1
src/cross_sections.cpp
View file

@ -247,7 +247,6 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
LibraryKey key {Library::Type::photon, element};
int idx = data::library_map[key];
std::string& filename = data::libraries[idx].path_;
int i_element = data::element_map[element];
write_message("Reading " + element + " from " + filename, 6);
// Open file and make sure version is sufficient

70
src/dagmc.cpp
View file

@ -4,10 +4,11 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/material.h"
#include "openmc/string_utils.h"
#include "openmc/settings.h"
#include "openmc/geometry.h"
#ifdef DAGMC
@ -88,13 +89,60 @@ bool write_uwuw_materials_xml() {
return found_uwuw_mats;
}
void legacy_assign_material(const std::string& mat_string, DAGCell* c)
{
bool mat_found_by_name = false;
// attempt to find a material with a matching name
for (const auto& m : model::materials) {
if (mat_string == m->name_) {
// assign the material with that name
if (!mat_found_by_name) {
mat_found_by_name = true;
c->material_.push_back(m->id_);
// report error if more than one material is found
} else {
std::stringstream err_msg;
err_msg << "More than one material found with name " << mat_string
<< ". Please ensure materials have unique names if using this"
<< " property to assign materials.";
fatal_error(err_msg);
}
}
}
// if no material was set using a name, assign by id
if (!mat_found_by_name) {
try {
auto id = std::stoi(mat_string);
c->material_.emplace_back(id);
} catch (const std::invalid_argument&) {
std::stringstream err_msg;
err_msg << "No material " << mat_string
<< " found for volume (cell) " << c->id_;
fatal_error(err_msg);
}
}
if (settings::verbosity >= 10) {
Material* m = model::materials[model::material_map[c->material_[0]]].get();
std::stringstream msg;
msg << "DAGMC material " << mat_string << " was assigned";
if (mat_found_by_name) {
msg << " using material name: " << m->name_;
} else {
msg << " using material id: " << m->id_;
}
write_message(msg.str(), 10);
}
}
void load_dagmc_geometry()
{
if (!model::DAG) {
model::DAG = new moab::DagMC();
}
/// Materials \\\
// --- Materials ---
// create uwuw instance
UWUW uwuw(DAGMC_FILENAME.c_str());
@ -128,7 +176,7 @@ void load_dagmc_geometry()
rval = model::DAG->parse_properties(keywords, dum, delimiters.c_str());
MB_CHK_ERR_CONT(rval);
/// Cells (Volumes) \\\
// --- Cells (Volumes) ---
// initialize cell objects
int n_cells = model::DAG->num_entities(3);
@ -138,7 +186,8 @@ void load_dagmc_geometry()
// set cell ids using global IDs
DAGCell* c = new DAGCell();
c->id_ = model::DAG->id_by_index(3, i+1);
c->dag_index_ = i+1;
c->id_ = model::DAG->id_by_index(3, c->dag_index_);
c->dagmc_ptr_ = model::DAG;
c->universe_ = dagmc_univ_id; // set to zero for now
c->fill_ = C_NONE; // no fill, single universe
@ -187,7 +236,7 @@ void load_dagmc_geometry()
size_t _comp_pos = mat_value.find(_comp);
if (_comp_pos != std::string::npos) { mat_value.erase(_comp_pos, _comp.length()); }
// assign IC material by id
c->material_.push_back(std::stoi(mat_value));
legacy_assign_material(mat_value, c);
}
} else {
// if no material is found, the implicit complement is void
@ -234,9 +283,7 @@ void load_dagmc_geometry()
fatal_error(err_msg);
}
} else {
// if not using UWUW materials, we'll find this material
// later in the materials.xml
c->material_.push_back(std::stoi(mat_value));
legacy_assign_material(mat_value, c);
}
}
}
@ -251,7 +298,7 @@ void load_dagmc_geometry()
"This may result in lost particles and rapid simulation failure.");
}
/// Surfaces \\\
// --- Surfaces ---
// initialize surface objects
int n_surfaces = model::DAG->num_entities(2);
@ -261,7 +308,8 @@ void load_dagmc_geometry()
// set cell ids using global IDs
DAGSurface* s = new DAGSurface();
s->id_ = model::DAG->id_by_index(2, i+1);
s->dag_index_ = i+1;
s->id_ = model::DAG->id_by_index(2, s->dag_index_);
s->dagmc_ptr_ = model::DAG;
// set BCs
@ -303,7 +351,7 @@ void load_dagmc_geometry()
// add to global array and map
model::surfaces.emplace_back(s);
model::surface_map[s->id_] = s->id_;
model::surface_map[s->id_] = i;
}
return;

3
src/distribution_energy.cpp
View file

@ -245,6 +245,9 @@ double ContinuousTabular::sample(double E) const
E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k +
2.0*frac*(r1 - c_k))) - p_l_k)/frac;
}
} else {
throw std::runtime_error{"Unexpected interpolation for continuous energy "
"distribution."};
}
// Now interpolate between incident energy bins i and i + 1

68
src/eigenvalue.cpp
View file

@ -20,9 +20,14 @@
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#include <algorithm> // for min
#include <array>
#include <cmath> // for sqrt, abs, pow
#include <iterator> // for back_inserter
#include <string>
namespace openmc {
@ -81,21 +86,22 @@ void synchronize_bank()
#ifdef OPENMC_MPI
int64_t start = 0;
MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
int64_t n_bank = simulation::fission_bank.size();
MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
// While we would expect the value of start on rank 0 to be 0, the MPI
// standard says that the receive buffer on rank 0 is undefined and not
// significant
if (mpi::rank == 0) start = 0;
int64_t finish = start + simulation::n_bank;
int64_t finish = start + simulation::fission_bank.size();
int64_t total = finish;
MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
#else
int64_t start = 0;
int64_t finish = simulation::n_bank;
int64_t total = simulation::n_bank;
int64_t finish = simulation::fission_bank.size();
int64_t total = simulation::fission_bank.size();
#endif
// If there are not that many particles per generation, it's possible that no
@ -103,7 +109,7 @@ void synchronize_bank()
// extra logic to treat this circumstance, we really want to ensure the user
// runs enough particles to avoid this in the first place.
if (simulation::n_bank == 0) {
if (simulation::fission_bank.empty()) {
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
@ -132,23 +138,23 @@ void synchronize_bank()
// Allocate temporary source bank
int64_t index_temp = 0;
std::vector<Particle::Bank> temp_sites(3*simulation::work);
std::vector<Particle::Bank> temp_sites(3*simulation::work_per_rank);
for (int64_t i = 0; i < simulation::n_bank; ++i) {
for (const auto& site : simulation::fission_bank) {
// If there are less than n_particles particles banked, automatically add
// int(n_particles/total) sites to temp_sites. For example, if you need
// 1000 and 300 were banked, this would add 3 source sites per banked site
// and the remaining 100 would be randomly sampled.
if (total < settings::n_particles) {
for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
temp_sites[index_temp] = simulation::fission_bank[i];
temp_sites[index_temp] = site;
++index_temp;
}
}
// Randomly sample sites needed
if (prn() < p_sample) {
temp_sites[index_temp] = simulation::fission_bank[i];
temp_sites[index_temp] = site;
++index_temp;
}
}
@ -189,7 +195,7 @@ void synchronize_bank()
// fission bank
sites_needed = settings::n_particles - finish;
for (int i = 0; i < sites_needed; ++i) {
int i_bank = simulation::n_bank - sites_needed + i;
int i_bank = simulation::fission_bank.size() - sites_needed + i;
temp_sites[index_temp] = simulation::fission_bank[i_bank];
++index_temp;
}
@ -346,38 +352,32 @@ void calculate_average_keff()
#ifdef _OPENMP
void join_bank_from_threads()
{
// Initialize the total number of fission bank sites
int64_t total = 0;
int n_threads = omp_get_max_threads();
#pragma omp parallel
#pragma omp parallel
{
// Copy thread fission bank sites to one shared copy
#pragma omp for ordered schedule(static)
for (int i = 0; i < simulation::n_threads; ++i) {
#pragma omp ordered
#pragma omp for ordered schedule(static)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
{
std::copy(
&simulation::fission_bank[0],
&simulation::fission_bank[0] + simulation::n_bank,
&simulation::master_fission_bank[total]
simulation::fission_bank.cbegin(),
simulation::fission_bank.cend(),
std::back_inserter(simulation::master_fission_bank)
);
total += simulation::n_bank;
}
}
// Make sure all threads have made it to this point
#pragma omp barrier
#pragma omp barrier
// Now copy the shared fission bank sites back to the master thread's copy.
if (simulation::thread_id == 0) {
simulation::n_bank = total;
std::copy(
&simulation::master_fission_bank[0],
&simulation::master_fission_bank[0] + simulation::n_bank,
&simulation::fission_bank[0]
);
if (omp_get_thread_num() == 0) {
simulation::fission_bank = simulation::master_fission_bank;
simulation::master_fission_bank.clear();
} else {
simulation::n_bank = 0;
simulation::fission_bank.clear();
}
}
}
@ -537,8 +537,8 @@ void shannon_entropy()
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
simulation::fission_bank.data(), 0, nullptr, &sites_outside);
xt::xtensor<double, 1> p = m->count_sites(simulation::fission_bank,
&sites_outside);
// display warning message if there were sites outside entropy box
if (sites_outside) {
@ -577,8 +577,8 @@ void ufs_count_sites()
} else {
// count number of source sites in each ufs mesh cell
bool sites_outside;
simulation::source_frac = m->count_sites(simulation::work,
simulation::source_bank.data(), 0, nullptr, &sites_outside);
simulation::source_frac = m->count_sites(simulation::source_bank,
&sites_outside);
// Check for sites outside of the mesh
if (mpi::master && sites_outside) {
@ -597,7 +597,7 @@ void ufs_count_sites()
// Since the total starting weight is not equal to n_particles, we need to
// renormalize the weight of the source sites
for (int i = 0; i < simulation::work; ++i) {
for (int i = 0; i < simulation::work_per_rank; ++i) {
simulation::source_bank[i].wgt *= settings::n_particles / total;
}
}

3
src/endf.cpp
View file

@ -164,9 +164,8 @@ double Tabulated1D::operator()(double x) const
Interpolation interp;
if (n_regions_ == 0) {
interp = Interpolation::lin_lin;
} else if (n_regions_ == 1) {
} else {
interp = int_[0];
} else if (n_regions_ > 1) {
for (int j = 0; j < n_regions_; ++j) {
if (i < nbt_[j]) {
interp = int_[j];

71
src/geometry.cpp
View file

@ -22,6 +22,7 @@ namespace openmc {
namespace model {
int root_universe {-1};
int n_coord_levels;
std::vector<int64_t> overlap_check_count;
@ -31,15 +32,13 @@ std::vector<int64_t> overlap_check_count;
// Non-member functions
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p)
bool check_cell_overlap(Particle* p)
{
int n_coord = p->n_coord_;
// Loop through each coordinate level
for (int j = 0; j < n_coord; j++) {
Universe& univ = *model::universes[p->coord_[j].universe];
int n = univ.cells_.size();
// Loop through each cell on this level
for (auto index_cell : univ.cells_) {
@ -90,7 +89,12 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
} else {
int i_universe = p->coord_[p->n_coord_-1].universe;
const auto& cells {model::universes[i_universe]->cells_};
const auto& univ {*model::universes[i_universe]};
const auto& cells {
!univ.partitioner_
? model::universes[i_universe]->cells_
: univ.partitioner_->get_cells(p->r_local(), p->u_local())
};
for (auto it = cells.cbegin(); it != cells.cend(); it++) {
i_cell = *it;
@ -245,7 +249,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//==============================================================================
extern "C" bool
bool
find_cell(Particle* p, bool use_neighbor_lists)
{
// Determine universe (if not yet set, use root universe).
@ -257,7 +261,7 @@ find_cell(Particle* p, bool use_neighbor_lists)
}
// Reset all the deeper coordinate levels.
for (int i = p->n_coord_; i < MAX_COORD; i++) {
for (int i = p->n_coord_; i < p->coord_.size(); i++) {
p->coord_[i].reset();
}
@ -287,8 +291,8 @@ find_cell(Particle* p, bool use_neighbor_lists)
//==============================================================================
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3])
void
cross_lattice(Particle* p, const BoundaryInfo& boundary)
{
auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]};
@ -302,9 +306,9 @@ cross_lattice(Particle* p, int lattice_translation[3])
}
// Set the lattice indices.
p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2];
p->coord_[p->n_coord_-1].lattice_x += boundary.lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += boundary.lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += boundary.lattice_translation[2];
std::array<int, 3> i_xyz {p->coord_[p->n_coord_-1].lattice_x,
p->coord_[p->n_coord_-1].lattice_y,
p->coord_[p->n_coord_-1].lattice_z};
@ -345,18 +349,13 @@ cross_lattice(Particle* p, int lattice_translation[3])
//==============================================================================
extern "C" void
distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
int lattice_translation[3], int* next_level)
BoundaryInfo distance_to_boundary(Particle* p)
{
*dist = INFINITY;
BoundaryInfo info;
double d_lat = INFINITY;
double d_surf = INFINITY;
lattice_translation[0] = 0;
lattice_translation[1] = 0;
lattice_translation[2] = 0;
int32_t level_surf_cross;
std::array<int, 3> level_lat_trans;
std::array<int, 3> level_lat_trans {};
// Loop over each coordinate level.
for (int i = 0; i < p->n_coord_; i++) {
@ -401,43 +400,43 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
// If the boundary on this coordinate level is coincident with a boundary on
// a higher level then we need to make sure that the higher level boundary
// is selected. This logic must consider floating point precision.
if (d_surf < d_lat) {
if (*dist == INFINITY || ((*dist) - d_surf)/(*dist) >= FP_REL_PRECISION) {
*dist = d_surf;
double& d = info.distance;
if (d_surf < d_lat - FP_COINCIDENT) {
if (d == INFINITY || (d - d_surf)/d >= FP_REL_PRECISION) {
d = d_surf;
// If the cell is not simple, it is possible that both the negative and
// positive half-space were given in the region specification. Thus, we
// have to explicitly check which half-space the particle would be
// traveling into if the surface is crossed
if (c.simple_) {
*surface_crossed = level_surf_cross;
info.surface_index = level_surf_cross;
} else {
Position r_hit = r + d_surf * u;
Surface& surf {*model::surfaces[std::abs(level_surf_cross)-1]};
Direction norm = surf.normal(r_hit);
if (u.dot(norm) > 0) {
*surface_crossed = std::abs(level_surf_cross);
info.surface_index = std::abs(level_surf_cross);
} else {
*surface_crossed = -std::abs(level_surf_cross);
info.surface_index = -std::abs(level_surf_cross);
}
}
lattice_translation[0] = 0;
lattice_translation[1] = 0;
lattice_translation[2] = 0;
*next_level = i + 1;
info.lattice_translation[0] = 0;
info.lattice_translation[1] = 0;
info.lattice_translation[2] = 0;
info.coord_level = i + 1;
}
} else {
if (*dist == INFINITY || ((*dist) - d_lat)/(*dist) >= FP_REL_PRECISION) {
*dist = d_lat;
*surface_crossed = F90_NONE;
lattice_translation[0] = level_lat_trans[0];
lattice_translation[1] = level_lat_trans[1];
lattice_translation[2] = level_lat_trans[2];
*next_level = i + 1;
if (d == INFINITY || (d - d_lat)/d >= FP_REL_PRECISION) {
d = d_lat;
info.surface_index = 0;
info.lattice_translation = level_lat_trans;
info.coord_level = i + 1;
}
}
}
return info;
}
//==============================================================================

45
src/geometry_aux.cpp
View file

@ -120,6 +120,40 @@ adjust_indices()
}
}
//==============================================================================
//! Partition some universes with many z-planes for faster find_cell searches.
void
partition_universes()
{
// Iterate over universes with more than 10 cells. (Fewer than 10 is likely
// not worth partitioning.)
for (const auto& univ : model::universes) {
if (univ->cells_.size() > 10) {
// Collect the set of surfaces in this universe.
std::unordered_set<int32_t> surf_inds;
for (auto i_cell : univ->cells_) {
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1);
}
}
// Partition the universe if there are more than 5 z-planes. (Fewer than
// 5 is likely not worth it.)
int n_zplanes = 0;
for (auto i_surf : surf_inds) {
if (dynamic_cast<const SurfaceZPlane*>(model::surfaces[i_surf].get())) {
++n_zplanes;
if (n_zplanes > 5) {
univ->partitioner_ = std::make_unique<UniversePartitioner>(*univ);
break;
}
}
}
}
}
}
//==============================================================================
void
@ -210,6 +244,7 @@ void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
// Perform some final operations to set up the geometry
adjust_indices();
count_cell_instances(model::root_universe);
partition_universes();
// Assign temperatures to cells that don't have temperatures already assigned
assign_temperatures();
@ -217,14 +252,8 @@ void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
// Determine desired temperatures for each nuclide and S(a,b) table
get_temperatures(nuc_temps, thermal_temps);
// Check to make sure there are not too many nested coordinate levels in the
// geometry since the coordinate list is statically allocated for performance
// reasons
if (maximum_levels(model::root_universe) > MAX_COORD) {
fatal_error("Too many nested coordinate levels in the geometry. "
"Try increasing the maximum number of coordinate levels by "
"providing the CMake -Dmaxcoord= option.");
}
// Determine number of nested coordinate levels in the geometry
model::n_coord_levels = maximum_levels(model::root_universe);
}
//==============================================================================

5
src/hdf5_interface.cpp
View file

@ -3,6 +3,7 @@
#include <array>
#include <cstring>
#include <sstream>
#include <stdexcept>
#include <string>
#include "xtensor/xtensor.hpp"
@ -47,6 +48,8 @@ get_shape(hid_t obj_id, hsize_t* dims)
dspace = H5Dget_space(obj_id);
} else if (type == H5I_ATTR) {
dspace = H5Aget_space(obj_id);
} else {
throw std::runtime_error{"Expected dataset or attribute in call to get_shape."};
}
H5Sget_simple_extent_dims(dspace, dims, nullptr);
H5Sclose(dspace);
@ -70,6 +73,8 @@ std::vector<hsize_t> object_shape(hid_t obj_id)
dspace = H5Dget_space(obj_id);
} else if (type == H5I_ATTR) {
dspace = H5Aget_space(obj_id);
} else {
throw std::runtime_error{"Expected dataset or attribute in call to object_shape."};
}
int n = H5Sget_simple_extent_ndims(dspace);

7
src/initialize.cpp
View file

@ -121,7 +121,6 @@ void initialize_mpi(MPI_Comm intracomm)
int
parse_command_line(int argc, char* argv[])
{
char buffer[256]; // buffer for reading attribute
int last_flag = 0;
for (int i=1; i < argc; ++i) {
std::string arg {argv[i]};
@ -196,13 +195,13 @@ parse_command_line(int argc, char* argv[])
#ifdef _OPENMP
// Read and set number of OpenMP threads
simulation::n_threads = std::stoi(argv[i]);
if (simulation::n_threads < 1) {
int n_threads = std::stoi(argv[i]);
if (n_threads < 1) {
std::string msg {"Number of threads must be positive."};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
omp_set_num_threads(simulation::n_threads);
omp_set_num_threads(n_threads);
#else
if (mpi::master)
warning("Ignoring number of threads specified on command line.");

87
src/lattice.cpp
View file

@ -6,6 +6,7 @@
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/string_utils.h"
@ -291,7 +292,7 @@ RectLattice::get_indices(Position r, Direction u) const
double ix_ {(r.x - lower_left_.x) / pitch_.x};
long ix_close {std::lround(ix_)};
int ix;
if (std::abs(ix_ - ix_close) < FP_COINCIDENT) {
if (coincident(ix_, ix_close)) {
ix = (u.x > 0) ? ix_close : ix_close - 1;
} else {
ix = std::floor(ix_);
@ -301,7 +302,7 @@ RectLattice::get_indices(Position r, Direction u) const
double iy_ {(r.y - lower_left_.y) / pitch_.y};
long iy_close {std::lround(iy_)};
int iy;
if (std::abs(iy_ - iy_close) < FP_COINCIDENT) {
if (coincident(iy_, iy_close)) {
iy = (u.y > 0) ? iy_close : iy_close - 1;
} else {
iy = std::floor(iy_);
@ -312,7 +313,7 @@ RectLattice::get_indices(Position r, Direction u) const
if (is_3d_) {
double iz_ {(r.z - lower_left_.z) / pitch_.z};
long iz_close {std::lround(iz_)};
if (std::abs(iz_ - iz_close) < FP_COINCIDENT) {
if (coincident(iz_, iz_close)) {
iz = (u.z > 0) ? iz_close : iz_close - 1;
} else {
iz = std::floor(iz_);
@ -710,69 +711,84 @@ const
std::array<int, 3>
HexLattice::get_indices(Position r, Direction u) const
{
// The implementation for HexLattice currently doesn't use direction
// information. As a result, we move the position slightly forward to
// determine what lattice index the particle is most likely to be in.
r += TINY_BIT * u;
// Offset the xyz by the lattice center.
Position r_o {r.x - center_.x, r.y - center_.y, r.z};
if (is_3d_) {r_o.z -= center_.z;}
// Index the z direction.
std::array<int, 3> out;
// Index the z direction, accounting for coincidence
int iz = 0;
if (is_3d_) {
out[2] = std::floor(r_o.z / pitch_[1] + 0.5 * n_axial_);
} else {
out[2] = 0;
double iz_ {r_o.z / pitch_[1] + 0.5 * n_axial_};
long iz_close {std::lround(iz_)};
if (coincident(iz_, iz_close)) {
iz = (u.z > 0) ? iz_close : iz_close - 1;
} else {
iz = std::floor(iz_);
}
}
// Convert coordinates into skewed bases. The (x, alpha) basis is used to
// find the index of the global coordinates to within 4 cells.
double alpha = r_o.y - r_o.x / std::sqrt(3.0);
out[0] = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
out[1] = std::floor(alpha / pitch_[0]);
int ix = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
int ia = std::floor(alpha / pitch_[0]);
// Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
// the array is offset so that the indices never go below 0).
out[0] += n_rings_-1;
out[1] += n_rings_-1;
ix += n_rings_-1;
ia += n_rings_-1;
// Calculate the (squared) distance between the particle and the centers of
// the four possible cells. Regular hexagonal tiles form a Voronoi
// tessellation so the xyz should be in the hexagonal cell that it is closest
// to the center of. This method is used over a method that uses the
// remainders of the floor divisions above because it provides better finite
// precision performance. Squared distances are used becasue they are more
// precision performance. Squared distances are used because they are more
// computationally efficient than normal distances.
int k {1};
int k_min {1};
// COINCIDENCE CHECK
// if a distance to center, d, is within the coincidence tolerance of the
// current minimum distance, d_min, the particle is on an edge or vertex.
// In this case, the dot product of the position vector and direction vector
// for the current indices, dp, and the dot product for the currently selected
// indices, dp_min, are compared. The cell which the particle is moving into
// is kept (i.e. the cell with the lowest dot product as the vectors will be
// completely opposed if the particle is moving directly toward the center of
// the cell).
int ix_chg {};
int ia_chg {};
double d_min {INFTY};
double dp_min {INFTY};
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
const std::array<int, 3> i_xyz {out[0] + j, out[1] + i, 0};
// get local coordinates
const std::array<int, 3> i_xyz {ix + j, ia + i, 0};
Position r_t = get_local_position(r, i_xyz);
// calculate distance
double d = r_t.x*r_t.x + r_t.y*r_t.y;
if (d < d_min) {
// check for coincidence
bool on_boundary = coincident(d, d_min);
if (d < d_min || on_boundary) {
// normalize r_t and find dot product
r_t /= std::sqrt(d);
double dp = u.x * r_t.x + u.y * r_t.y;
// do not update values if particle is on a
// boundary and not moving into this cell
if (on_boundary && dp > dp_min) continue;
// update values
d_min = d;
k_min = k;
ix_chg = j;
ia_chg = i;
dp_min = dp;
}
k++;
}
}
// Select the minimum squared distance which corresponds to the cell the
// coordinates are in.
if (k_min == 2) {
++out[0];
} else if (k_min == 3) {
++out[1];
} else if (k_min == 4) {
++out[0];
++out[1];
}
// update outgoing indices
ix += ix_chg;
ia += ia_chg;
return out;
return {ix, ia, iz};
}
//==============================================================================
@ -800,7 +816,6 @@ HexLattice::is_valid_index(int indx) const
{
int nx {2*n_rings_ - 1};
int ny {2*n_rings_ - 1};
int nz {n_axial_};
int iz = indx / (nx * ny);
int iy = (indx - nx*ny*iz) / nx;
int ix = indx - nx*ny*iz - nx*iy;

295
src/material.cpp
View file

@ -451,6 +451,101 @@ void Material::init_thermal()
thermal_tables_ = tables;
}
void Material::collision_stopping_power(double* s_col, bool positron)
{
// Average electron number and average atomic weight
double electron_density = 0.0;
double mass_density = 0.0;
// Log of the mean excitation energy of the material
double log_I = 0.0;
// Effective number of conduction electrons in the material
double n_conduction = 0.0;
// Oscillator strength and square of the binding energy for each oscillator
// in material
std::vector<double> f;
std::vector<double> e_b_sq;
for (int i = 0; i < element_.size(); ++i) {
const auto& elm = data::elements[element_[i]];
double awr = data::nuclides[nuclide_[i]]->awr_;
// Get atomic density of nuclide given atom/weight percent
double atom_density = (atom_density_[0] > 0.0) ?
atom_density_[i] : -atom_density_[i] / awr;
electron_density += atom_density * elm.Z_;
mass_density += atom_density * awr * MASS_NEUTRON;
log_I += atom_density * elm.Z_ * std::log(elm.I_);
for (int j = 0; j < elm.n_electrons_.size(); ++j) {
if (elm.n_electrons_[j] < 0) {
n_conduction -= elm.n_electrons_[j] * atom_density;
continue;
}
e_b_sq.push_back(elm.ionization_energy_[j] * elm.ionization_energy_[j]);
f.push_back(elm.n_electrons_[j] * atom_density);
}
}
log_I /= electron_density;
n_conduction /= electron_density;
for (auto& f_i : f) f_i /= electron_density;
// Get density in g/cm^3 if it is given in atom/b-cm
double density = (density_ < 0.0) ? -density_ : mass_density / N_AVOGADRO;
// Calculate the square of the plasma energy
double e_p_sq = PLANCK_C * PLANCK_C * PLANCK_C * N_AVOGADRO *
electron_density * density / (2.0 * PI * PI * FINE_STRUCTURE *
MASS_ELECTRON_EV * mass_density);
// Get the Sternheimer adjustment factor
double rho = sternheimer_adjustment(f, e_b_sq, e_p_sq, n_conduction, log_I,
1.0e-6, 100);
// Classical electron radius in cm
constexpr double CM_PER_ANGSTROM {1.0e-8};
constexpr double r_e = CM_PER_ANGSTROM * PLANCK_C / (2.0 * PI *
FINE_STRUCTURE * MASS_ELECTRON_EV);
// Constant in expression for collision stopping power
constexpr double BARN_PER_CM_SQ {1.0e24};
double c = BARN_PER_CM_SQ * 2.0 * PI * r_e * r_e * MASS_ELECTRON_EV *
electron_density;
// Loop over incident charged particle energies
for (int i = 0; i < data::ttb_e_grid.size(); ++i) {
double E = data::ttb_e_grid(i);
// Get the density effect correction
double delta = density_effect(f, e_b_sq, e_p_sq, n_conduction, rho, E,
1.0e-6, 100);
// Square of the ratio of the speed of light to the velocity of the charged
// particle
double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) / ((E + MASS_ELECTRON_EV)
* (E + MASS_ELECTRON_EV));
double tau = E / MASS_ELECTRON_EV;
double F;
if (positron) {
double t = tau + 2.0;
F = std::log(4.0) - (beta_sq / 12.0) * (23.0 + 14.0 / t + 10.0 / (t * t)
+ 4.0 / (t * t * t));
} else {
F = (1.0 - beta_sq) * (1.0 + tau * tau / 8.0 - (2.0 * tau + 1.0) *
std::log(2.0));
}
// Calculate the collision stopping power for this energy
s_col[i] = c / beta_sq * (2.0 * (std::log(E) - log_I) + std::log(1.0 + tau
/ 2.0) + F - delta);
}
}
void Material::init_bremsstrahlung()
{
// Create new object
@ -481,16 +576,11 @@ void Material::init_bremsstrahlung()
double Z_eq_sq = 0.0;
double sum_density = 0.0;
// Calculate the molecular DCS and the molecular total stopping power using
// Get the collision stopping power of the material
this->collision_stopping_power(stopping_power_collision.data(), positron);
// Calculate the molecular DCS and the molecular radiative stopping power using
// Bragg's additivity rule.
// TODO: The collision stopping power cannot be accurately calculated using
// Bragg's additivity rule since the mean excitation energies and the
// density effect corrections cannot simply be summed together. Bragg's
// additivity rule fails especially when a higher-density compound is
// composed of elements that are in lower-density form at normal temperature
// and pressure (at which the NIST stopping powers are given). It will be
// used to approximate the collision stopping powers for now, but should be
// fixed in the future.
for (int i = 0; i < n; ++i) {
// Get pointer to current element
const auto& elm = data::elements[element_[i]];
@ -499,7 +589,6 @@ void Material::init_bremsstrahlung()
// Get atomic density and mass density of nuclide given atom/weight percent
double atom_density = (atom_density_[0] > 0.0) ?
atom_density_[i] : -atom_density_[i] / awr;
double mass_density = atom_density * awr;
// Calculate the "equivalent" atomic number Zeq of the material
Z_eq_sq += atom_density * elm.Z_ * elm.Z_;
@ -508,13 +597,8 @@ void Material::init_bremsstrahlung()
// Accumulate material DCS
dcs += (atom_density * elm.Z_ * elm.Z_) * elm.dcs_;
// Accumulate material collision stopping power
stopping_power_collision += (mass_density * MASS_NEUTRON / N_AVOGADRO)
* elm.stopping_power_collision_;
// Accumulate material radiative stopping power
stopping_power_radiative += (mass_density * MASS_NEUTRON / N_AVOGADRO)
* elm.stopping_power_radiative_;
stopping_power_radiative += atom_density * elm.stopping_power_radiative_;
}
Z_eq_sq /= sum_density;
@ -569,12 +653,13 @@ void Material::init_bremsstrahlung()
// photon energy k
double x = x_l + (k - k_l)*(x_r - x_l)/(k_r - k_l);
// Ratio of the velocity of the charged particle to the speed of light
double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) /
(e + MASS_ELECTRON_EV);
// Square of the ratio of the speed of light to the velocity of the
// charged particle
double beta_sq = e * (e + 2.0 * MASS_ELECTRON_EV) / ((e +
MASS_ELECTRON_EV) * (e + MASS_ELECTRON_EV));
// Compute the integrand of the PDF
f(j) = x / (beta*beta * stopping_power(j) * w);
f(j) = x / (beta_sq * stopping_power(j) * w);
}
// Number of points to integrate
@ -644,13 +729,13 @@ void Material::init_nuclide_index()
}
}
void Material::calculate_xs(const Particle& p) const
void Material::calculate_xs(Particle& p) const
{
// Set all material macroscopic cross sections to zero
simulation::material_xs.total = 0.0;
simulation::material_xs.absorption = 0.0;
simulation::material_xs.fission = 0.0;
simulation::material_xs.nu_fission = 0.0;
p.macro_xs_.total = 0.0;
p.macro_xs_.absorption = 0.0;
p.macro_xs_.fission = 0.0;
p.macro_xs_.nu_fission = 0.0;
if (p.type_ == Particle::Type::neutron) {
this->calculate_neutron_xs(p);
@ -659,7 +744,7 @@ void Material::calculate_xs(const Particle& p) const
}
}
void Material::calculate_neutron_xs(const Particle& p) const
void Material::calculate_neutron_xs(Particle& p) const
{
// Find energy index on energy grid
int neutron = static_cast<int>(Particle::Type::neutron);
@ -707,13 +792,12 @@ void Material::calculate_neutron_xs(const Particle& p) const
int i_nuclide = nuclide_[i];
// Calculate microscopic cross section for this nuclide
const auto& micro {simulation::micro_xs[i_nuclide]};
const auto& micro {p.neutron_xs_[i_nuclide]};
if (p.E_ != micro.last_E
|| p.sqrtkT_ != micro.last_sqrtkT
|| i_sab != micro.index_sab
|| sab_frac != micro.sab_frac) {
data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E_, i_grid,
p.sqrtkT_, sab_frac);
data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p);
}
// ======================================================================
@ -723,19 +807,19 @@ void Material::calculate_neutron_xs(const Particle& p) const
double atom_density = atom_density_(i);
// Add contributions to cross sections
simulation::material_xs.total += atom_density * micro.total;
simulation::material_xs.absorption += atom_density * micro.absorption;
simulation::material_xs.fission += atom_density * micro.fission;
simulation::material_xs.nu_fission += atom_density * micro.nu_fission;
p.macro_xs_.total += atom_density * micro.total;
p.macro_xs_.absorption += atom_density * micro.absorption;
p.macro_xs_.fission += atom_density * micro.fission;
p.macro_xs_.nu_fission += atom_density * micro.nu_fission;
}
}
void Material::calculate_photon_xs(const Particle& p) const
void Material::calculate_photon_xs(Particle& p) const
{
simulation::material_xs.coherent = 0.0;
simulation::material_xs.incoherent = 0.0;
simulation::material_xs.photoelectric = 0.0;
simulation::material_xs.pair_production = 0.0;
p.macro_xs_.coherent = 0.0;
p.macro_xs_.incoherent = 0.0;
p.macro_xs_.photoelectric = 0.0;
p.macro_xs_.pair_production = 0.0;
// Add contribution from each nuclide in material
for (int i = 0; i < nuclide_.size(); ++i) {
@ -746,9 +830,9 @@ void Material::calculate_photon_xs(const Particle& p) const
int i_element = element_[i];
// Calculate microscopic cross section for this nuclide
const auto& micro {simulation::micro_photon_xs[i_element]};
const auto& micro {p.photon_xs_[i_element]};
if (p.E_ != micro.last_E) {
data::elements[i_element].calculate_xs(p.E_);
data::elements[i_element].calculate_xs(p);
}
// ========================================================================
@ -758,11 +842,11 @@ void Material::calculate_photon_xs(const Particle& p) const
double atom_density = atom_density_(i);
// Add contributions to material macroscopic cross sections
simulation::material_xs.total += atom_density * micro.total;
simulation::material_xs.coherent += atom_density * micro.coherent;
simulation::material_xs.incoherent += atom_density * micro.incoherent;
simulation::material_xs.photoelectric += atom_density * micro.photoelectric;
simulation::material_xs.pair_production += atom_density * micro.pair_production;
p.macro_xs_.total += atom_density * micro.total;
p.macro_xs_.coherent += atom_density * micro.coherent;
p.macro_xs_.incoherent += atom_density * micro.incoherent;
p.macro_xs_.photoelectric += atom_density * micro.photoelectric;
p.macro_xs_.pair_production += atom_density * micro.pair_production;
}
}
@ -866,6 +950,124 @@ void Material::to_hdf5(hid_t group) const
// Non-method functions
//==============================================================================
double sternheimer_adjustment(const std::vector<double>& f, const
std::vector<double>& e_b_sq, double e_p_sq, double n_conduction, double
log_I, double tol, int max_iter)
{
// Get the total number of oscillators
int n = f.size();
// Calculate the Sternheimer adjustment factor using Newton's method
double rho = 2.0;
int iter;
for (iter = 0; iter < max_iter; ++iter) {
double rho_0 = rho;
// Function to find the root of and its derivative
double g = 0.0;
double gp = 0.0;
for (int i = 0; i < n; ++i) {
// Square of resonance energy of a bound-shell oscillator
double e_r_sq = e_b_sq[i] * rho * rho + 2.0 / 3.0 * f[i] * e_p_sq;
g += f[i] * std::log(e_r_sq);
gp += e_b_sq[i] * f[i] * rho / e_r_sq;
}
// Include conduction electrons
if (n_conduction > 0.0) {
g += n_conduction * std::log(n_conduction * e_p_sq);
}
// Set the next guess: rho_n+1 = rho_n - g(rho_n)/g'(rho_n)
rho -= (g - 2.0 * log_I) / (2.0 * gp);
// If the initial guess is too large, rho can be negative
if (rho < 0.0) rho = rho_0 / 2.0;
// Check for convergence
if (std::abs(rho - rho_0) / rho_0 < tol) break;
}
// Did not converge
if (iter >= max_iter) {
warning("Maximum Newton-Raphson iterations exceeded.");
rho = 1.0e-6;
}
return rho;
}
double density_effect(const std::vector<double>& f, const std::vector<double>&
e_b_sq, double e_p_sq, double n_conduction, double rho, double E, double tol,
int max_iter)
{
// Get the total number of oscillators
int n = f.size();
// Square of the ratio of the speed of light to the velocity of the charged
// particle
double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) / ((E + MASS_ELECTRON_EV) *
(E + MASS_ELECTRON_EV));
// For nonmetals, delta = 0 for beta < beta_0, where beta_0 is obtained by
// setting the frequency w = 0.
double beta_0_sq = 0.0;
if (n_conduction == 0.0) {
for (int i = 0; i < n; ++i) {
beta_0_sq += f[i] * e_p_sq / (e_b_sq[i] * rho * rho);
}
beta_0_sq = 1.0 / (1.0 + beta_0_sq);
}
double delta = 0.0;
if (beta_sq < beta_0_sq) return delta;
// Compute the square of the frequency w^2 using Newton's method, with the
// initial guess of w^2 equal to beta^2 * gamma^2
double w_sq = E / MASS_ELECTRON_EV * (E / MASS_ELECTRON_EV + 2);
int iter;
for (iter = 0; iter < max_iter; ++iter) {
double w_sq_0 = w_sq;
// Function to find the root of and its derivative
double g = 0.0;
double gp = 0.0;
for (int i = 0; i < n; ++i) {
double c = e_b_sq[i] * rho * rho / e_p_sq + w_sq;
g += f[i] / c;
gp -= f[i] / (c * c);
}
// Include conduction electrons
g += n_conduction / w_sq;
gp -= n_conduction / (w_sq * w_sq);
// Set the next guess: w_n+1 = w_n - g(w_n)/g'(w_n)
w_sq -= (g + 1.0 - 1.0 / beta_sq) / gp;
// If the initial guess is too large, w can be negative
if (w_sq < 0.0) w_sq = w_sq_0 / 2.0;
// Check for convergence
if (std::abs(w_sq - w_sq_0) / w_sq_0 < tol) break;
}
// Did not converge
if (iter >= max_iter) {
warning("Maximum Newton-Raphson iterations exceeded: setting density "
"effect correction to zero.");
return delta;
}
// Solve for the density effect correction
for (int i = 0; i < n; ++i) {
double l_sq = e_b_sq[i] * rho * rho / e_p_sq + 2.0 / 3.0 * f[i];
delta += f[i] * std::log((l_sq + w_sq)/l_sq);
}
// Include conduction electrons
if (n_conduction > 0.0) {
delta += n_conduction * std::log((n_conduction + w_sq) / n_conduction);
}
return delta - w_sq * (1.0 - beta_sq);
}
void read_materials_xml()
{
write_message("Reading materials XML file...", 5);
@ -1081,12 +1283,11 @@ openmc_material_set_densities(int32_t index, int n, const char** name, const dou
}
// Set total density to the sum of the vector
int err = mat->set_density(sum_density, "atom/b-cm");
// Assign S(a,b) tables
mat->init_thermal();
return 0;
return err;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;

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