Merge is inevitable because of new test created and test suite was changed

Merge branch 'develop' into test_surface_tallies
This commit is contained in:
guillaume 2017-08-06 14:15:48 -04:00
commit c0efbf13c5
155 changed files with 15153 additions and 12882 deletions

View file

@ -9,12 +9,23 @@ addons:
packages:
- gfortran
- g++
- mpich
- libmpich-dev
cache:
directories:
- $HOME/mpich_install
- $HOME/hdf5_install
- $HOME/phdf5_install
- $HOME/nndc_hdf5
env:
- OPENMC_CONFIG="check_source"
- OPENMC_CONFIG="^hdf5-debug$"
- OPENMC_CONFIG="^omp-hdf5-debug$"
- OPENMC_CONFIG="^mpi-hdf5-debug$"
- OPENMC_CONFIG="^phdf5-debug$"
# We aren't testing the check_source script so just run it with Python 3.
matrix:
exclude:
- python: "2.7"
env: OPENMC_CONFIG="check_source"
before_install:
# ============== Handle Python third-party packages ==============
@ -29,30 +40,37 @@ before_install:
- conda config --set always_yes yes --set changeps1 no
- conda update -q conda
- conda info -a
- conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION six numpy scipy h5py=2.5 pandas
- source activate test-environment
# Install GCC, MPICH, HDF5, PHDF5
- ./tests/travis_install.sh
- export FC=gfortran
- export MPI_DIR=$HOME/mpich_install
- export PHDF5_DIR=$HOME/phdf5_install
- export HDF5_DIR=$HOME/hdf5_install
- if [[ $OPENMC_CONFIG != "check_source" ]]; then
conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION six numpy scipy h5py=2.5 pandas;
source activate test-environment;
sudo add-apt-repository ppa:nschloe/hdf5-backports -y;
sudo apt-get update -q;
sudo apt-get install libhdf5-serial-dev libhdf5-mpich-dev -y;
export FC=gfortran;
export MPI_DIR=/usr;
export PHDF5_DIR=/usr;
export HDF5_DIR=/usr;
fi
install: true
before_script:
- if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/a6sw2cep34wlz6b9i9jwiotaqoayxcxt.xz -O - | tar -C $HOME -xvJ;
- if [[ $OPENMC_CONFIG != "check_source" ]]; then
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/a6sw2cep34wlz6b9i9jwiotaqoayxcxt.xz -O - | tar -C $HOME -xvJ;
fi;
export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml;
git clone --branch=master git://github.com/smharper/windowed_multipole_library.git wmp_lib;
tar xzvf wmp_lib/multipole_lib.tar.gz;
export OPENMC_MULTIPOLE_LIBRARY=$PWD/multipole_lib;
fi
- export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
- git clone --branch=master git://github.com/smharper/windowed_multipole_library.git wmp_lib
- tar xzvf wmp_lib/multipole_lib.tar.gz
- export OPENMC_MULTIPOLE_LIBRARY=$PWD/multipole_lib
script:
- cd tests
- export OMP_NUM_THREADS=2
- ./travis.sh
- if [[ $OPENMC_CONFIG == "check_source" ]]; then
./check_source.py;
else
./run_tests.py -C $OPENMC_CONFIG -j 2;
fi
- cd ..

View file

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 2.8 FATAL_ERROR)
cmake_minimum_required(VERSION 2.8.12 FATAL_ERROR)
project(openmc Fortran C CXX)
# Setup output directories
@ -21,6 +21,11 @@ if (${UNIX})
add_definitions(-DUNIX)
endif()
# Set MACOSX_RPATH
if(POLICY CMP0042)
cmake_policy(SET CMP0042 NEW)
endif()
#===============================================================================
# Command line options
#===============================================================================
@ -62,10 +67,12 @@ endif()
# this, we check for the environment variable HDF5_ROOT and if it exists, use it
# to check whether its a parallel version.
if(DEFINED ENV{HDF5_ROOT} AND EXISTS $ENV{HDF5_ROOT}/bin/h5pcc)
set(HDF5_PREFER_PARALLEL TRUE)
else()
set(HDF5_PREFER_PARALLEL FALSE)
if(NOT DEFINED HDF5_PREFER_PARALLEL)
if(DEFINED ENV{HDF5_ROOT} AND EXISTS $ENV{HDF5_ROOT}/bin/h5pcc)
set(HDF5_PREFER_PARALLEL TRUE)
else()
set(HDF5_PREFER_PARALLEL FALSE)
endif()
endif()
find_package(HDF5 COMPONENTS Fortran_HL)
@ -98,6 +105,8 @@ endif()
# endif()
#endif()
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
# Make sure version is sufficient
execute_process(COMMAND ${CMAKE_Fortran_COMPILER} -dumpversion
@ -107,60 +116,47 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
endif()
# GCC compiler options
list(APPEND f90flags -cpp -std=f2008 -fbacktrace -O2)
list(APPEND cflags -cpp -std=c99 -O2)
list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2)
if(debug)
list(REMOVE_ITEM f90flags -O2)
list(REMOVE_ITEM cflags -O2)
list(APPEND f90flags -g -Wall -pedantic -fbounds-check
-ffpe-trap=invalid,overflow,underflow)
list(APPEND cflags -g -Wall -pedantic -fbounds-check)
list(APPEND ldflags -g)
endif()
if(profile)
list(APPEND f90flags -pg)
list(APPEND cflags -pg)
list(APPEND ldflags -pg)
endif()
if(optimize)
list(REMOVE_ITEM f90flags -O2)
list(REMOVE_ITEM cflags -O2)
list(APPEND f90flags -O3)
list(APPEND cflags -O3)
endif()
if(openmp)
list(APPEND f90flags -fopenmp)
list(APPEND cflags -fopenmp)
list(APPEND ldflags -fopenmp)
endif()
if(coverage)
list(APPEND f90flags -coverage)
list(APPEND cflags -coverage)
list(APPEND ldflags -coverage)
endif()
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
# Intel compiler options
list(APPEND f90flags -fpp -std08 -assume byterecl -traceback)
list(APPEND cflags -std=c99)
if(debug)
list(APPEND f90flags -g -warn -ftrapuv -fp-stack-check
"-check all" -fpe0 -O0)
list(APPEND cflags -g -w3 -ftrapuv -fp-stack-check -O0)
list(APPEND ldflags -g)
endif()
if(profile)
list(APPEND f90flags -pg)
list(APPEND cflags -pg)
list(APPEND ldflags -pg)
endif()
if(optimize)
list(APPEND f90flags -O3)
list(APPEND cflags -O3)
endif()
if(openmp)
list(APPEND f90flags -qopenmp)
list(APPEND cflags -qopenmp)
list(APPEND ldflags -qopenmp)
endif()
@ -212,6 +208,49 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Cray)
endif()
if(CMAKE_C_COMPILER_ID STREQUAL GNU)
# GCC compiler options
list(APPEND cflags -std=c99 -O2)
if(debug)
list(REMOVE_ITEM cflags -O2)
list(APPEND cflags -g -Wall -pedantic -fbounds-check)
endif()
if(profile)
list(APPEND cflags -pg)
endif()
if(optimize)
list(REMOVE_ITEM cflags -O2)
list(APPEND cflags -O3)
endif()
if(coverage)
list(APPEND cflags -coverage)
endif()
elseif(CMAKE_C_COMPILER_ID STREQUAL Intel)
# Intel compiler options
list(APPEND cflags -std=c99)
if(debug)
list(APPEND cflags -g -w3 -ftrapuv -fp-stack-check -O0)
endif()
if(profile)
list(APPEND cflags -pg)
endif()
if(optimize)
list(APPEND cflags -O3)
endif()
elseif(CMAKE_C_COMPILER_ID MATCHES Clang)
# Clang options
list(APPEND cflags -std=c99)
if(debug)
list(APPEND cflags -g -O0 -ftrapv)
endif()
if(optimize)
list(APPEND cflags -O3)
endif()
endif()
# Show flags being used
message(STATUS "Fortran flags: ${f90flags}")
message(STATUS "C flags: ${cflags}")
@ -262,6 +301,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/angle_distribution.F90
src/angleenergy_header.F90
src/bank_header.F90
src/api.F90
src/cmfd_data.F90
src/cmfd_execute.F90
src/cmfd_header.F90
@ -279,9 +319,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/endf.F90
src/endf_header.F90
src/energy_distribution.F90
src/energy_grid.F90
src/error.F90
src/finalize.F90
src/geometry.F90
src/geometry_header.F90
src/global.F90
@ -344,34 +382,18 @@ set(LIBOPENMC_FORTRAN_SRC
src/volume_calc.F90
src/volume_header.F90
src/xml_interface.F90)
add_library(libopenmc STATIC ${LIBOPENMC_FORTRAN_SRC})
add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC})
set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc)
add_executable(${program} src/main.F90)
set_property(TARGET ${program} libopenmc pugixml_fortran
PROPERTY LINKER_LANGUAGE Fortran)
# target_include_directories was added in CMake 2.8.11 and is the recommended
# way to set include directories. For lesser versions, we revert to set_property
if(CMAKE_VERSION VERSION_LESS 2.8.11)
include_directories(${HDF5_INCLUDE_DIRS})
else()
target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS})
endif()
target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS})
# target_compile_options was added in CMake 2.8.12 and is the recommended way to
# set compile flags. Note that this sets the COMPILE_OPTIONS property (also
# available only in 2.8.12+) rather than the COMPILE_FLAGS property, which is
# deprecated. The former can handle lists whereas the latter cannot.
if (CMAKE_VERSION VERSION_LESS 2.8.12)
string(REPLACE ";" " " f90flags "${f90flags}")
string(REPLACE ";" " " cflags "${cflags}")
set_property(TARGET ${program} PROPERTY COMPILE_FLAGS "${f90flags}")
set_property(TARGET faddeeva PROPERTY COMPILE_FLAGS "${cflags}")
else()
target_compile_options(${program} PUBLIC ${f90flags})
target_compile_options(libopenmc PUBLIC ${f90flags})
target_compile_options(faddeeva PRIVATE ${cflags})
endif()
# set compile flags via target_compile_options
target_compile_options(${program} PUBLIC ${f90flags})
target_compile_options(libopenmc PUBLIC ${f90flags})
target_compile_options(faddeeva PRIVATE ${cflags})
# Add HDF5 library directories to link line with -L
foreach(LIBDIR ${HDF5_LIBRARY_DIRS})
@ -383,6 +405,16 @@ endforeach()
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml_fortran faddeeva)
target_link_libraries(${program} ${ldflags} libopenmc)
#===============================================================================
# Python package
#===============================================================================
add_custom_command(TARGET libopenmc POST_BUILD
COMMAND ${CMAKE_COMMAND} -E copy
$<TARGET_FILE:libopenmc>
${CMAKE_CURRENT_SOURCE_DIR}/openmc/capi/_$<TARGET_FILE_NAME:libopenmc>
COMMENT "Copying libopenmc to Python module directory")
#===============================================================================
# Install executable, scripts, manpage, license
#===============================================================================

261
docs/source/capi/index.rst Normal file
View file

@ -0,0 +1,261 @@
.. _capi:
=====
C API
=====
.. c:function:: void openmc_calculate_voumes()
Run a stochastic volume calculation
.. c:function:: int openmc_cell_get_id(int32_t index, int32_t* id)
Get the ID of a cell
:param index: Index in the cells array
:type index: int32_t
:param id: ID of the cell
:type id: int32_t*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_set_temperature(index index, double T, int32_t* instance)
Set the temperature of a cell.
:param index: Index in the cells array
:type index: int32_t
:param T: Temperature in Kelvin
:type T: double
:param instance: Which instance of the cell. To set the temperature for all
instances, pass a null pointer.
:type instance: int32_t*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: void openmc_finalize()
Finalize a simulation
.. c:function:: void openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance)
Determine the ID of the cell/material containing a given point
:param xyz: Cartesian coordinates
:type xyz: double[3]
:param rtype: Which ID to return (1=cell, 2=material)
:type rtype: int
:param id: ID of the cell/material found. If a material is requested and the
point is in a void, the ID is 0. If an error occurs, the ID is -1.
:type id: int32_t*
:param instance: If a cell is repetaed in the geometry, the instance of the
cell that was found and zero otherwise.
:type instance: int32_t*
.. c:function:: int openmc_get_cell(int32_t id, int32_t* index)
Get the index in the cells array for a cell with a given ID
:param id: ID of the cell
:type id: int32_t
:param index: Index in the cells array
:type index: int32_t*
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_keff(double k_combined[])
:param k_combined: Combined estimate of k-effective
:type k_combined: double[2]
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_nuclide(char name[], int* index)
Get the index in the nuclides array for a nuclide with a given name
:param name: Name of the nuclide
:type name: char[]
:param index: Index in the nuclides array
:type index: int*
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_tally(int32_t id, int32_t* index)
Get the index in the tallies array for a tally with a given ID
:param id: ID of the tally
:type id: int32_t
:param index: Index in the tallies array
:type index: int32_t*
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_material(int32_t id, int32_t* index)
Get the index in the materials array for a material with a given ID
:param id: ID of the material
:type id: int32_t
:param index: Index in the materials array
:type index: int32_t*
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_hard_reset()
Reset tallies, timers, and pseudo-random number generator state
.. c:function:: void openmc_init(int intracomm)
Initialize OpenMC
:param intracomm: MPI intracommunicator
:type intracomm: int
.. c:function:: int openmc_load_nuclide(char name[])
Load data for a nuclide from the HDF5 data library.
:param name: Name of the nuclide.
:type name: char[]
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_material_add_nuclide(int32_t index, char name[], double density)
Add a nuclide to an existing material. If the nuclide already exists, the
density is overwritten.
:param index: Index in the materials array
:type index: int32_t
:param name: Name of the nuclide
:type name: char[]
:param density: Density in atom/b-cm
:type density: double
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_material_get_densities(int32_t index, int* nuclides[], double* densities[])
Get density for each nuclide in a material.
:param index: Index in the materials array
:type index: int32_t
:param nuclides: Pointer to array of nuclide indices
:type nuclides: int**
:param densities: Pointer to the array of densities
:type densities: double**
:param n: Length of the array
:type n: int
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_material_get_id(int32_t index, int32_t* id)
Get the ID of a material
:param index: Index in the materials array
:type index: int32_t
:param id: ID of the material
:type id: int32_t*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_material_set_density(int32_t index, double density)
Set the density of a material.
:param index: Index in the materials array
:type index: int32_t
:param density: Density of the material in atom/b-cm
:type density: double
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_material_set_densities(int32_t, n, char* name[], double density[])
:param index: Index in the materials array
:type index: int32_t
:param n: Length of name/density
:type n: int
:param name: Array of nuclide names
:type name: char**
:param density: Array of densities
:type density: double[]
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_nuclide_name(int index, char* name[])
Get name of a nuclide
:param index: Index in the nuclides array
:type index: int
:param name: Name of the nuclide
:type name: char**
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_plot_geometry()
Run plotting mode.
.. c:function:: void openmc_reset()
Resets all tally scores
.. c:function:: void openmc_run()
Run a simulation
.. c:function:: int openmc_tally_get_id(int32_t index, int32_t* id)
Get the ID of a tally
:param index: Index in the tallies array
:type index: int32_t
:param id: ID of the tally
:type id: int32_t*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_tally_get_nuclides(int32_t index, int* nuclides[], int* n)
Get nuclides specified in a tally
:param index: Index in the tallies array
:type index: int32_t
:param nuclides: Array of nuclide indices
:type nuclides: int**
:param n: Number of nuclides
:type n: int*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_tally_results(int32_t index, double** ptr, int shape_[3])
Get a pointer to tally results array.
:param index: Index in the tallies array
:type index: int32_t
:param ptr: Pointer to the results array
:type ptr: double**
:param shape_: Shape of the results array
:type shape_: int[3]
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_tally_set_nuclides(int32_t index, int n, char* nuclides[])
Set the nuclides for a tally
:param index: Index in the tallies array
:type index: int32_t
:param n: Number of nuclides
:type n: int
:param nuclides: Array of nuclide names
:type nuclides: char**
:return: Return status (negative if an error occurred)
:rtype: int

View file

@ -25,8 +25,9 @@ except ImportError:
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'h5py', 'pandas', 'uncertainties', 'openmoc',
'openmc.data.reconstruct']
'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.interpolate',
'scipy.integrate', 'scipy.optimize', 'scipy.special', 'h5py',
'pandas', 'uncertainties', 'openmoc', 'openmc.data.reconstruct']
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
import numpy as np

View file

@ -36,6 +36,7 @@ free to send a message to the User's Group `mailing list`_.
usersguide/index
devguide/index
pythonapi/index
capi/index
io_formats/index
publications
license

View file

@ -107,9 +107,13 @@ Each ``material`` element can have the following attributes or sub-elements:
multi-group :ref:`energy_mode`.
:sab:
Associates an S(a,b) table with the material. This element has one
Associates an S(a,b) table with the material. This element has an
attribute/sub-element called ``name``. The ``name`` attribute
is the name of the S(a,b) table that should be associated with the material.
There is also an optional ``fraction`` element which indicates what fraction
of the relevant nuclides will be affected by the S(a,b) table (e.g. which
fraction of a material is crystalline versus amorphous). ``fraction``
defaults to unity.
*Default*: None

View file

@ -661,6 +661,17 @@ methods like "nearest" and "interpolation" in the resolved resonance range.
*Default*: False
-------------------------------
``<temperature_range>`` Element
-------------------------------
The ``<temperature_range>`` element specifies a minimum and maximum temperature
in Kelvin above and below which cross sections should be loaded for all nuclides
and thermal scattering tables. This can be used for multi-physics simulations
where the temperatures might change from one iteration to the next.
*Default*: None
.. _temperature_tolerance:
-----------------------------------

View file

@ -0,0 +1,39 @@
---------------------------------------------------
:data:`openmc.capi` -- Python bindings to the C API
---------------------------------------------------
.. automodule:: openmc.capi
Functions
---------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.capi.calculate_volumes
openmc.capi.finalize
openmc.capi.find_cell
openmc.capi.find_material
openmc.capi.hard_reset
openmc.capi.init
openmc.capi.keff
openmc.capi.load_nuclide
openmc.capi.plot_geometry
openmc.capi.reset
openmc.capi.run
openmc.capi.run_in_memory
Classes
-------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.capi.CellView
openmc.capi.MaterialView
openmc.capi.NuclideView
openmc.capi.TallyView

View file

@ -47,5 +47,6 @@ Modules
mgxs
model
data
capi
examples
openmoc

View file

@ -19,7 +19,7 @@ surface is a locus of zeros of a function of Cartesian coordinates
:math:`x,y,z`, e.g.
- A plane perpendicular to the :math:`x` axis: :math:`x - x_0 = 0`
- A cylinder perpendicular to the :math:`z` axis: :math:`(x - x_0)^2 + (y -
- A cylinder parallel to the :math:`z` axis: :math:`(x - x_0)^2 + (y -
y_0)^2 - R^2 = 0`
- A sphere: :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0`

View file

@ -401,14 +401,6 @@ distributions.
NumPy is used extensively within the Python API for its powerful
N-dimensional array.
`h5py <http://www.h5py.org/>`_
h5py provides Python bindings to the HDF5 library. Since OpenMC outputs
various HDF5 files, h5py is needed to provide access to data within these
files from Python.
.. admonition:: Optional
:class: note
`SciPy <https://www.scipy.org/>`_
SciPy's special functions, sparse matrices, and spatial data structures
are used for several optional features in the API.
@ -417,6 +409,14 @@ distributions.
Pandas is used to generate tally DataFrames as demonstrated in
:ref:`examples_pandas` example notebook.
`h5py <http://www.h5py.org/>`_
h5py provides Python bindings to the HDF5 library. Since OpenMC outputs
various HDF5 files, h5py is needed to provide access to data within these
files from Python.
.. admonition:: Optional
:class: note
`Matplotlib <http://matplotlib.org/>`_
Matplotlib is used to providing plotting functionality in the API like the
:meth:`Universe.plot` method and the :func:`openmc.plot_xs` function.

View file

@ -678,7 +678,7 @@
"fission_rates.shape = mesh.dimension\n",
"\n",
"# Normalize to the average pin power\n",
"fission_rates /= np.mean(fission_rates)\n",
"fission_rates /= np.mean(fission_rates[fission_rates > 0.])\n",
"\n",
"# Force zeros to be NaNs so their values are not included when matplotlib calculates\n",
"# the color scale\n",

View file

@ -1187,7 +1187,7 @@
"mg_fission_rates.shape = (17,17)\n",
"\n",
"# Normalize to the average pin power\n",
"mg_fission_rates /= np.mean(mg_fission_rates)"
"mg_fission_rates /= np.mean(mg_fission_rates[mg_fission_rates > 0.])"
]
},
{
@ -1213,7 +1213,7 @@
"ce_fission_rates.shape = (17,17)\n",
"\n",
"# Normalize to the average pin power\n",
"ce_fission_rates /= np.mean(ce_fission_rates)"
"ce_fission_rates /= np.mean(ce_fission_rates[ce_fission_rates > 0.])"
]
},
{

View file

@ -1413,7 +1413,7 @@
" fiss_rates[-1].shape = mesh.dimension\n",
" \n",
" # Normalize the fission rates\n",
" fiss_rates[-1] /= np.mean(fiss_rates[-1])\n",
" fiss_rates[-1] /= np.mean(fiss_rates[-1][fiss_rates[-1] > 0.])\n",
" \n",
" # Set 0s to NaNs so they show as white\n",
" fiss_rates[-1][fiss_rates[-1] == 0.] = np.nan\n",

View file

@ -89,7 +89,7 @@
"\n",
"$$\\sigma_{n,x,k,g} = \\frac{\\int_{E_{g}}^{E_{g-1}}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\sigma_{n,x}(\\mathbf{r},E')\\Phi(\\mathbf{r},E')}{\\int_{E_{g}}^{E_{g-1}}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\Phi(\\mathbf{r},E')}$$\n",
"\n",
"This scalar flux-weighted average microscopic cross section is computed by `openmc.mgxs` for most multi-group cross sections, including total, absorption, and fission reaction types. These double integrals are stochastically computed with OpenMC's tally system - in particular, [filters](https://mit-crpg.github.io/openmc/pythonapi/filter.html) on the energy range and spatial zone (material, cell or universe) define the bounds of integration for both numerator and denominator."
"This scalar flux-weighted average microscopic cross section is computed by `openmc.mgxs` for most multi-group cross sections, including total, absorption, and fission reaction types. These double integrals are stochastically computed with OpenMC's tally system - in particular, [filters](http://openmc.readthedocs.io/en/latest/usersguide/tallies.html#filters) on the energy range and spatial zone (material, cell or universe) define the bounds of integration for both numerator and denominator."
]
},
{
@ -764,7 +764,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the `openmc.mgxs` module uses [tally arithmetic](https://mit-crpg.github.io/openmc/pythonapi/examples/tally-arithmetic.html) under-the-hood, the cross section is stored as a \"derived\" `Tally` object. This means that it can be queried and manipulated using all of the same methods supported for the `Tally` class in the OpenMC Python API. For example, we can construct a [Pandas](http://pandas.pydata.org/) `DataFrame` of the multi-group cross section data."
"Since the `openmc.mgxs` module uses [tally arithmetic](http://openmc.readthedocs.io/en/latest/examples/tally-arithmetic.html) under-the-hood, the cross section is stored as a \"derived\" `Tally` object. This means that it can be queried and manipulated using all of the same methods supported for the `Tally` class in the OpenMC Python API. For example, we can construct a [Pandas](http://pandas.pydata.org/) `DataFrame` of the multi-group cross section data."
]
},
{
@ -875,7 +875,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we illustrate how one can leverage OpenMC's [tally arithmetic](https://mit-crpg.github.io/openmc/pythonapi/examples/tally-arithmetic.html) data processing feature with `MGXS` objects. The `openmc.mgxs` module uses tally arithmetic to compute multi-group cross sections with automated uncertainty propagation. Each `MGXS` object includes an `xs_tally` attribute which is a \"derived\" `Tally` based on the tallies needed to compute the cross section type of interest. These derived tallies can be used in subsequent tally arithmetic operations. For example, we can use tally artithmetic to confirm that the `TotalXS` is equal to the sum of the `AbsorptionXS` and `ScatterXS` objects."
"Finally, we illustrate how one can leverage OpenMC's [tally arithmetic](http://openmc.readthedocs.io/en/latest/examples/tally-arithmetic.html) data processing feature with `MGXS` objects. The `openmc.mgxs` module uses tally arithmetic to compute multi-group cross sections with automated uncertainty propagation. Each `MGXS` object includes an `xs_tally` attribute which is a \"derived\" `Tally` based on the tallies needed to compute the cross section type of interest. These derived tallies can be used in subsequent tally arithmetic operations. For example, we can use tally artithmetic to confirm that the `TotalXS` is equal to the sum of the `AbsorptionXS` and `ScatterXS` objects."
]
},
{
@ -1176,7 +1176,7 @@
"# Use tally arithmetic to ensure that the absorption- and scattering-to-total MGXS ratios sum to unity\n",
"sum_ratio = absorption_to_total + scattering_to_total\n",
"\n",
"# The scattering-to-total ratio is a derived tally which can generate Pandas DataFrames for inspection\n",
"# The sum ratio is a derived tally which can generate Pandas DataFrames for inspection\n",
"sum_ratio.get_pandas_dataframe()"
]
}

View file

@ -8,7 +8,7 @@
"\n",
"* Creation of multi-group cross sections on a **heterogeneous geometry**\n",
"* Calculation of cross sections on a **nuclide-by-nuclide basis**\n",
"* The use of **[tally precision triggers](http://openmc.readthedocs.io/en/latest/usersguide/input.html#trigger-element)** with multi-group cross sections\n",
"* The use of **[tally precision triggers](http://openmc.readthedocs.io/en/latest/io_formats/settings.html#trigger-element)** with multi-group cross sections\n",
"* Built-in features for **energy condensation** in downstream data processing\n",
"* The use of the **`openmc.data`** module to plot continuous-energy vs. multi-group cross sections\n",
"* **Validation** of multi-group cross sections with **[OpenMOC](https://mit-crpg.github.io/OpenMOC/)**\n",
@ -254,7 +254,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we must define simulation parameters. In this case, we will use 10 inactive batches and 190 active batches each with 10,000 particles."
"Next, we must define simulation parameters. In this case, we will use 10 inactive batches and 40 active batches each with 10,000 particles."
]
},
{
@ -350,7 +350,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we showcase the use of OpenMC's [tally precision trigger](https://mit-crpg.github.io/openmc/usersguide/input.html#trigger-element) feature in conjunction with the `openmc.mgxs` module. In particular, we will assign a tally trigger of 1E-2 on the standard deviation for each of the tallies used to compute multi-group cross sections."
"Next, we showcase the use of OpenMC's [tally precision trigger](http://openmc.readthedocs.io/en/latest/io_formats/settings.html#trigger-element) feature in conjunction with the `openmc.mgxs` module. In particular, we will assign a tally trigger of 1E-2 on the standard deviation for each of the tallies used to compute multi-group cross sections."
]
},
{
@ -962,7 +962,7 @@
},
"outputs": [],
"source": [
"# Extract the 16-group transport cross section for the fuel\n",
"# Extract the 8-group transport cross section for the fuel\n",
"fine_xs = xs_library[fuel_cell.id]['transport']\n",
"\n",
"# Condense to the 2-group structure\n",
@ -973,7 +973,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Group condensation is as simple as that! We now have a new coarse 2-group `TransportXS` in addition to our original 16-group `TransportXS`. Let's inspect the 2-group `TransportXS` by printing it to the screen and extracting a Pandas `DataFrame` as we have already learned how to do."
"Group condensation is as simple as that! We now have a new coarse 2-group `TransportXS` in addition to our original 8-group `TransportXS`. Let's inspect the 2-group `TransportXS` by printing it to the screen and extracting a Pandas `DataFrame` as we have already learned how to do."
]
},
{

View file

@ -299,7 +299,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"OpenMC requires that there is a \"root\" universe. Let us create a root cell that is filled by the pin cell universe and then assign it to the root universe."
"OpenMC requires that there is a \"root\" universe. Let us create a root cell that is filled by the assembly and then assign it to the root universe."
]
},
{
@ -513,7 +513,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we will instantiate an `openmc.mgxs.Library` for the energy groups with our the fuel assembly geometry."
"Next, we will instantiate an `openmc.mgxs.Library` for the energy groups with the fuel assembly geometry."
]
},
{
@ -524,7 +524,7 @@
},
"outputs": [],
"source": [
"# Initialize an 2-group MGXS Library for OpenMOC\n",
"# Initialize a 2-group MGXS Library for OpenMOC\n",
"mgxs_lib = openmc.mgxs.Library(geometry)\n",
"mgxs_lib.energy_groups = groups"
]
@ -1575,7 +1575,7 @@
"openmc_fission_rates.shape = (17,17)\n",
"\n",
"# Normalize to the average pin power\n",
"openmc_fission_rates /= np.mean(openmc_fission_rates)"
"openmc_fission_rates /= np.mean(openmc_fission_rates[openmc_fission_rates > 0.])"
]
},
{
@ -1607,7 +1607,7 @@
"openmoc_fission_rates = np.fliplr(openmoc_fission_rates)\n",
"\n",
"# Normalize to the average pin fission rate\n",
"openmoc_fission_rates /= np.mean(openmoc_fission_rates)"
"openmoc_fission_rates /= np.mean(openmoc_fission_rates[openmoc_fission_rates > 0.])"
]
},
{

File diff suppressed because one or more lines are too long

View file

@ -240,7 +240,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"So far you've seen the \"hard\" way to create a material. The \"easy\" way is to just pass strings to `add_nuclide()` and `add_element()` -- they are implicitly coverted to `Nuclide` and `Element` objects. For example, we could have created our UO2 material as follows:"
"So far you've seen the \"hard\" way to create a material. The \"easy\" way is to just pass strings to `add_nuclide()` and `add_element()` -- they are implicitly converted to `Nuclide` and `Element` objects. For example, we could have created our UO2 material as follows:"
]
},
{
@ -509,12 +509,12 @@
"At this point, we have three materials defined, exported to XML, and ready to be used in our model. To finish our model, we need to define the geometric arrangement of materials. OpenMC represents physical volumes using constructive solid geometry (CSG), also known as combinatorial geometry. The object that allows us to assign a material to a region of space is called a `Cell` (same concept in MCNP, for those familiar). In order to define a region that we can assign to a cell, we must first define surfaces which bound the region. A *surface* is a locus of zeros of a function of Cartesian coordinates $x$, $y$, and $z$, e.g.\n",
"\n",
"- A plane perpendicular to the x axis: $x - x_0 = 0$\n",
"- A cylinder perpendicular to the z axis: $(x - x_0)^2 + (y - y_0)^2 - R^2 = 0$\n",
"- A cylinder parallel to the z axis: $(x - x_0)^2 + (y - y_0)^2 - R^2 = 0$\n",
"- A sphere: $(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0$\n",
"\n",
"Between those three classes of surfaces (planes, cylinders, spheres), one can construct a wide variety of models. It is also possible to define cones and general second-order surfaces (torii are not currently supported).\n",
"Between those three classes of surfaces (planes, cylinders, spheres), one can construct a wide variety of models. It is also possible to define cones and general second-order surfaces (tori are not currently supported).\n",
"\n",
"Note that defining a surface is not sufficient to specify a volume -- in order to define an actual volume, one must reference the half-space of a surface. A surface *half-space* is the region whose points satisfy a positive of negative inequality of the surface equation. For example, for a sphere of radius one centered at the origin, the surface equation is $f(x,y,z) = x^2 + y^2 + z^2 - 1 = 0$. Thus, we say that the negative half-space of the sphere, is defined as the collection of points satisfying $f(x,y,z) < 0$, which one can reason is the inside of the sphere. Conversely, the positive half-space of the sphere would correspond to all points outside of the sphere.\n",
"Note that defining a surface is not sufficient to specify a volume -- in order to define an actual volume, one must reference the half-space of a surface. A surface *half-space* is the region whose points satisfy a positive or negative inequality of the surface equation. For example, for a sphere of radius one centered at the origin, the surface equation is $f(x,y,z) = x^2 + y^2 + z^2 - 1 = 0$. Thus, we say that the negative half-space of the sphere, is defined as the collection of points satisfying $f(x,y,z) < 0$, which one can reason is the inside of the sphere. Conversely, the positive half-space of the sphere would correspond to all points outside of the sphere.\n",
"\n",
"Let's go ahead and create a sphere and confirm that what we've told you is true."
]
@ -792,7 +792,7 @@
"\n",
"We now have enough knowledge to create our pin-cell. We need three surfaces to define the fuel and clad:\n",
"\n",
"1. The outer surface of the fuel -- a cylinder perpendicular to the z axis\n",
"1. The outer surface of the fuel -- a cylinder parallel to the z axis\n",
"2. The inner surface of the clad -- same as above\n",
"3. The outer surface of the clad -- same as above\n",
"\n",

View file

@ -131,7 +131,6 @@
"clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720)\n",
"\n",
"# Create boundary planes to surround the geometry\n",
"# Use both reflective and vacuum boundaries to make life interesting\n",
"min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective')\n",
"max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')\n",
"min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')\n",
@ -292,7 +291,7 @@
"plot.origin = [0, 0, 0]\n",
"plot.width = [1.26, 1.26]\n",
"plot.pixels = [250, 250]\n",
"plot.color = 'mat'\n",
"plot.color_by = 'material'\n",
"\n",
"# Instantiate a Plots collection and export to \"plots.xml\"\n",
"plot_file = openmc.Plots([plot])\n",
@ -686,8 +685,7 @@
"Tally\n",
"\tID =\t10000\n",
"\tName =\tflux\n",
"\tFilters =\t\n",
" \t\tMeshFilter\t[10000]\n",
"\tFilters =\tMeshFilter\t[10000]\n",
"\tNuclides =\ttotal \n",
"\tScores =\t['flux', 'fission']\n",
"\tEstimator =\ttracklength\n",
@ -704,7 +702,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"The statepoint file actually stores the sum and sum-of-squares for each tally bin from which the mean and variance can be calculated as described [here](http://mit-crpg.github.io/openmc/methods/tallies.html#variance). The sum and sum-of-squares can be accessed using the ``sum`` and ``sum_sq`` properties:"
"The statepoint file actually stores the sum and sum-of-squares for each tally bin from which the mean and variance can be calculated as described [here](http://openmc.readthedocs.io/en/latest/methods/tallies.html#variance). The sum and sum-of-squares can be accessed using the ``sum`` and ``sum_sq`` properties:"
]
},
{
@ -821,8 +819,7 @@
"Tally\n",
"\tID =\t10001\n",
"\tName =\tflux\n",
"\tFilters =\t\n",
" \t\tMeshFilter\t[10000]\n",
"\tFilters =\tMeshFilter\t[10000]\n",
"\tNuclides =\ttotal \n",
"\tScores =\t['flux']\n",
"\tEstimator =\ttracklength\n",
@ -1043,7 +1040,7 @@
}
],
"source": [
"# Create log-spaced energy bins from 1 keV to 100 MeV\n",
"# Create log-spaced energy bins from 1 keV to 10 MeV\n",
"energy_bins = np.logspace(3,7)\n",
"\n",
"# Calculate pdf for source energies\n",

View file

@ -37,7 +37,7 @@
"\n",
"To perform the search we will use the `openmc.search_for_keff` function. This function requires a different function be defined which creates an parametrized model to analyze. This model is required to be stored in an `openmc.model.Model` object. The first parameter of this function will be modified during the search process for our critical eigenvalue.\n",
"\n",
"Our model will be a pin-cell from the [Multi-Group Mode Part II](./mg-mode-part-ii.rst) assembly, except this time the entire model building process will be contained within a function, and the Boron concentration will be parametrized."
"Our model will be a pin-cell from the [Multi-Group Mode Part II](http://openmc.readthedocs.io/en/latest/examples/mg-mode-part-ii.html) assembly, except this time the entire model building process will be contained within a function, and the Boron concentration will be parametrized."
]
},
{

View file

@ -115,7 +115,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's move on to the geometry. Our problem will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces -- in this case two cylinders and six reflective planes."
"Now let's move on to the geometry. Our problem will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces -- in this case two cylinders and six planes."
]
},
{
@ -293,7 +293,7 @@
"plot.origin = [0, 0, 0]\n",
"plot.width = [1.26, 1.26]\n",
"plot.pixels = [250, 250]\n",
"plot.color = 'mat'\n",
"plot.color_by = 'material'\n",
"\n",
"# Instantiate a Plots collection and export to \"plots.xml\"\n",
"plot_file = openmc.Plots([plot])\n",
@ -412,7 +412,7 @@
"tally.nuclides = [o16, h1]\n",
"tallies_file.append(tally)\n",
"\n",
"# Instantiate a tally mesh \n",
"# Instantiate a tally mesh\n",
"mesh = openmc.Mesh(mesh_id=1)\n",
"mesh.type = 'regular'\n",
"mesh.dimension = [1, 1, 1]\n",

View file

@ -1,21 +1,30 @@
<?xml version="1.0"?>
<tallies>
<filter id="1" type="cell">
<bins>100</bins>
</filter>
<filter id="2" type="energy">
<bins>0 20.0e6</bins>
</filter>
<filter id="3" type="energyout">
<bins>0 20.0e6</bins>
</filter>
<tally id="1">
<filter type="cell" bins="100" />
<filters>1</filters>
<scores>total scatter nu-scatter absorption fission nu-fission</scores>
</tally>
<tally id="2">
<filter type="cell" bins="100" />
<filter type="energy" bins="0 20.0e6" />
<filters>1 2</filters>
<scores>total scatter nu-scatter absorption fission nu-fission</scores>
</tally>
<tally id="3">
<filter type="cell" bins="100" />
<filter type="energy" bins="0 20.0e6" />
<filter type="energyout" bins="0 20.0e6" />
<filters>1 2 3</filters>
<scores>scatter nu-scatter nu-fission</scores>
</tally>

View file

@ -8,8 +8,12 @@
<width>1.0 1.0</width>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1">
<filter type="mesh" bins="1" />
<filters>1</filters>
<scores>total</scores>
</tally>

View file

@ -8,8 +8,12 @@
<width>1.0 1.0</width>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1">
<filter type="mesh" bins="1" />
<filters>1</filters>
<scores>total</scores>
</tally>

View file

@ -7,9 +7,16 @@
<upper_right>0.62992 0.62992 1.e50</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>0. 4. 20.0e6</bins>
</filter>
<tally id="1">
<filter type="mesh" bins="1" />
<filter type="energy" bins="0. 4. 20.0e6" />
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
</tally>

View file

@ -5,9 +5,14 @@
<lower_left>-0.63 -0.63 -1e+50</lower_left>
<upper_right>0.63 0.63 1e+50</upper_right>
</mesh>
<filter id="1" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
</filter>
<filter id="2" type="mesh">
<bins>1</bins>
</filter>
<tally id="1" name="tally 1">
<filter bins="1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0" type="energy" />
<filter bins="1" type="mesh" />
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
</tally>
</tallies>

View file

@ -1,7 +1,6 @@
from openmc.arithmetic import *
from openmc.cell import *
from openmc.mesh import *
from openmc.lattice import *
from openmc.element import *
from openmc.geometry import *
from openmc.nuclide import *
@ -14,6 +13,7 @@ from openmc.source import *
from openmc.settings import *
from openmc.surface import *
from openmc.universe import *
from openmc.lattice import *
from openmc.filter import *
from openmc.trigger import *
from openmc.tally_derivative import *

View file

@ -4,6 +4,7 @@ from collections import Iterable
from six import string_types
import numpy as np
import pandas as pd
import openmc
from openmc.filter import _FILTER_TYPES
@ -786,9 +787,6 @@ class AggregateFilter(object):
CrossFilter.get_pandas_dataframe()
"""
import pandas as pd
# Create NumPy array of the bin tuples for repeating / tiling
filter_bins = np.empty(self.num_bins, dtype=tuple)
for i, bin in enumerate(self.bins):

38
openmc/capi/__init__.py Normal file
View file

@ -0,0 +1,38 @@
"""
This module provides bindings to C functions defined by OpenMC shared library.
When the :mod:`openmc` package is imported, the OpenMC shared library is
automatically loaded. Calls to the OpenMC library can then be via functions or
objects in the :mod:`openmc.capi` subpackage, for example:
.. code-block:: python
openmc.capi.init()
openmc.capi.run()
openmc.capi.finalize()
"""
from ctypes import CDLL
import sys
from warnings import warn
import pkg_resources
# Determine shared-library suffix
if sys.platform == 'darwin':
_suffix = 'dylib'
else:
_suffix = 'so'
# Open shared library
_filename = pkg_resources.resource_filename(
__name__, '_libopenmc.{}'.format(_suffix))
_dll = CDLL(_filename)
from .error import *
from .core import *
from .nuclide import *
from .material import *
from .cell import *
from .tally import *

87
openmc/capi/cell.py Normal file
View file

@ -0,0 +1,87 @@
from collections import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
import numpy as np
from . import _dll
from .error import _error_handler
__all__ = ['CellView', 'cells']
# Cell functions
_dll.openmc_cell_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_cell_get_id.restype = c_int
_dll.openmc_cell_get_id.errcheck = _error_handler
_dll.openmc_cell_set_temperature.argtypes = [
c_int32, c_double, POINTER(c_int32)]
_dll.openmc_cell_set_temperature.restype = c_int
_dll.openmc_cell_set_temperature.errcheck = _error_handler
_dll.openmc_get_cell.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_cell.restype = c_int
_dll.openmc_get_cell.errcheck = _error_handler
class CellView(object):
"""View of a cell.
This class exposes a cell that is stored internally in the OpenMC solver. To
obtain a view of a cell with a given ID, use the
:data:`openmc.capi.nuclides` mapping.
Parameters
----------
index : int
Index in the `cells` array.
Attributes
----------
id : int
ID of the cell
"""
__instances = WeakValueDictionary()
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __init__(self, index):
self._index = index
@property
def id(self):
cell_id = c_int32()
_dll.openmc_cell_get_id(self._index, cell_id)
return cell_id.value
def set_temperature(self, T, instance=None):
"""Set the temperature of a cell
Parameters
----------
T : float
Temperature in K
instance : int or None
Which instance of the cell
"""
_dll.openmc_cell_set_temperature(self._index, T, instance)
class _CellMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_cell(key, index)
return CellView(index.value)
def __iter__(self):
for i in range(len(self)):
yield CellView(i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_cells').value
cells = _CellMapping()

157
openmc/capi/core.py Normal file
View file

@ -0,0 +1,157 @@
from contextlib import contextmanager
from ctypes import CDLL, c_int, c_int32, c_double, POINTER
from warnings import warn
from . import _dll
from .error import _error_handler
_dll.openmc_calculate_volumes.restype = None
_dll.openmc_finalize.restype = None
_dll.openmc_find.argtypes = [POINTER(c_double*3), c_int, POINTER(c_int32),
POINTER(c_int32)]
_dll.openmc_find.restype = c_int
_dll.openmc_find.errcheck = _error_handler
_dll.openmc_hard_reset.restype = None
_dll.openmc_init.argtypes = [POINTER(c_int)]
_dll.openmc_init.restype = None
_dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
_dll.openmc_get_keff.restype = c_int
_dll.openmc_get_keff.errcheck = _error_handler
_dll.openmc_plot_geometry.restype = None
_dll.openmc_run.restype = None
_dll.openmc_reset.restype = None
def calculate_volumes():
"""Run stochastic volume calculation"""
_dll.openmc_calculate_volumes()
def finalize():
"""Finalize simulation and free memory"""
_dll.openmc_finalize()
def find_cell(xyz):
"""Find the cell at a given point
Parameters
----------
xyz : iterable of float
Cartesian coordinates of position
Returns
-------
int
ID of the cell.
int
If the cell at the given point is repeated in the geometry, this
indicates which instance it is, i.e., 0 would be the first instance.
"""
uid = c_int32()
instance = c_int32()
_dll.openmc_find((c_double*3)(*xyz), 1, uid, instance)
return uid.value, instance.value
def find_material(xyz):
"""Find the material at a given point
Parameters
----------
xyz : iterable of float
Cartesian coordinates of position
Returns
-------
int or None
ID of the material or None is no material is found
"""
uid = c_int32()
instance = c_int32()
_dll.openmc_find((c_double*3)(*xyz), 2, uid, instance)
return uid.value if uid != 0 else None
def hard_reset():
"""Reset tallies, timers, and pseudo-random number generator state."""
_dll.openmc_hard_reset()
def init(intracomm=None):
"""Initialize OpenMC
Parameters
----------
intracomm : mpi4py.MPI.Intracomm or None
MPI intracommunicator
"""
if intracomm is not None:
# If an mpi4py communicator was passed, convert it to an integer to
# be passed to openmc_init
try:
intracomm = intracomm.py2f()
except AttributeError:
pass
_dll.openmc_init(c_int(intracomm))
else:
_dll.openmc_init(None)
def keff():
"""Return the calculated k-eigenvalue and its standard deviation.
Returns
-------
tuple
Mean k-eigenvalue and standard deviation of the mean
"""
k = (c_double*2)()
_dll.openmc_get_keff(k)
return tuple(k)
def plot_geometry():
"""Plot geometry"""
_dll.openmc_plot_geometry()
def reset():
"""Reset tallies and timers."""
_dll.openmc_reset()
def run():
"""Run simulation"""
_dll.openmc_run()
@contextmanager
def run_in_memory(intracomm=None):
"""Provides context manager for calling OpenMC shared library functions.
This function is intended to be used in a 'with' statement and ensures that
OpenMC is properly initialized/finalized. At the completion of the 'with'
block, all memory that was allocated during the block is freed. For
example::
with openmc.capi.run_in_memory():
for i in range(n_iters):
openmc.capi.reset()
do_stuff()
openmc.capi.run()
Parameters
----------
intracomm : mpi4py.MPI.Intracomm or None
MPI intracommunicator
"""
init(intracomm)
yield
finalize()

66
openmc/capi/error.py Normal file
View file

@ -0,0 +1,66 @@
from ctypes import c_int
from . import _dll
class GeometryError(Exception):
pass
def _error_handler(err, func, args):
"""Raise exception according to error code."""
# Get error code corresponding to global constant.
def errcode(s):
return c_int.in_dll(_dll, s).value
if err == errcode('e_out_of_bounds'):
raise IndexError('Array index out of bounds.')
elif err == errcode('e_cell_not_allocated'):
raise MemoryError("Memory has not been allocated for cells.")
elif err == errcode('e_cell_invalid_id'):
raise KeyError("No cell exists with ID={}.".format(args[0]))
elif err == errcode('e_cell_not_found'):
raise GeometryError("Could not find cell at position ({}, {}, {})"
.format(*args[0]))
elif err == errcode('e_nuclide_not_allocated'):
raise MemoryError("Memory has not been allocated for nuclides.")
elif err == errcode('e_nuclide_not_loaded'):
raise KeyError("No nuclide named '{}' has been loaded.")
elif err == errcode('e_nuclide_not_in_library'):
raise KeyError("Specified nuclide doesn't exist in the cross "
"section library.")
elif err == errcode('e_material_not_allocated'):
raise MemoryError("Memory has not been allocated for materials.")
elif err == errcode('e_material_invalid_id'):
raise KeyError("No material exists with ID={}.".format(args[0]))
elif err == errcode('e_tally_not_allocated'):
raise MemoryError("Memory has not been allocated for tallies.")
elif err == errcode('e_tally_invalid_id'):
raise KeyError("No tally exists with ID={}.".format(args[0]))
elif err == errcode('e_invalid_size'):
raise MemoryError("Array size mismatch with memory allocated.")
elif err == errcode('e_cell_no_material'):
raise GeometryError("Operation on cell requires that it be filled"
" with a material.")
elif err == errcode('w_below_min_bound'):
warn("Data has not been loaded beyond lower bound of {}.".format(args[0]))
elif err == errcode('w_above_max_bound'):
warn("Data has not been loaded beyond upper bound of {}.".format(args[0]))
elif err < 0:
raise Exception("Unknown error encountered (code {}).".format(err))

170
openmc/capi/material.py Normal file
View file

@ -0,0 +1,170 @@
from collections import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from . import _dll, NuclideView
from .error import _error_handler
__all__ = ['MaterialView', 'materials']
# Material functions
_dll.openmc_get_material.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_material.restype = c_int
_dll.openmc_get_material.errcheck = _error_handler
_dll.openmc_material_add_nuclide.argtypes = [
c_int32, c_char_p, c_double]
_dll.openmc_material_add_nuclide.restype = c_int
_dll.openmc_material_add_nuclide.errcheck = _error_handler
_dll.openmc_material_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_material_get_id.restype = c_int
_dll.openmc_material_get_id.errcheck = _error_handler
_dll.openmc_material_get_densities.argtypes = [
c_int32, POINTER(POINTER(c_int)), POINTER(POINTER(c_double)),
POINTER(c_int)]
_dll.openmc_material_get_densities.restype = c_int
_dll.openmc_material_get_densities.errcheck = _error_handler
_dll.openmc_material_set_density.argtypes = [c_int32, c_double]
_dll.openmc_material_set_density.restype = c_int
_dll.openmc_material_set_density.errcheck = _error_handler
_dll.openmc_material_set_densities.argtypes = [
c_int32, c_int, POINTER(c_char_p), POINTER(c_double)]
_dll.openmc_material_set_densities.restype = c_int
_dll.openmc_material_set_densities.errcheck = _error_handler
class MaterialView(object):
"""View of a material.
This class exposes a material that is stored internally in the OpenMC
solver. To obtain a view of a material with a given ID, use the
:data:`openmc.capi.materials` mapping.
Parameters
----------
index : int
Index in the `materials` array.
Attributes
----------
id : int
ID of the material
nuclides : list of str
List of nuclides in the material
densities : numpy.ndarray
Array of densities in atom/b-cm
"""
__instances = WeakValueDictionary()
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __init__(self, index):
self._index = index
@property
def id(self):
mat_id = c_int32()
_dll.openmc_material_get_id(self._index, mat_id)
return mat_id.value
@property
def nuclides(self):
return self._get_densities()[0]
return nuclides
@property
def densities(self):
return self._get_densities()[1]
def _get_densities(self):
"""Get atom densities in a material.
Returns
-------
list of string
List of nuclide names
numpy.ndarray
Array of densities in atom/b-cm
"""
# Allocate memory for arguments that are written to
nuclides = POINTER(c_int)()
densities = POINTER(c_double)()
n = c_int()
# Get nuclide names and densities
_dll.openmc_material_get_densities(self._index, nuclides, densities, n)
# Convert to appropriate types and return
nuclide_list = [NuclideView(nuclides[i]).name for i in range(n.value)]
density_array = as_array(densities, (n.value,))
return nuclide_list, density_array
def add_nuclide(name, density):
"""Add a nuclide to a material.
Parameters
----------
name : str
Name of nuclide, e.g. 'U235'
density : float
Density in atom/b-cm
"""
_dll.openmc_material_add_nuclide(self._index, name.encode(), density)
def set_density(self, density):
"""Set density of a material.
Parameters
----------
density : float
Density in atom/b-cm
"""
_dll.openmc_material_set_density(self._index, density)
def set_densities(self, nuclides, densities):
"""Set the densities of a list of nuclides in a material
Parameters
----------
nuclides : iterable of str
Nuclide names
densities : iterable of float
Corresponding densities in atom/b-cm
"""
# Convert strings to an array of char*
nucs = (c_char_p * len(nuclides))()
nucs[:] = [x.encode() for x in nuclides]
# Get numpy array as a double*
d = np.asarray(densities)
dp = d.ctypes.data_as(POINTER(c_double))
_dll.openmc_material_set_densities(self._index, len(nuclides), nucs, dp)
class _MaterialMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_material(key, index)
return MaterialView(index.value)
def __iter__(self):
for i in range(len(self)):
yield MaterialView(i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_materials').value
materials = _MaterialMapping()

93
openmc/capi/nuclide.py Normal file
View file

@ -0,0 +1,93 @@
from collections import Mapping
from ctypes import c_int, c_char_p, POINTER
from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from . import _dll
from .error import _error_handler
__all__ = ['NuclideView', 'nuclides', 'load_nuclide']
# Nuclide functions
_dll.openmc_get_nuclide.argtypes = [c_char_p, POINTER(c_int)]
_dll.openmc_get_nuclide.restype = c_int
_dll.openmc_get_nuclide.errcheck = _error_handler
_dll.openmc_load_nuclide.argtypes = [c_char_p]
_dll.openmc_load_nuclide.restype = c_int
_dll.openmc_load_nuclide.errcheck = _error_handler
_dll.openmc_nuclide_name.argtypes = [c_int, POINTER(c_char_p)]
_dll.openmc_nuclide_name.restype = c_int
_dll.openmc_nuclide_name.errcheck = _error_handler
def load_nuclide(name):
"""Load cross section data for a nuclide.
Parameters
----------
name : str
Name of the nuclide, e.g. 'U235'
"""
_dll.openmc_load_nuclide(name.encode())
class NuclideView(object):
"""View of a nuclide.
This class exposes a nuclide that is stored internally in the OpenMC
solver. To obtain a view of a nuclide with a given name, use the
:data:`openmc.capi.nuclides` mapping.
Parameters
----------
index : int
Index in the `nuclides` array.
Attributes
----------
name : str
Name of the nuclide, e.g. 'U235'
"""
__instances = WeakValueDictionary()
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __init__(self, index):
self._index = index
@property
def name(self):
name = c_char_p()
_dll.openmc_nuclide_name(self._index, name)
# Find blank in name
i = 0
while name.value[i:i+1] != b' ':
i += 1
return name.value[:i].decode()
class _NuclideMapping(Mapping):
"""Provide mapping from nuclide name to index in nuclides array."""
def __getitem__(self, key):
index = c_int()
_dll.openmc_get_nuclide(key.encode(), index)
return NuclideView(index)
def __iter__(self):
for i in range(len(self)):
yield NuclideView(i + 1).name
def __len__(self):
return c_int.in_dll(_dll, 'n_nuclides').value
nuclides = _NuclideMapping()

107
openmc/capi/tally.py Normal file
View file

@ -0,0 +1,107 @@
from collections import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
from numpy.ctypeslib import as_array
from . import _dll, NuclideView
from .error import _error_handler
__all__ = ['TallyView', 'tallies']
# Tally functions
_dll.openmc_get_tally.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_tally.restype = c_int
_dll.openmc_get_tally.errcheck = _error_handler
_dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_id.restype = c_int
_dll.openmc_tally_get_id.errcheck = _error_handler
_dll.openmc_tally_get_nuclides.argtypes = [
c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_tally_get_nuclides.restype = c_int
_dll.openmc_tally_get_nuclides.errcheck = _error_handler
_dll.openmc_tally_results.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)]
_dll.openmc_tally_results.restype = c_int
_dll.openmc_tally_results.errcheck = _error_handler
_dll.openmc_tally_set_nuclides.argtypes = [c_int32, c_int, POINTER(c_char_p)]
_dll.openmc_tally_set_nuclides.restype = c_int
_dll.openmc_tally_set_nuclides.errcheck = _error_handler
class TallyView(object):
"""View of a tally.
This class exposes a tally that is stored internally in the OpenMC
solver. To obtain a view of a tally with a given ID, use the
:data:`openmc.capi.tallies` mapping.
Parameters
----------
index : int
Index in the `tallys` array.
Attributes
----------
id : int
ID of the tally
nuclides : list of str
List of nuclides to score results for
results : numpy.ndarray
Array of tally results
"""
__instances = WeakValueDictionary()
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __init__(self, index):
self._index = index
@property
def id(self):
tally_id = c_int32()
_dll.openmc_tally_get_id(self._index, tally_id)
return tally_id.value
@property
def nuclides(self):
nucs = POINTER(c_int)()
n = c_int()
_dll.openmc_tally_get_nuclides(self._index, nucs, n)
return [NuclideView(nucs[i]).name if nucs[i] > 0 else 'total'
for i in range(n.value)]
@property
def results(self):
data = POINTER(c_double)()
shape = (c_int*3)()
_dll.openmc_tally_results(self._index, data, shape)
return as_array(data, tuple(shape[::-1]))
@nuclides.setter
def nuclides(self, nuclides):
nucs = (c_char_p * len(nuclides))()
nucs[:] = [x.encode() for x in nuclides]
_dll.openmc_tally_set_nuclides(self._index, len(nuclides), nucs)
class _TallyMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_tally(key, index)
return TallyView(index.value)
def __iter__(self):
for i in range(len(self)):
yield TallyView(i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_tallies').value
tallies = _TallyMapping()

View file

@ -13,19 +13,10 @@ import openmc
import openmc.checkvalue as cv
from openmc.surface import Halfspace
from openmc.region import Region, Intersection, Complement
from .mixin import IDManagerMixin
# A static variable for auto-generated Cell IDs
AUTO_CELL_ID = 10000
def reset_auto_cell_id():
"""Reset counter for auto-generated cell IDs."""
global AUTO_CELL_ID
AUTO_CELL_ID = 10000
class Cell(object):
class Cell(IDManagerMixin):
r"""A region of space defined as the intersection of half-space created by
quadric surfaces.
@ -93,6 +84,9 @@ class Cell(object):
"""
next_id = 1
used_ids = set()
def __init__(self, cell_id=None, name='', fill=None, region=None):
# Initialize Cell class attributes
self.id = cell_id
@ -103,7 +97,8 @@ class Cell(object):
self._rotation_matrix = None
self._temperature = None
self._translation = None
self._paths = []
self._paths = None
self._num_instances = None
self._volume = None
self._atoms = None
@ -163,10 +158,6 @@ class Cell(object):
return string
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -214,7 +205,7 @@ class Cell(object):
@property
def paths(self):
if not self._paths:
if self._paths is None:
raise ValueError('Cell instance paths have not been determined. '
'Call the Geometry.determine_paths() method.')
return self._paths
@ -229,18 +220,11 @@ class Cell(object):
@property
def num_instances(self):
return len(self.paths)
@id.setter
def id(self, cell_id):
if cell_id is None:
global AUTO_CELL_ID
self._id = AUTO_CELL_ID
AUTO_CELL_ID += 1
else:
cv.check_type('cell ID', cell_id, Integral)
cv.check_greater_than('cell ID', cell_id, 0, equality=True)
self._id = cell_id
if self._num_instances is None:
raise ValueError(
'Number of cell instances have not been determined. Call the '
'Geometry.determine_paths() method.')
return self._num_instances
@name.setter
def name(self, name):
@ -373,7 +357,7 @@ class Cell(object):
self.region = region
else:
if isinstance(self.region, Intersection):
self.region.nodes.append(region)
self.region &= region
else:
self.region = Intersection(self.region, region)
@ -527,8 +511,17 @@ class Cell(object):
# If no nemoize'd clone exists, instantiate one
if self not in memo:
# Temporarily remove paths
paths = self._paths
self._paths = None
clone = deepcopy(self)
clone.id = None
clone._num_instances = None
# Restore paths on original instance
self._paths = paths
if self.region is not None:
clone.region = self.region.clone(memo)
if self.fill is not None:
@ -542,7 +535,7 @@ class Cell(object):
memo[self] = clone
return memo[self]
def create_xml_subelement(self, xml_element):
element = ET.Element("cell")
element.set("id", str(self.id))
@ -586,7 +579,7 @@ class Cell(object):
elif isinstance(node, Complement):
create_surface_elements(node.node, element)
else:
for subnode in node.nodes:
for subnode in node:
create_surface_elements(subnode, element)
# Call the recursive function from the top node

View file

@ -45,8 +45,8 @@ def sort_xml_elements(tree):
key = element.get('id')
# If this element has an "ID" tag, append it to the list to sort
if not key is None:
tagged_data.append((key, element))
if key is not None:
tagged_data.append((int(key), element))
# Sort the elements according to the IDs for this tag
tagged_data.sort()
@ -65,7 +65,7 @@ def sort_xml_elements(tree):
tree.extend(sorted_elements)
def clean_xml_indentation(element, level=0, spaces_per_level=4):
def clean_xml_indentation(element, level=0, spaces_per_level=2):
"""
copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
it basically walks your tree and adds spaces and newlines so the tree is

View file

@ -353,8 +353,9 @@ class Evaluation(object):
# Information about target/projectile
items = get_head_record(file_obj)
self.target['atomic_number'] = items[0] // 1000
self.target['mass_number'] = items[0] % 1000
Z, A = divmod(items[0], 1000)
self.target['atomic_number'] = Z
self.target['mass_number'] = A
self.target['mass'] = items[1]
self._LRP = items[2]
self.target['fissionable'] = (items[3] == 1)
@ -370,10 +371,14 @@ class Evaluation(object):
self.target['excitation_energy'] = items[0]
self.target['stable'] = (int(items[1]) == 0)
self.target['state'] = items[2]
self.target['isomeric_state'] = items[3]
self.target['isomeric_state'] = m = items[3]
self.info['format'] = items[5]
assert self.info['format'] == 6
# Set correct excited state for Am242_m1, which is wrong in ENDF/B-VII.1
if Z == 95 and A == 242 and m == 1:
self.target['state'] = 2
# Control record 2
items = get_cont_record(file_obj)
self.projectile['mass'] = items[0]

View file

@ -866,12 +866,13 @@ class IncidentNeutron(EqualityMixin):
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)
# Add 0K elastic scattering cross section
pendf = Evaluation(pendf_file)
file_obj = StringIO(pendf.section[3, 2])
get_head_record(file_obj)
params, xs = get_tab1_record(file_obj)
data.energy['0K'] = xs.x
data[2].xs['0K'] = xs
if '0K' not in data.energy:
pendf = Evaluation(pendf_file)
file_obj = StringIO(pendf.section[3, 2])
get_head_record(file_obj)
params, xs = get_tab1_record(file_obj)
data.energy['0K'] = xs.x
data[2].xs['0K'] = xs
finally:
# Get rid of temporary files

View file

@ -50,48 +50,45 @@ _THERMAL_DATA = {
75: ThermalTuple('ouo2', [8016, 8017, 8018], 1),
}
_PENDF_TEMPLATE = """
_TEMPLATE_RECONR = """
reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%%
20 22
{nendf} {npendf}
'{library} PENDF for {zsymam}'/
{mat} 2/
0.001 0.0 0.003/ err tempr errmax
{error}/ err
'{library}: {zsymam}'/
'Processed by NJOY'/
0/
stop
"""
_ACE_TEMPLATE = """
reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%%
20 21
'{library} PENDF for {zsymam}'/
{mat} 2/
0.001 0.0 0.003/ err tempr errmax
'{library}: {zsymam}'/
'Processed by NJOY'/
0/
_TEMPLATE_BROADR = """
broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
20 21 22
{nendf} {npendf} {nbroadr}
{mat} {num_temp} 0 0 0. /
0.001 1.0e6 0.003 /
{error}/ errthn
{temps}
0/
"""
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
20 22 23 /
{nendf} {nheatr_in} {nheatr} /
{mat} 3 /
302 318 402 /
"""
_TEMPLATE_PURR = """
purr / %%%%%%%%%%%%%%%%%%%%%%%% Add probability tables %%%%%%%%%%%%%%%%%%%%%%%%%
20 23 24
{nendf} {npurr_in} {npurr} /
{mat} {num_temp} 1 20 64 /
{temps}
1.e10
0/
"""
_ACE_TEMPLATE_ACER = """acer /
20 24 0 {nace} {ndir}
_TEMPLATE_ACER = """
acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nacer_in} 0 {nace} {ndir}
1 0 1 .{ext} /
'{library}: {zsymam} at {temperature}'/
{mat} {temperature}
@ -99,35 +96,22 @@ _ACE_TEMPLATE_ACER = """acer /
/
"""
_ACE_THERMAL_TEMPLATE = """
reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%%
20 22
'{library} PENDF for {zsymam}'/
{mat} 2/
0.001 0. 0.001/ err tempr errmax
'{library}: PENDF for {zsymam}'/
'Processed by NJOY'/
0/
broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
20 22 23
{mat} {num_temp} 0 0 0./
0.001 2.0e+6 0.001/ errthn thnmax errmax
{temps}
0/
_THERMAL_TEMPLATE_THERMR = """
thermr / %%%%%%%%%%%%%%%% Add thermal scattering data (free gas) %%%%%%%%%%%%%%%
0 23 62
0 {nthermr1_in} {nthermr1}
0 {mat} 12 {num_temp} 1 0 {iform} 1 221 1/
{temps}
0.001 4.0
{error} {energy_max}
thermr / %%%%%%%%%%%%%%%% Add thermal scattering data (bound) %%%%%%%%%%%%%%%%%%
60 62 27
{nthermal_endf} {nthermr2_in} {nthermr2}
{mat_thermal} {mat} 16 {num_temp} {inelastic} {elastic} {iform} {natom} 222 1/
{temps}
0.001 4.0
{error} {energy_max}
"""
_ACE_THERMAL_TEMPLATE_ACER = """acer /
20 27 0 {nace} {ndir}
_THERMAL_TEMPLATE_ACER = """
acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nthermal_acer_in} 0 {nace} {ndir}
2 0 1 .{ext}/
'{library}: {zsymam_thermal} processed by NJOY'/
{mat} {temperature} '{data.name}' /
@ -136,7 +120,8 @@ _ACE_THERMAL_TEMPLATE_ACER = """acer /
"""
def run(commands, tapein, tapeout, stdout=False, njoy_exec='njoy'):
def run(commands, tapein, tapeout, input_filename=None, stdout=False,
njoy_exec='njoy'):
"""Run NJOY with given commands
Parameters
@ -147,6 +132,8 @@ def run(commands, tapein, tapeout, stdout=False, njoy_exec='njoy'):
Dictionary mapping tape numbers to paths for any input files
tapeout : dict
Dictionary mapping tape numbers to paths for any output files
input_filename : str, optional
File name to write out NJOY input commands
stdout : bool, optional
Whether to display output when running NJOY
njoy_exec : str, optional
@ -159,6 +146,10 @@ def run(commands, tapein, tapeout, stdout=False, njoy_exec='njoy'):
"""
if input_filename is not None:
with open(input_filename, 'w') as f:
f.write(commands)
# Create temporary directory -- it would be preferable to use
# TemporaryDirectory(), but it is only available in Python 3.2
tmpdir = tempfile.mkdtemp()
@ -195,7 +186,7 @@ def run(commands, tapein, tapeout, stdout=False, njoy_exec='njoy'):
return njoy.returncode
def make_pendf(filename, pendf='pendf', stdout=False):
def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
"""Generate ACE file from an ENDF file
Parameters
@ -204,6 +195,8 @@ def make_pendf(filename, pendf='pendf', stdout=False):
Path to ENDF file
pendf : str, optional
Path of pointwise ENDF file to write
error : float, optional
Fractional error tolerance for NJOY processing
stdout : bool
Whether to display NJOY standard output
@ -213,21 +206,14 @@ def make_pendf(filename, pendf='pendf', stdout=False):
Return code of NJOY process
"""
ev = endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
# Determine name of library
library = '{}-{}.{}'.format(*ev.info['library'])
commands = _PENDF_TEMPLATE.format(**locals())
tapein = {20: filename}
tapeout = {22: pendf}
return run(commands, tapein, tapeout, stdout)
return make_ace(filename, pendf=pendf, error=error, broadr=False,
heatr=False, purr=False, acer=False, stdout=stdout)
def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir',
pendf=None, **kwargs):
def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
error=0.001, broadr=True, heatr=True, purr=True, acer=True,
**kwargs):
"""Generate incident neutron ACE file from an ENDF file
Parameters
@ -243,6 +229,16 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir',
Path of xsdir file to write
pendf : str, optional
Path of pendf file to write. If omitted, the pendf file is not saved.
error : float, optional
Fractional error tolerance for NJOY processing
broadr : bool, optional
Indicating whether to Doppler broaden XS when running NJOY
heatr : bool, optional
Indicating whether to add heating kerma 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
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -264,26 +260,59 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir',
num_temp = len(temperatures)
temps = ' '.join(str(i) for i in temperatures)
commands = _ACE_TEMPLATE.format(**locals())
tapein = {20: filename}
# Create njoy commands by modules
commands = ""
nendf, npendf = 20, 21
tapein = {nendf: filename}
tapeout = {}
if pendf is not None:
tapeout[21] = pendf
fname = '{}_{:.1f}'
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = 25 + 2*i
ndir = 25 + 2*i + 1
ext = '{:02}'.format(i + 1)
commands += _ACE_TEMPLATE_ACER.format(**locals())
tapeout[npendf] = pendf
# Indicate tapes to save for each ACER run
tapeout[nace] = fname.format(ace, temperature)
tapeout[ndir] = fname.format(xsdir, temperature)
# reconr
commands += _TEMPLATE_RECONR
nlast = npendf
# broadr
if broadr:
nbroadr = nlast + 1
commands += _TEMPLATE_BROADR
nlast = nbroadr
# heatr
if heatr:
nheatr_in = nlast
nheatr = nheatr_in + 1
commands += _TEMPLATE_HEATR
nlast = nheatr
# purr
if purr:
npurr_in = nlast
npurr = npurr_in + 1
commands += _TEMPLATE_PURR
nlast = npurr
commands = commands.format(**locals())
# acer
if acer:
nacer_in = nlast
fname = '{}_{:.1f}'
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = nacer_in + 1 + 2*i
ndir = nace + 1
ext = '{:02}'.format(i + 1)
commands += _TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = fname.format(ace, temperature)
tapeout[ndir] = fname.format(xsdir, temperature)
commands += 'stop\n'
retcode = run(commands, tapein, tapeout, **kwargs)
if retcode == 0:
if acer and retcode == 0:
with open(ace, 'w') as ace_file, open(xsdir, 'w') as xsdir_file:
for temperature in temperatures:
# Get contents of ACE file
@ -311,8 +340,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir',
return retcode
def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir='xsdir', **kwargs):
def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir='xsdir', error=0.001, **kwargs):
"""Generate thermal scattering ACE file from ENDF files
Parameters
@ -328,6 +357,8 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
Path of ACE file to write
xsdir : str, optional
Path of xsdir file to write
error : float, optional
Fractional error tolerance for NJOY processing
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -347,7 +378,6 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
data = _THERMAL_DATA[mat_thermal]
zaids = ' '.join(str(zaid) for zaid in data.zaids[:3])
energy_max = ev_thermal.info['energy_max']
# Determine name of library
library = '{}-{}.{}'.format(*ev_thermal.info['library'])
@ -359,7 +389,7 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
# Determine whether elastic is incoherent (0) or coherent (1)
file_obj = StringIO(ev_thermal.section[7, 2])
elastic_type = endf.get_head_record(file_obj)[2]
elastic_type = endf.get_head_record(file_obj)[2] - 1
else:
elastic = 0
mt_elastic = 0
@ -369,6 +399,7 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
file_obj = StringIO(ev_thermal.section[7, 4])
items = endf.get_head_record(file_obj)
items, values = endf.get_list_record(file_obj)
energy_max = values[3]
natom = int(values[5])
# Note that the 'iform' parameter is omitted in NJOY 99. We assume that the
@ -390,16 +421,41 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
num_temp = len(temperatures)
temps = ' '.join(str(i) for i in temperatures)
commands = _ACE_THERMAL_TEMPLATE.format(**locals())
tapein = {20: filename, 60: filename_thermal}
# Create njoy commands by modules
commands = ""
nendf, nthermal_endf, npendf = 20, 21, 22
tapein = {nendf: filename, nthermal_endf:filename_thermal}
tapeout = {}
# reconr
commands += _TEMPLATE_RECONR
nlast = npendf
# broadr
nbroadr = nlast + 1
commands += _TEMPLATE_BROADR
nlast = nbroadr
# thermr
nthermr1_in = nlast
nthermr1 = nthermr1_in + 1
nthermr2_in = nthermr1
nthermr2 = nthermr2_in + 1
commands += _THERMAL_TEMPLATE_THERMR
nlast = nthermr2
commands = commands.format(**locals())
# acer
nthermal_acer_in = nlast
fname = '{}_{:.1f}'
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = 28 + 2*i
ndir = 28 + 2*i + 1
nace = nthermal_acer_in + 1 + 2*i
ndir = nace + 1
ext = '{:02}'.format(i + 1)
commands += _ACE_THERMAL_TEMPLATE_ACER.format(**locals())
commands += _THERMAL_TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = fname.format(ace, temperature)

View file

@ -94,7 +94,7 @@ def _get_products(ev, mt):
# Get yield for this product
params, yield_ = get_tab1_record(file_obj)
za = params[0]
za = int(params[0])
awr = params[1]
lip = params[2]
law = params[3]
@ -106,9 +106,8 @@ def _get_products(ev, mt):
elif za == 1000:
p = Product('electron')
else:
z = za // 1000
a = za % 1000
p = Product('{}{}'.format(ATOMIC_SYMBOL[z], a))
Z, A = divmod(za, 1000)
p = Product('{}{}'.format(ATOMIC_SYMBOL[Z], A))
p.yield_ = yield_
@ -467,6 +466,93 @@ def _get_fission_products_endf(ev):
return products, derived_products
def _get_activation_products(ev, rx):
"""Generate activation products from an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
The ENDF evaluation
rx : openmc.data.Reaction
Reaction which generates activation products
Returns
-------
products : list of openmc.data.Product
Activation products
"""
file_obj = StringIO(ev.section[8, rx.mt])
# Determine total number of states and whether decay chain is given in a
# decay sublibrary
items = get_head_record(file_obj)
n_states = items[4]
decay_sublib = (items[5] == 1)
# Determine if file 9/10 are present
present = {9: False, 10: False}
for i in range(n_states):
if decay_sublib:
items = get_cont_record(file_obj)
else:
items, values = get_list_record(file_obj)
lmf = items[2]
if lmf == 9:
present[9] = True
elif lmf == 10:
present[10] = True
products = []
for mf in (9, 10):
if not present[mf]:
continue
file_obj = StringIO(ev.section[mf, rx.mt])
items = get_head_record(file_obj)
n_states = items[4]
for i in range(n_states):
# Determine what the product is
items, xs = get_tab1_record(file_obj)
Z, A = divmod(items[2], 1000)
excited_state = items[3]
# Get GND name for product
symbol = ATOMIC_SYMBOL[Z]
if excited_state > 0:
name = '{}{}_e{}'.format(symbol, A, excited_state)
else:
name = '{}{}'.format(symbol, A)
p = Product(name)
if mf == 9:
p.yield_ = xs
else:
# Re-interpolate production cross section and neutron cross
# section to union energy grid
energy = np.union1d(xs.x, rx.xs['0K'].x)
prod_xs = xs(energy)
neutron_xs = rx.xs['0K'](energy)
idx = np.where(neutron_xs > 0)
# Calculate yield as ratio
yield_ = np.zeros_like(energy)
yield_[idx] = prod_xs[idx] / neutron_xs[idx]
p.yield_ = Tabulated1D(energy, yield_)
# Check if product already exists from MF=6 and if it does, just
# overwrite the existing yield.
for product in rx.products:
if name == product.particle:
product.yield_ = p.yield_
break
else:
products.append(p)
return products
def _get_photon_products_ace(ace, rx):
"""Generate photon products from an ACE table
@ -653,7 +739,7 @@ def _get_photon_products_endf(ev, rx):
photon = Product('photon')
items, xs = get_tab1_record(file_obj)
# Re-interpolation photon production cross section and neutron cross
# Re-interpolate photon production cross section and neutron cross
# section to union energy grid
energy = np.union1d(xs.x, rx.xs['0K'].x)
photon_prod_xs = xs(energy)
@ -1086,6 +1172,9 @@ class Reaction(EqualityMixin):
else:
rx.products.append(neutron)
if (8, mt) in ev.section:
rx.products += _get_activation_products(ev, rx)
if (12, mt) in ev.section or (13, mt) in ev.section:
rx.products += _get_photon_products_endf(ev, rx)

View file

@ -8,14 +8,13 @@ from xml.etree import ElementTree as ET
from six import add_metaclass
import numpy as np
import pandas as pd
import openmc
import openmc.checkvalue as cv
from .mixin import IDManagerMixin
# "Static" variable for auto-generated Filter IDs
AUTO_FILTER_ID = 10000
_FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface',
'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal',
'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom']
@ -28,12 +27,6 @@ _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in',
11: 'z-max out', 12: 'z-max in'}
def reset_auto_filter_id():
"""Reset counter for auto-generated filter IDs."""
global AUTO_FILTER_ID
AUTO_FILTER_ID = 10000
class FilterMeta(ABCMeta):
def __new__(cls, name, bases, namespace, **kwargs):
# Check the class name.
@ -73,7 +66,7 @@ class FilterMeta(ABCMeta):
@add_metaclass(FilterMeta)
class Filter(object):
class Filter(IDManagerMixin):
"""Tally modifier that describes phase-space and other characteristics.
Parameters
@ -99,6 +92,9 @@ class Filter(object):
"""
next_id = 1
used_ids = set()
def __init__(self, bins, filter_id=None):
self.bins = bins
self.id = filter_id
@ -194,10 +190,6 @@ class Filter(object):
def bins(self):
return self._bins
@property
def id(self):
return self._id
@property
def num_bins(self):
return self._num_bins
@ -216,17 +208,6 @@ class Filter(object):
self._bins = bins
@id.setter
def id(self, filter_id):
if filter_id is None:
global AUTO_FILTER_ID
self._id = AUTO_FILTER_ID
AUTO_FILTER_ID += 1
else:
cv.check_type('filter ID', filter_id, Integral)
cv.check_greater_than('filter ID', filter_id, 0, equality=True)
self._id = filter_id
@num_bins.setter
def num_bins(self, num_bins):
cv.check_type('filter num_bins', num_bins, Integral)
@ -460,9 +441,7 @@ class Filter(object):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
filter_bins = np.repeat(self.bins, self.stride)
@ -748,9 +727,7 @@ class SurfaceFilter(Filter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
filter_bins = np.repeat(self.bins, self.stride)
@ -908,9 +885,7 @@ class MeshFilter(Filter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Initialize dictionary to build Pandas Multi-index column
@ -1156,9 +1131,7 @@ class EnergyFilter(RealFilter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Extract the lower and upper energy bounds, then repeat and tile
@ -1368,9 +1341,7 @@ class DistribcellFilter(Filter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
level_df = None
@ -1564,9 +1535,7 @@ class MuFilter(RealFilter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Extract the lower and upper energy bounds, then repeat and tile
@ -1671,9 +1640,7 @@ class PolarFilter(RealFilter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Extract the lower and upper angle bounds, then repeat and tile
@ -1778,9 +1745,7 @@ class AzimuthalFilter(RealFilter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Extract the lower and upper angle bounds, then repeat and tile
@ -2072,8 +2037,6 @@ class EnergyFunctionFilter(Filter):
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
import pandas as pd
df = pd.DataFrame()
# There is no clean way of sticking all the energy, y data into a

View file

@ -9,14 +9,6 @@ from openmc.clean_xml import sort_xml_elements, clean_xml_indentation
from openmc.checkvalue import check_type
def reset_auto_ids():
"""Reset counters for all auto-generated IDs"""
openmc.reset_auto_material_id()
openmc.reset_auto_surface_id()
openmc.reset_auto_cell_id()
openmc.reset_auto_universe_id()
class Geometry(object):
"""Geometry representing a collection of surfaces, cells, and universes.
@ -258,7 +250,7 @@ class Geometry(object):
lattices[cell.fill.id] = cell.fill
return lattices
def get_all_surfaces(self):
"""
Return all surfaces used in the geometry
@ -270,11 +262,11 @@ class Geometry(object):
"""
surfaces = OrderedDict()
for cell in self.get_all_cells().values():
surfaces = cell.region.get_surfaces(surfaces)
return surfaces
def get_materials_by_name(self, name, case_sensitive=False, matching=False):
"""Return a list of materials with matching names.
@ -482,22 +474,30 @@ class Geometry(object):
lattices.sort(key=lambda x: x.id)
return lattices
def determine_paths(self):
def determine_paths(self, instances_only=False):
"""Determine paths through CSG tree for cells and materials.
This method recursively traverses the CSG tree to determine each unique
path that reaches every cell and material. The paths are stored in the
:attr:`Cell.paths` and :attr:`Material.paths` attributes.
Parameters
----------
instances_only : bool, optional
If true, this method will only determine the number of instances of
each cell and material.
"""
# (Re-)initialize all cell instances to 0
for cell in self.get_all_cells().values():
cell._paths = []
cell._num_instances = 0
for material in self.get_all_materials().values():
material._paths = []
material._num_instances = 0
# Recursively traverse the CSG tree to count all cell instances
self.root_universe._determine_paths()
self.root_universe._determine_paths(instances_only=instances_only)
def clone(self):
"""Create a copy of this geometry with new unique IDs for all of its

View file

@ -12,10 +12,11 @@ import numpy as np
import openmc.checkvalue as cv
import openmc
from openmc.mixin import IDManagerMixin
@add_metaclass(ABCMeta)
class Lattice(object):
class Lattice(IDManagerMixin):
"""A repeating structure wherein each element is a universe.
Parameters
@ -41,6 +42,10 @@ class Lattice(object):
of the lattice
"""
next_id = 1
used_ids = openmc.Universe.used_ids
def __init__(self, lattice_id=None, name=''):
# Initialize Lattice class attributes
self.id = lattice_id
@ -68,10 +73,6 @@ class Lattice(object):
def __ne__(self, other):
return not self == other
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -88,16 +89,6 @@ class Lattice(object):
def universes(self):
return self._universes
@id.setter
def id(self, lattice_id):
if lattice_id is None:
self._id = openmc.universe.AUTO_UNIVERSE_ID
openmc.universe.AUTO_UNIVERSE_ID += 1
else:
cv.check_type('lattice ID', lattice_id, Integral)
cv.check_greater_than('lattice ID', lattice_id, 0, equality=True)
self._id = lattice_id
@name.setter
def name(self, name):
if name is not None:

View file

@ -11,16 +11,7 @@ import openmc
import openmc.data
import openmc.checkvalue as cv
from openmc.clean_xml import sort_xml_elements, clean_xml_indentation
# A static variable for auto-generated Material IDs
AUTO_MATERIAL_ID = 10000
def reset_auto_material_id():
"""Reset counter for auto-generated material IDs."""
global AUTO_MATERIAL_ID
AUTO_MATERIAL_ID = 10000
from .mixin import IDManagerMixin
# Units for density supported by OpenMC
@ -28,7 +19,7 @@ DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum',
'macro']
class Material(object):
class Material(IDManagerMixin):
"""A material composed of a collection of nuclides/elements.
To create a material, one should create an instance of this class, add
@ -89,6 +80,9 @@ class Material(object):
"""
next_id = 1
used_ids = set()
def __init__(self, material_id=None, name='', temperature=None):
# Initialize class attributes
self.id = material_id
@ -97,7 +91,8 @@ class Material(object):
self._density = None
self._density_units = ''
self._depletable = False
self._paths = []
self._paths = None
self._num_instances = None
self._volume = None
self._atoms = {}
@ -183,10 +178,6 @@ class Material(object):
return string
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -209,14 +200,18 @@ class Material(object):
@property
def paths(self):
if not self._paths:
if self._paths is None:
raise ValueError('Material instance paths have not been determined. '
'Call the Geometry.determine_paths() method.')
return self._paths
@property
def num_instances(self):
return len(self.paths)
if self._num_instances is None:
raise ValueError(
'Number of material instances have not been determined. Call '
'the Geometry.determine_paths() method.')
return self._num_instances
@property
def elements(self):
@ -259,18 +254,6 @@ class Material(object):
def volume(self):
return self._volume
@id.setter
def id(self, material_id):
if material_id is None:
global AUTO_MATERIAL_ID
self._id = AUTO_MATERIAL_ID
AUTO_MATERIAL_ID += 1
else:
cv.check_type('material ID', material_id, Integral)
cv.check_greater_than('material ID', material_id, 0, equality=True)
self._id = material_id
@name.setter
def name(self, name):
if name is not None:
@ -298,11 +281,6 @@ class Material(object):
cv.check_type('material volume', volume, Real)
self._volume = volume
@num_instances.setter
def num_instances(self, num_instances):
cv.check_type('num_instances', num_instances, Integral)
self._num_instances = num_instances
@classmethod
def from_hdf5(cls, group):
"""Create material from HDF5 group
@ -640,13 +618,18 @@ class Material(object):
if element == elm[0]:
self._elements.remove(elm)
def add_s_alpha_beta(self, name):
def add_s_alpha_beta(self, name, fraction=1.0):
r"""Add an :math:`S(\alpha,\beta)` table to the material
Parameters
----------
name : str
Name of the :math:`S(\alpha,\beta)` table
fraction : float
The fraction of relevant nuclei that are affected by the
:math:`S(\alpha,\beta)` table. For example, if the material is a
block of carbon that is 60% graphite and 40% amorphous then add a
graphite :math:`S(\alpha,\beta)` table with fraction=0.6.
"""
@ -660,13 +643,17 @@ class Material(object):
'non-string table name "{}"'.format(self._id, name)
raise ValueError(msg)
cv.check_type('S(a,b) fraction', fraction, Real)
cv.check_greater_than('S(a,b) fraction', fraction, 0.0, True)
cv.check_less_than('S(a,b) fraction', fraction, 1.0, True)
new_name = openmc.data.get_thermal_name(name)
if new_name != name:
msg = 'OpenMC S(a,b) tables follow the GND naming convention. ' \
'Table "{}" is being renamed as "{}".'.format(name, new_name)
warnings.warn(msg)
self._sab.append(new_name)
self._sab.append((new_name, fraction))
def make_isotropic_in_lab(self):
for nuclide, percent, percent_type in self._nuclides:
@ -820,8 +807,18 @@ class Material(object):
# If no nemoize'd clone exists, instantiate one
if self not in memo:
# Temporarily remove paths -- this is done so that when the clone is
# made, it doesn't create a copy of the paths (which are specific to
# an instance)
paths = self._paths
self._paths = None
clone = deepcopy(self)
clone.id = None
clone._num_instances = None
# Restore paths on original instance
self._paths = paths
# Memoize the clone
memo[self] = clone
@ -984,7 +981,9 @@ class Material(object):
if len(self._sab) > 0:
for sab in self._sab:
subelement = ET.SubElement(element, "sab")
subelement.set("name", sab)
subelement.set("name", sab[0])
if sab[1] != 1.0:
subelement.set("fraction", str(sab[1]))
return element

View file

@ -8,20 +8,10 @@ import numpy as np
import openmc.checkvalue as cv
import openmc
from openmc.mixin import EqualityMixin
from openmc.mixin import EqualityMixin, IDManagerMixin
# "Static" variable for auto-generated and Mesh IDs
AUTO_MESH_ID = 10000
def reset_auto_mesh_id():
"""Reset counter for auto-generated mesh IDs."""
global AUTO_MESH_ID
AUTO_MESH_ID = 10000
class Mesh(EqualityMixin):
class Mesh(EqualityMixin, IDManagerMixin):
"""A structured Cartesian mesh in one, two, or three dimensions
Parameters
@ -52,6 +42,9 @@ class Mesh(EqualityMixin):
"""
next_id = 1
used_ids = set()
def __init__(self, mesh_id=None, name=''):
# Initialize Mesh class attributes
self.id = mesh_id
@ -62,10 +55,6 @@ class Mesh(EqualityMixin):
self._upper_right = None
self._width = None
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -94,17 +83,6 @@ class Mesh(EqualityMixin):
def num_mesh_cells(self):
return np.prod(self._dimension)
@id.setter
def id(self, mesh_id):
if mesh_id is None:
global AUTO_MESH_ID
self._id = AUTO_MESH_ID
AUTO_MESH_ID += 1
else:
cv.check_type('mesh ID', mesh_id, Integral)
cv.check_greater_than('mesh ID', mesh_id, 0, equality=True)
self._id = mesh_id
@name.setter
def name(self, name):
if name is not None:

View file

@ -5,6 +5,9 @@ import os
from six import string_types
import numpy as np
import h5py
from scipy.interpolate import interp1d
from scipy.integrate import simps
from scipy.special import eval_legendre
import openmc
import openmc.mgxs
@ -1808,10 +1811,6 @@ class XSdata(object):
"""
from scipy.interpolate import interp1d
from scipy.integrate import simps
from scipy.special import eval_legendre
check_value('target_format', target_format, _SCATTER_TYPES)
check_type('target_order', target_order, Integral)
if target_format == 'legendre':

View file

@ -1,5 +1,10 @@
from numbers import Integral
from warnings import warn
import numpy as np
import openmc.checkvalue as cv
class EqualityMixin(object):
"""A Class which provides generic __eq__ and __ne__ functionality which
@ -18,3 +23,81 @@ class EqualityMixin(object):
def __ne__(self, other):
return not self.__eq__(other)
class IDWarning(UserWarning):
pass
class IDManagerMixin(object):
"""A Class which automatically manages unique IDs.
This mixin gives any subclass the ability to assign unique IDs through an
'id' property and keeps track of which ones have already been
assigned. Crucially, each subclass must define class variables 'next_id' and
'used_ids' as they are used in the 'id' property that is supplied here.
"""
@property
def id(self):
return self._id
@id.setter
def id(self, uid):
cls = type(self)
name = cls.__name__
if uid is None:
while cls.next_id in cls.used_ids:
cls.next_id += 1
self._id = cls.next_id
cls.used_ids.add(cls.next_id)
else:
cv.check_type('{} ID'.format(name), uid, Integral)
cv.check_greater_than('{} ID'.format(name), uid, 0, equality=True)
if uid in cls.used_ids:
msg = 'Another {} instance already exists with id={}.'.format(
name, uid)
warn(msg, IDWarning)
else:
cls.used_ids.add(uid)
self._id = uid
def reset_auto_ids():
"""Reset counters for all auto-generated IDs"""
for cls in IDManagerMixin.__subclasses__():
cls.used_ids.clear()
cls.next_id = 1
def reserve_ids(ids, cls=None):
"""Reserve a set of IDs that won't be used for auto-generated IDs.
Parameters
----------
ids : iterable of int
IDs to reserve
cls : type or None
Class for which IDs should be reserved (e.g., :class:`openmc.Cell`). If
None, all classes that have auto-generated IDs will be used.
"""
if cls is None:
for cls in IDManagerMixin.__subclasses__():
cls.used_ids |= set(ids)
else:
cls.used_ids |= set(ids)
def set_auto_id(next_id):
"""Set the next ID for auto-generated IDs.
Parameters
----------
next_id : int
The next ID to assign to objects with auto-generated IDs.
"""
for cls in IDManagerMixin.__subclasses__():
cls.next_id = next_id

View file

@ -513,7 +513,7 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
universes = np.empty(shape[::-1], dtype=openmc.Universe)
for idx, triso_list in sorted(triso_locations.items()):
if len(triso_list) > 0:
outside_trisos = openmc.Intersection(*[~t.region for t in triso_list])
outside_trisos = openmc.Intersection(~t.region for t in triso_list)
background_cell = openmc.Cell(fill=background, region=outside_trisos)
else:
background_cell = openmc.Cell(fill=background)

View file

@ -11,12 +11,12 @@ class Nuclide(object):
Parameters
----------
name : str
Name of the nuclide, e.g. U235
Name of the nuclide, e.g. 'U235'
Attributes
----------
name : str
Name of the nuclide, e.g. U235
Name of the nuclide, e.g. 'U235'
scattering : 'data' or 'iso-in-lab' or None
The type of angular scattering distribution to use

View file

@ -371,11 +371,11 @@ def get_openmoc_region(openmc_region):
openmoc.Halfspace(halfspace, get_openmoc_surface(surface))
elif isinstance(openmc_region, openmc.Intersection):
openmoc_region = openmoc.Intersection()
for openmc_node in openmc_region.nodes:
for openmc_node in openmc_region:
openmoc_region.addNode(get_openmoc_region(openmc_node))
elif isinstance(openmc_region, openmc.Union):
openmoc_region = openmoc.Union()
for openmc_node in openmc_region.nodes:
for openmc_node in openmc_region:
openmoc_region.addNode(get_openmoc_region(openmc_node))
elif isinstance(openmc_region, openmc.Complement):
openmoc_region = openmoc.Complement()
@ -411,12 +411,12 @@ def get_openmc_region(openmoc_region):
openmc_region = openmc.Intersection()
for openmoc_node in openmoc_region.getNodes():
openmc_node = get_openmc_region(openmoc_node)
openmc_region.nodes.append(openmc_node)
openmc_region.append(openmc_node)
elif openmoc_region.getRegionType() == openmoc.UNION:
openmc_region = openmc.Union()
for openmoc_node in openmoc_region.getNodes():
openmc_node = get_openmc_region(openmoc_node)
openmc_region.nodes.append(openmc_node)
openmc_region.append(openmc_node)
elif openmoc_region.getRegionType() == openmoc.COMPLEMENT:
openmoc_nodes = openmoc_region.getNodes()
openmc_node = get_openmc_region(openmoc_nodes[0])

View file

@ -10,16 +10,7 @@ import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.clean_xml import clean_xml_indentation
# A static variable for auto-generated Plot IDs
AUTO_PLOT_ID = 10000
def reset_auto_plot_id():
"""Reset counter for auto-generated plot IDs."""
global AUTO_PLOT_ID
AUTO_PLOT_ID = 10000
from openmc.mixin import IDManagerMixin
_BASES = ['xy', 'xz', 'yz']
@ -175,7 +166,7 @@ _SVG_COLORS = {
}
class Plot(object):
class Plot(IDManagerMixin):
"""Definition of a finite region of space to be plotted.
OpenMC is capable of generating two-dimensional slice plots and
@ -227,6 +218,9 @@ class Plot(object):
"""
next_id = 1
used_ids = set()
def __init__(self, plot_id=None, name=''):
# Initialize Plot class attributes
self.id = plot_id
@ -245,10 +239,6 @@ class Plot(object):
self._level = None
self._meshlines = None
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -305,17 +295,6 @@ class Plot(object):
def meshlines(self):
return self._meshlines
@id.setter
def id(self, plot_id):
if plot_id is None:
global AUTO_PLOT_ID
self._id = AUTO_PLOT_ID
AUTO_PLOT_ID += 1
else:
cv.check_type('plot ID', plot_id, Integral)
cv.check_greater_than('plot ID', plot_id, 0, equality=True)
self._id = plot_id
@name.setter
def name(self, name):
cv.check_type('plot name', name, string_types)

View file

@ -1,5 +1,5 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable, OrderedDict
from collections import Iterable, OrderedDict, MutableSequence
from copy import deepcopy
from six import add_metaclass
@ -20,21 +20,21 @@ class Region(object):
"""
def __and__(self, other):
return Intersection(self, other)
return Intersection((self, other))
def __or__(self, other):
return Union(self, other)
return Union((self, other))
def __invert__(self):
return Complement(self)
@abstractmethod
def __contains__(self, point):
return False
pass
@abstractmethod
def __str__(self):
return ''
pass
def __eq__(self, other):
if not isinstance(other, type(self)):
@ -55,7 +55,7 @@ class Region(object):
----------
surfaces: collections.OrderedDict, optional
Dictionary mapping surface IDs to :class:`openmc.Surface` instances
Returns
-------
surfaces: collections.OrderedDict
@ -150,32 +150,26 @@ class Region(object):
r2 = output.pop()
if operator == ' ':
r1 = output.pop()
if isinstance(r1, Intersection) and can_be_combined(r2):
r1.nodes.append(r2)
if isinstance(r1, Intersection):
r1 &= r2
output.append(r1)
elif isinstance(r2, Intersection) and can_be_combined(r1):
r2.nodes.insert(0, r1)
r2.insert(0, r1)
output.append(r2)
elif isinstance(r1, Intersection) and isinstance(r2, Intersection):
r1.nodes += r2.nodes
output.append(r1)
else:
output.append(Intersection(r1, r2))
output.append(r1 & r2)
elif operator == '|':
r1 = output.pop()
if isinstance(r1, Union) and can_be_combined(r2):
r1.nodes.append(r2)
if isinstance(r1, Union):
r1 |= r2
output.append(r1)
elif isinstance(r2, Union) and can_be_combined(r1):
r2.nodes.insert(0, r1)
r2.insert(0, r1)
output.append(r2)
elif isinstance(r1, Union) and isinstance(r2, Union):
r1.nodes += r2.nodes
output.append(r1)
else:
output.append(Union(r1, r2))
output.append(r1 | r2)
elif operator == '~':
output.append(Complement(r2))
output.append(~r2)
# The following is an implementation of the shunting yard algorithm to
# generate an abstract syntax tree for the region expression.
@ -248,12 +242,12 @@ class Region(object):
'the abstract region class.')
class Intersection(Region):
class Intersection(Region, MutableSequence):
r"""Intersection of two or more regions.
Instances of Intersection are generally created via the __and__ operator
applied to two instances of :class:`openmc.Region`. This is illustrated in
the following example:
Instances of Intersection are generally created via the & operator applied
to two instances of :class:`openmc.Region`. This is illustrated in the
following example:
>>> equator = openmc.ZPlane(z0=0.0)
>>> earth = openmc.Sphere(R=637.1e6)
@ -262,31 +256,51 @@ class Intersection(Region):
>>> type(northern_hemisphere)
<class 'openmc.region.Intersection'>
Instances of this class behave like a mutable sequence, e.g., they can be
indexed and have an append() method.
Parameters
----------
\*nodes
nodes : iterable of openmc.Region
Regions to take the intersection of
Attributes
----------
nodes : list of openmc.Region
Regions to take the intersection of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
def __init__(self, *nodes):
self.nodes = list(nodes)
def __init__(self, nodes):
self._nodes = list(nodes)
def __and__(self, other):
new = Intersection(*self.nodes)
new.nodes.append(other)
new = Intersection(self)
new &= other
return new
def __iter__(self):
for n in self.nodes:
yield n
def __iand__(self, other):
if isinstance(other, Intersection):
self.extend(other)
else:
self.append(other)
return self
# Implement mutable sequence protocol by delegating to list
def __getitem__(self, key):
return self._nodes[key]
def __setitem__(self, key, value):
self._nodes[key] = value
def __delitem__(self, key):
del self._nodes[key]
def __len__(self):
return len(self._nodes)
def insert(self, index, value):
self._nodes.insert(index, value)
def __contains__(self, point):
"""Check whether a point is contained in the region.
@ -302,30 +316,21 @@ class Intersection(Region):
Whether the point is in the region
"""
return all(point in n for n in self.nodes)
return all(point in n for n in self)
def __str__(self):
return '(' + ' '.join(map(str, self.nodes)) + ')'
@property
def nodes(self):
return self._nodes
return '(' + ' '.join(map(str, self)) + ')'
@property
def bounding_box(self):
lower_left = np.array([-np.inf, -np.inf, -np.inf])
upper_right = np.array([np.inf, np.inf, np.inf])
for n in self.nodes:
for n in self:
lower_left_n, upper_right_n = n.bounding_box
lower_left[:] = np.maximum(lower_left, lower_left_n)
upper_right[:] = np.minimum(upper_right, upper_right_n)
return lower_left, upper_right
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
self._nodes = nodes
def clone(self, memo=None):
"""Create a copy of this region - each of the surfaces in the
intersection's nodes will be cloned and will have new unique IDs.
@ -347,47 +352,67 @@ class Intersection(Region):
memo = {}
clone = deepcopy(self)
clone.nodes = [n.clone(memo) for n in self.nodes]
clone[:] = [n.clone(memo) for n in self]
return clone
class Union(Region):
class Union(Region, MutableSequence):
r"""Union of two or more regions.
Instances of Union are generally created via the __or__ operator applied to
two instances of :class:`openmc.Region`. This is illustrated in the
following example:
Instances of Union are generally created via the | operator applied to two
instances of :class:`openmc.Region`. This is illustrated in the following
example:
>>> s1 = openmc.ZPlane(z0=0.0)
>>> s2 = openmc.Sphere(R=637.1e6)
>>> type(-s2 | +s1)
<class 'openmc.region.Union'>
Instances of this class behave like a mutable sequence, e.g., they can be
indexed and have an append() method.
Parameters
----------
\*nodes
nodes : iterable of openmc.Region
Regions to take the union of
Attributes
----------
nodes : tuple of openmc.Region
Regions to take the union of
bounding_box : 2-tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
def __init__(self, *nodes):
self.nodes = list(nodes)
def __init__(self, nodes):
self._nodes = list(nodes)
def __or__(self, other):
new = Union(*self.nodes)
new.nodes.append(other)
new = Union(self)
new |= other
return new
def __iter__(self):
for n in self.nodes:
yield n
def __ior__(self, other):
if isinstance(other, Union):
self.extend(other)
else:
self.append(other)
return self
# Implement mutable sequence protocol by delegating to list
def __getitem__(self, key):
return self._nodes[key]
def __setitem__(self, key, value):
self._nodes[key] = value
def __delitem__(self, key):
del self._nodes[key]
def __len__(self):
return len(self._nodes)
def insert(self, index, value):
self._nodes.insert(index, value)
def __contains__(self, point):
"""Check whether a point is contained in the region.
@ -403,30 +428,21 @@ class Union(Region):
Whether the point is in the region
"""
return any(point in n for n in self.nodes)
return any(point in n for n in self)
def __str__(self):
return '(' + ' | '.join(map(str, self.nodes)) + ')'
@property
def nodes(self):
return self._nodes
return '(' + ' | '.join(map(str, self)) + ')'
@property
def bounding_box(self):
lower_left = np.array([np.inf, np.inf, np.inf])
upper_right = np.array([-np.inf, -np.inf, -np.inf])
for n in self.nodes:
for n in self:
lower_left_n, upper_right_n = n.bounding_box
lower_left[:] = np.minimum(lower_left, lower_left_n)
upper_right[:] = np.maximum(upper_right, upper_right_n)
return lower_left, upper_right
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
self._nodes = nodes
def clone(self, memo=None):
"""Create a copy of this region - each of the surfaces in the
union's nodes will be cloned and will have new unique IDs.
@ -448,7 +464,7 @@ class Union(Region):
memo = {}
clone = copy.deepcopy(self)
clone.nodes = [n.clone(memo) for n in self.nodes]
clone[:] = [n.clone(memo) for n in self]
return clone
@ -456,7 +472,7 @@ class Complement(Region):
"""Complement of a region.
The Complement of an existing :class:`openmc.Region` can be created by using
the __invert__ operator as the following example demonstrates:
the ~ operator as the following example demonstrates:
>>> xl = openmc.XPlane(x0=-10.0)
>>> xr = openmc.XPlane(x0=10.0)
@ -518,9 +534,9 @@ class Complement(Region):
# only applies to surface half-spaces, thus allowing us to calculate the
# bounding box in the usual recursive manner.
if isinstance(self.node, Union):
temp_region = Intersection(*[~n for n in self.node.nodes])
temp_region = Intersection(~n for n in self.node)
elif isinstance(self.node, Intersection):
temp_region = Union(*[~n for n in self.node.nodes])
temp_region = Union(~n for n in self.node)
elif isinstance(self.node, Complement):
temp_region = self.node.node
else:

View file

@ -1,6 +1,8 @@
from collections import Callable
from numbers import Real
import scipy.optimize as sopt
import openmc
import openmc.model
import openmc.checkvalue as cv
@ -144,8 +146,6 @@ def search_for_keff(model_builder, initial_guess=None, target=1.0,
model = model_builder(initial_guess, **model_args)
cv.check_type('model_builder return', model, openmc.model.Model)
import scipy.optimize as sopt
# Set the iteration data storage variables
guesses = []
results = []

View file

@ -125,13 +125,16 @@ class Settings(object):
temperature : dict
Defines a default temperature and method for treating intermediate
temperatures at which nuclear data doesn't exist. Accepted keys are
'default', 'method', 'tolerance', and 'multipole'. The value for
'default' should be a float representing the default temperature in
'default', 'method', 'range', 'tolerance', and 'multipole'. The value
for 'default' should be a float representing the default temperature in
Kelvin. The value for 'method' should be 'nearest' or 'interpolation'.
If the method is 'nearest', 'tolerance' indicates a range of temperature
within which cross sections may be used. 'multipole' is a boolean
indicating whether or not the windowed multipole method should be used
to evaluate resolved resonance cross sections.
within which cross sections may be used. The value for 'range' should be
a pair a minimum and maximum temperatures which are used to indicate
that cross sections be loaded at all temperatures within the
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
@ -639,7 +642,8 @@ class Settings(object):
cv.check_type('temperature settings', temperature, Mapping)
for key, value in temperature.items():
cv.check_value('temperature key', key,
['default', 'method', 'tolerance', 'multipole'])
['default', 'method', 'tolerance', 'multipole',
'range'])
if key == 'default':
cv.check_type('default temperature', value, Real)
elif key == 'method':
@ -649,6 +653,11 @@ class Settings(object):
cv.check_type('temperature tolerance', value, Real)
elif key == 'multipole':
cv.check_type('temperature multipole', value, bool)
elif key == 'range':
cv.check_length('temperature range', value, 2)
for T in value:
cv.check_type('temperature', T, Real)
self._temperature = temperature
@threads.setter
@ -999,6 +1008,8 @@ class Settings(object):
"temperature_{}".format(key))
if isinstance(value, bool):
element.text = str(value).lower()
elif key == 'range':
element.text = ' '.join(str(T) for T in value)
else:
element.text = str(value)

View file

@ -356,71 +356,76 @@ class StatePoint(object):
else:
tally_ids = []
# Iterate over all tallies
for tally_id in tally_ids:
group = tallies_group['tally {}'.format(tally_id)]
# Ignore warnings about duplicate IDs
with warnings.catch_warnings():
warnings.simplefilter('ignore', openmc.IDWarning)
# Read the number of realizations
n_realizations = group['n_realizations'].value
# Iterate over all tallies
for tally_id in tally_ids:
group = tallies_group['tally {}'.format(tally_id)]
# 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.num_realizations = n_realizations
# Read the number of realizations
n_realizations = group['n_realizations'].value
# Read derivative information.
if 'derivative' in group:
deriv_id = group['derivative'].value
tally.derivative = self.tally_derivatives[deriv_id]
# 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.num_realizations = n_realizations
# Read all filters
n_filters = group['n_filters'].value
if n_filters > 0:
filter_ids = group['filters'].value
filters_group = self._f['tallies/filters']
for filter_id in filter_ids:
filter_group = filters_group['filter {}'.format(filter_id)]
new_filter = openmc.Filter.from_hdf5(filter_group,
meshes=self.meshes)
tally.filters.append(new_filter)
# Read derivative information.
if 'derivative' in group:
deriv_id = group['derivative'].value
tally.derivative = self.tally_derivatives[deriv_id]
# Read nuclide bins
nuclide_names = group['nuclides'].value
# Read all filters
n_filters = group['n_filters'].value
if n_filters > 0:
filter_ids = group['filters'].value
filters_group = self._f['tallies/filters']
for filter_id in filter_ids:
filter_group = filters_group['filter {}'.format(
filter_id)]
new_filter = openmc.Filter.from_hdf5(
filter_group, meshes=self.meshes)
tally.filters.append(new_filter)
# Add all nuclides to the Tally
for name in nuclide_names:
nuclide = openmc.Nuclide(name.decode().strip())
tally.nuclides.append(nuclide)
# Read nuclide bins
nuclide_names = group['nuclides'].value
scores = group['score_bins'].value
n_score_bins = group['n_score_bins'].value
# Add all nuclides to the Tally
for name in nuclide_names:
nuclide = openmc.Nuclide(name.decode().strip())
tally.nuclides.append(nuclide)
# Compute and set the filter strides
for i in range(n_filters):
tally_filter = tally.filters[i]
tally_filter.stride = n_score_bins * len(nuclide_names)
scores = group['score_bins'].value
n_score_bins = group['n_score_bins'].value
for j in range(i+1, n_filters):
tally_filter.stride *= tally.filters[j].num_bins
# Compute and set the filter strides
for i in range(n_filters):
tally_filter = tally.filters[i]
tally_filter.stride = n_score_bins * len(nuclide_names)
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = group['moment_orders'].value
for j in range(i+1, n_filters):
tally_filter.stride *= tally.filters[j].num_bins
# Add the scores to the Tally
for j, score in enumerate(scores):
score = score.decode()
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = group['moment_orders'].value
# If this is a moment, use generic moment order
pattern = r'-n$|-pn$|-yn$'
score = re.sub(pattern, '-' + moments[j].decode(), score)
# Add the scores to the Tally
for j, score in enumerate(scores):
score = score.decode()
tally.scores.append(score)
# If this is a moment, use generic moment order
pattern = r'-n$|-pn$|-yn$'
score = re.sub(pattern, '-' + moments[j].decode(), score)
# Add Tally to the global dictionary of all Tallies
tally.sparse = self.sparse
self._tallies[tally_id] = tally
tally.scores.append(score)
# Add Tally to the global dictionary of all Tallies
tally.sparse = self.sparse
self._tallies[tally_id] = tally
self._tallies_read = True

View file

@ -1,5 +1,6 @@
from collections import Iterable
import re
import warnings
import numpy as np
import h5py
@ -31,8 +32,6 @@ class Summary(object):
"""
def __init__(self, filename):
openmc.reset_auto_ids()
if not filename.endswith(('.h5', '.hdf5')):
msg = 'Unable to open "{0}" which is not an HDF5 summary file'
raise ValueError(msg)
@ -52,7 +51,9 @@ class Summary(object):
self._nuclides = {}
self._read_nuclides()
self._read_geometry()
with warnings.catch_warnings():
warnings.simplefilter("ignore", openmc.IDWarning)
self._read_geometry()
@property
def date_and_time(self):
@ -69,7 +70,7 @@ class Summary(object):
@property
def nuclides(self):
return self._nuclides
@property
def version(self):
return tuple(self._f.attrs['openmc_version'])

View file

@ -12,6 +12,7 @@ import numpy as np
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.region import Region, Intersection, Union
from openmc.mixin import IDManagerMixin
# A static variable for auto-generated Surface IDs
@ -20,13 +21,7 @@ AUTO_SURFACE_ID = 10000
_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic']
def reset_auto_surface_id():
"""Reset counters for all auto-generated surface IDs"""
global AUTO_SURFACE_ID
AUTO_SURFACE_ID = 10000
class Surface(object):
class Surface(IDManagerMixin):
"""An implicit surface with an associated boundary condition.
An implicit surface is defined as the set of zeros of a function of the
@ -64,6 +59,9 @@ class Surface(object):
"""
next_id = 1
used_ids = set()
def __init__(self, surface_id=None, boundary_type='transmission', name=''):
self.id = surface_id
self.name = name
@ -105,10 +103,6 @@ class Surface(object):
return string
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -125,17 +119,6 @@ class Surface(object):
def coefficients(self):
return self._coefficients
@id.setter
def id(self, surface_id):
if surface_id is None:
global AUTO_SURFACE_ID
self._id = AUTO_SURFACE_ID
AUTO_SURFACE_ID += 1
else:
check_type('surface ID', surface_id, Integral)
check_greater_than('surface ID', surface_id, 0, equality=True)
self._id = surface_id
@name.setter
def name(self, name):
if name is not None:
@ -1860,15 +1843,15 @@ class Halfspace(Region):
def __and__(self, other):
if isinstance(other, Intersection):
return Intersection(self, *other.nodes)
return Intersection([self] + other[:])
else:
return Intersection(self, other)
return Intersection((self, other))
def __or__(self, other):
if isinstance(other, Union):
return Union(self, *other.nodes)
return Union([self] + other[:])
else:
return Union(self, other)
return Union((self, other))
def __invert__(self):
return -self.surface if self.side == '+' else +self.surface

View file

@ -11,16 +11,16 @@ from xml.etree import ElementTree as ET
from six import string_types
import numpy as np
import pandas as pd
import scipy.sparse as sps
import h5py
import openmc
import openmc.checkvalue as cv
from openmc.clean_xml import clean_xml_indentation
from .mixin import IDManagerMixin
# "Static" variable for auto-generated Tally IDs
AUTO_TALLY_ID = 10000
# The tally arithmetic product types. The tensor product performs the full
# cross product of the data in two tallies with respect to a specified axis
# (filters, nuclides, or scores). The entrywise product performs the arithmetic
@ -39,13 +39,7 @@ _FILTER_CLASSES = (openmc.Filter, openmc.CrossFilter, openmc.AggregateFilter)
ESTIMATOR_TYPES = ['tracklength', 'collision', 'analog']
def reset_auto_tally_id():
"""Reset counter for auto-generated tally IDs."""
global AUTO_TALLY_ID
AUTO_TALLY_ID = 10000
class Tally(object):
class Tally(IDManagerMixin):
"""A tally defined by a set of scores that are accumulated for a list of
nuclides given a set of filters.
@ -107,6 +101,9 @@ class Tally(object):
"""
next_id = 1
used_ids = set()
def __init__(self, tally_id=None, name=''):
# Initialize Tally class attributes
self.id = tally_id
@ -202,10 +199,6 @@ class Tally(object):
return string
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -295,8 +288,6 @@ class Tally(object):
# Convert NumPy arrays to SciPy sparse LIL matrices
if self.sparse:
import scipy.sparse as sps
self._sum = \
sps.lil_matrix(self._sum.flatten(), self._sum.shape)
self._sum_sq = \
@ -337,8 +328,6 @@ class Tally(object):
# Convert NumPy array to SciPy sparse LIL matrix
if self.sparse:
import scipy.sparse as sps
self._mean = \
sps.lil_matrix(self._mean.flatten(), self._mean.shape)
@ -361,8 +350,6 @@ class Tally(object):
# Convert NumPy array to SciPy sparse LIL matrix
if self.sparse:
import scipy.sparse as sps
self._std_dev = \
sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
@ -420,17 +407,6 @@ class Tally(object):
DeprecationWarning)
self.triggers.append(trigger)
@id.setter
def id(self, tally_id):
if tally_id is None:
global AUTO_TALLY_ID
self._id = AUTO_TALLY_ID
AUTO_TALLY_ID += 1
else:
cv.check_type('tally ID', tally_id, Integral)
cv.check_greater_than('tally ID', tally_id, 0, equality=True)
self._id = tally_id
@name.setter
def name(self, name):
if name is not None:
@ -608,8 +584,6 @@ class Tally(object):
# Convert NumPy arrays to SciPy sparse LIL matrices
if sparse and not self.sparse:
import scipy.sparse as sps
if self._sum is not None:
self._sum = \
sps.lil_matrix(self._sum.flatten(), self._sum.shape)
@ -1610,7 +1584,6 @@ class Tally(object):
raise KeyError(msg)
# Initialize a pandas dataframe for the tally data
import pandas as pd
df = pd.DataFrame()
# Find the total length of the tally data array

View file

@ -7,18 +7,10 @@ from xml.etree import ElementTree as ET
from six import string_types
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.mixin import EqualityMixin, IDManagerMixin
# "Static" variable for auto-generated TallyDerivative IDs
AUTO_TALLY_DERIV_ID = 10000
def reset_auto_tally_deriv_id():
global AUTO_TALLY_ID
AUTO_TALLY_DERIV_ID = 10000
class TallyDerivative(EqualityMixin):
class TallyDerivative(EqualityMixin, IDManagerMixin):
"""A material perturbation derivative to apply to a tally.
Parameters
@ -48,6 +40,9 @@ class TallyDerivative(EqualityMixin):
"""
next_id = 1
used_ids = set()
def __init__(self, derivative_id=None, variable=None, material=None,
nuclide=None):
# Initialize Tally class attributes
@ -74,10 +69,6 @@ class TallyDerivative(EqualityMixin):
return string
@property
def id(self):
return self._id
@property
def variable(self):
return self._variable
@ -90,18 +81,6 @@ class TallyDerivative(EqualityMixin):
def nuclide(self):
return self._nuclide
@id.setter
def id(self, deriv_id):
if deriv_id is None:
global AUTO_TALLY_DERIV_ID
self._id = AUTO_TALLY_DERIV_ID
AUTO_TALLY_DERIV_ID += 1
else:
cv.check_type('tally derivative ID', deriv_id, Integral)
cv.check_greater_than('tally derivative ID', deriv_id, 0,
equality=True)
self._id = deriv_id
@variable.setter
def variable(self, var):
if var is not None:

View file

@ -11,19 +11,10 @@ import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.plots import _SVG_COLORS
from openmc.mixin import IDManagerMixin
# A static variable for auto-generated Lattice (Universe) IDs
AUTO_UNIVERSE_ID = 10000
def reset_auto_universe_id():
"""Reset counter for auto-generated universe IDs."""
global AUTO_UNIVERSE_ID
AUTO_UNIVERSE_ID = 10000
class Universe(object):
class Universe(IDManagerMixin):
"""A collection of cells that can be repeated.
Parameters
@ -55,6 +46,9 @@ class Universe(object):
"""
next_id = 1
used_ids = set()
def __init__(self, universe_id=None, name='', cells=None):
# Initialize Cell class attributes
self.id = universe_id
@ -95,10 +89,6 @@ class Universe(object):
list(self._cells.keys()))
return string
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@ -116,22 +106,11 @@ class Universe(object):
regions = [c.region for c in self.cells.values()
if c.region is not None]
if regions:
return openmc.Union(*regions).bounding_box
return openmc.Union(regions).bounding_box
else:
# Infinite bounding box
return openmc.Intersection().bounding_box
@id.setter
def id(self, universe_id):
if universe_id is None:
global AUTO_UNIVERSE_ID
self._id = AUTO_UNIVERSE_ID
AUTO_UNIVERSE_ID += 1
else:
cv.check_type('universe ID', universe_id, Integral)
cv.check_greater_than('universe ID', universe_id, 0, equality=True)
self._id = universe_id
@name.setter
def name(self, name):
if name is not None:
@ -566,7 +545,7 @@ class Universe(object):
cell_element.set("universe", str(self._id))
xml_element.append(cell_element)
def _determine_paths(self, path=''):
def _determine_paths(self, path='', instances_only=False):
"""Count the number of instances for each cell in the universe, and
record the count in the :attr:`Cell.num_instances` properties."""
@ -577,7 +556,7 @@ class Universe(object):
# If universe-filled, recursively count cells in filling universe
if cell.fill_type == 'universe':
cell.fill._determine_paths(cell_path + '->')
cell.fill._determine_paths(cell_path + '->', instances_only)
# If lattice-filled, recursively call for all universes in lattice
elif cell.fill_type == 'lattice':
@ -588,18 +567,22 @@ class Universe(object):
latt_path = '{}->l{}({})->'.format(
cell_path, latt.id, ",".join(str(x) for x in index))
univ = latt.get_universe(index)
univ._determine_paths(latt_path)
univ._determine_paths(latt_path, instances_only)
else:
if cell.fill_type == 'material':
mat = cell.fill
elif cell.fill_type == 'distribmat':
mat = cell.fill[len(cell._paths)]
mat = cell.fill[cell._num_instances]
else:
mat = None
if mat is not None:
mat._paths.append('{}->m{}'.format(cell_path, mat.id))
mat._num_instances += 1
if not instances_only:
mat._paths.append('{}->m{}'.format(cell_path, mat.id))
# Append current path
cell._paths.append(cell_path)
cell._num_instances += 1
if not instances_only:
cell._paths.append(cell_path)

View file

@ -1,7 +1,7 @@
from collections import Iterable, Mapping, OrderedDict
from numbers import Real, Integral
from xml.etree import ElementTree as ET
from warnings import warn
import warnings
import numpy as np
import pandas as pd
@ -85,15 +85,16 @@ class VolumeCalculation(object):
continue
if (np.any(np.asarray(lower_left) > ll) or
np.any(np.asarray(upper_right) < ur)):
warn("Specified bounding box is smaller than computed "
"bounding box for cell {}. Volume calculation may "
"be incorrect!".format(c.id))
warnings.warn(
"Specified bounding box is smaller than computed "
"bounding box for cell {}. Volume calculation may "
"be incorrect!".format(c.id))
self.lower_left = lower_left
self.upper_right = upper_right
else:
if self.domain_type == 'cell':
ll, ur = openmc.Union(*[c.region for c in domains]).bounding_box
ll, ur = openmc.Union(c.region for c in domains).bounding_box
if np.any(np.isinf(ll)) or np.any(np.isinf(ur)):
raise ValueError('Could not automatically determine bounding '
'box for stochastic volume calculation.')
@ -221,12 +222,14 @@ class VolumeCalculation(object):
# Instantiate some throw-away domains that are used by the constructor
# to assign IDs
if domain_type == 'cell':
domains = [openmc.Cell(uid) for uid in ids]
elif domain_type == 'material':
domains = [openmc.Material(uid) for uid in ids]
elif domain_type == 'universe':
domains = [openmc.Universe(uid) for uid in ids]
with warnings.catch_warnings():
warnings.simplefilter('ignore', openmc.IDWarning)
if domain_type == 'cell':
domains = [openmc.Cell(uid) for uid in ids]
elif domain_type == 'material':
domains = [openmc.Material(uid) for uid in ids]
elif domain_type == 'universe':
domains = [openmc.Universe(uid) for uid in ids]
# Instantiate the class and assign results
vol = cls(domains, samples, lower_left, upper_right)

View file

@ -247,6 +247,11 @@ def update_materials(root):
if 'name' in nuclide.attrib:
nucname = nuclide.attrib['name']
nucname = nucname.replace('-', '')
# If a nuclide name is in the ZAID notation (e.g., a number),
# convert it to the proper nuclide name.
if nucname.strip().isnumeric():
nucname = \
openmc.data.neutron._get_metadata(int(nucname))[0]
nucname = nucname.replace('Nat', '0')
if nucname.endswith('m'):
nucname = nucname[:-1] + '_m1'

View file

@ -1,6 +1,7 @@
#!/usr/bin/env python
import glob
import sys
import numpy as np
try:
from setuptools import setup
@ -16,6 +17,12 @@ except ImportError:
have_cython = False
# Determine shared library suffix
if sys.platform == 'darwin':
suffix = 'dylib'
else:
suffix = 'so'
# Get version information from __init__.py. This is ugly, but more reliable than
# using an import.
with open('openmc/__init__.py', 'r') as f:
@ -27,6 +34,12 @@ kwargs = {'name': 'openmc',
'openmc.stats'],
'scripts': glob.glob('scripts/openmc-*'),
# Data files and librarries
'package_data': {
'openmc': ['_libopenmc.{}'.format(suffix)],
'openmc.data': ['mass.mas12', 'fission_Q_data_endfb71.h5']
},
# Metadata
'author': 'Will Boyd',
'author_email': 'wbinventor@gmail.com',
@ -45,22 +58,16 @@ kwargs = {'name': 'openmc',
if have_setuptools:
kwargs.update({
# Required dependencies
'install_requires': ['six', 'numpy>=1.9', 'h5py'],
'install_requires': ['six', 'numpy>=1.9', 'h5py', 'scipy', 'pandas>=0.17.0'],
# Optional dependencies
'extras_require': {
'decay': ['uncertainties'],
'pandas': ['pandas>=0.17.0'],
'plot': ['matplotlib', 'ipython'],
'sparse' : ['scipy'],
'vtk': ['vtk', 'silomesh'],
'validate': ['lxml']
},
# Data files
'package_data': {
'openmc.data': ['mass.mas12', 'fission_Q_data_endfb71.h5']
},
})
# If Cython is present, add resonance reconstruction capability

895
src/api.F90 Normal file
View file

@ -0,0 +1,895 @@
module openmc_api
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, h5tclose_f, h5close_f
use constants, only: K_BOLTZMANN
use eigenvalue, only: k_sum, openmc_get_keff
use geometry, only: find_cell
use geometry_header, only: root_universe
use global
use hdf5_interface
use message_passing
use initialize, only: openmc_init
use input_xml, only: assign_0K_elastic_scattering, check_data_version
use particle_header, only: Particle
use plot, only: openmc_plot_geometry
use random_lcg, only: seed, initialize_prng
use simulation, only: openmc_run
use volume_calc, only: openmc_calculate_volumes
implicit none
private
public :: openmc_calculate_volumes
public :: openmc_cell_get_id
public :: openmc_cell_set_temperature
public :: openmc_finalize
public :: openmc_find
public :: openmc_get_cell
public :: openmc_get_keff
public :: openmc_get_material
public :: openmc_get_nuclide
public :: openmc_get_tally
public :: openmc_hard_reset
public :: openmc_init
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
public :: openmc_material_get_densities
public :: openmc_material_set_density
public :: openmc_material_set_densities
public :: openmc_nuclide_name
public :: openmc_plot_geometry
public :: openmc_reset
public :: openmc_run
public :: openmc_tally_get_id
public :: openmc_tally_get_nuclides
public :: openmc_tally_results
public :: openmc_tally_set_nuclides
! Error codes
integer(C_INT), public, bind(C) :: E_UNASSIGNED = -1
integer(C_INT), public, bind(C) :: E_OUT_OF_BOUNDS = -2
integer(C_INT), public, bind(C) :: E_CELL_NOT_ALLOCATED = -3
integer(C_INT), public, bind(C) :: E_CELL_INVALID_ID = -4
integer(C_INT), public, bind(C) :: E_CELL_NOT_FOUND = -5
integer(C_INT), public, bind(C) :: E_NUCLIDE_NOT_ALLOCATED = -6
integer(C_INT), public, bind(C) :: E_NUCLIDE_NOT_LOADED = -7
integer(C_INT), public, bind(C) :: E_NUCLIDE_NOT_IN_LIBRARY = -8
integer(C_INT), public, bind(C) :: E_MATERIAL_NOT_ALLOCATED = -9
integer(C_INT), public, bind(C) :: E_MATERIAL_INVALID_ID = -10
integer(C_INT), public, bind(C) :: E_TALLY_NOT_ALLOCATED = -11
integer(C_INT), public, bind(C) :: E_TALLY_INVALID_ID = -12
integer(C_INT), public, bind(C) :: E_INVALID_SIZE = -13
integer(C_INT), public, bind(C) :: E_CELL_NO_MATERIAL = -14
! Warning codes
integer(C_INT), public, bind(C) :: W_BELOW_MIN_BOUND = 1
integer(C_INT), public, bind(C) :: W_ABOVE_MAX_BOUND = 2
contains
!===============================================================================
! OPENMC_CELL_GET_ID returns the ID of a cell
!===============================================================================
function openmc_cell_get_id(index, id) result(err) bind(C)
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(cells)) then
id = cells(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_cell_get_id
!===============================================================================
! OPENMC_CELL_SET_TEMPERATURE sets the temperature of a cell
!===============================================================================
function openmc_cell_set_temperature(index, T, instance) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index ! cell index in cells
real(C_DOUBLE), value, intent(in) :: T ! temperature
integer(C_INT32_T), optional, intent(in) :: instance ! cell instance
integer(C_INT) :: err ! error code
integer :: j ! looping variable
integer :: n ! number of cell instances
integer :: material_ID ! material associated with cell
integer :: material_index ! material index in materials array
integer :: num_nuclides ! num nuclides in material
integer :: nuclide_index ! index of nuclide in nuclides array
real(8) :: min_temp ! min common-denominator avail temp
real(8) :: max_temp ! max common-denominator avail temp
real(8) :: temp ! actual temp we'll assign
logical :: outside_low ! lower than available data
logical :: outside_high ! higher than available data
outside_low = .false.
outside_high = .false.
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(cells)) then
! error if the cell is filled with another universe
if (cells(index) % fill /= NONE) then
err = E_CELL_NO_MATERIAL
else
! find which material is associated with this cell
if (present(instance)) then
material_ID = cells(index) % material(instance)
else
material_ID = cells(index) % material(1)
end if
! index of that material into the materials array
material_index = material_dict % get_key(material_ID)
! number of nuclides associated with this material
num_nuclides = size(materials(material_index) % nuclide)
min_temp = ZERO
max_temp = INFINITY
do j = 1, num_nuclides
nuclide_index = materials(material_index) % nuclide(j)
min_temp = max(min_temp, minval(nuclides(nuclide_index) % kTs))
max_temp = min(max_temp, maxval(nuclides(nuclide_index) % kTs))
end do
! adjust the temperature to be within bounds if necessary
if (K_BOLTZMANN * T < min_temp) then
outside_low = .true.
temp = min_temp / K_BOLTZMANN
else if (K_BOLTZMANN * T > max_temp) then
outside_high = .true.
temp = max_temp / K_BOLTZMANN
else
temp = T
end if
associate (c => cells(index))
if (allocated(c % sqrtkT)) then
n = size(c % sqrtkT)
if (present(instance) .and. n > 1) then
if (instance >= 0 .and. instance < n) then
c % sqrtkT(instance + 1) = sqrt(K_BOLTZMANN * temp)
err = 0
end if
else
c % sqrtkT(:) = sqrt(K_BOLTZMANN * temp)
err = 0
end if
end if
end associate
! Assign error codes for outside of temperature bounds provided the
! temperature was changed correctly. This needs to be done after
! changing the temperature based on the logical structure above.
if (err == 0) then
if (outside_low) err = W_BELOW_MIN_BOUND
if (outside_high) err = W_ABOVE_MAX_BOUND
end if
end if
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_cell_set_temperature
!===============================================================================
! OPENMC_FINALIZE frees up memory by deallocating arrays and resetting global
! variables
!===============================================================================
subroutine openmc_finalize() bind(C)
integer :: err
! Clear results
call openmc_reset()
! Reset global variables
assume_separate = .false.
check_overlaps = .false.
confidence_intervals = .false.
create_fission_neutrons = .true.
energy_cutoff = ZERO
energy_max_neutron = INFINITY
energy_min_neutron = ZERO
entropy_on = .false.
gen_per_batch = 1
i_user_tallies = -1
i_cmfd_tallies = -1
keff = ONE
legendre_to_tabular = .true.
legendre_to_tabular_points = 33
n_batch_interval = 1
n_filters = 0
n_meshes = 0
n_particles = 0
n_source_points = 0
n_state_points = 0
n_tallies = 0
n_user_filters = 0
n_user_meshes = 0
n_user_tallies = 0
output_summary = .true.
output_tallies = .true.
particle_restart_run = .false.
pred_batches = .false.
reduce_tallies = .true.
res_scat_on = .false.
res_scat_method = RES_SCAT_ARES
res_scat_energy_min = 0.01_8
res_scat_energy_max = 1000.0_8
restart_run = .false.
root_universe = -1
run_CE = .true.
run_mode = NONE
satisfy_triggers = .false.
seed = 1_8
source_latest = .false.
source_separate = .false.
source_write = .true.
survival_biasing = .false.
temperature_default = 293.6_8
temperature_method = TEMPERATURE_NEAREST
temperature_multipole = .false.
temperature_range = [ZERO, ZERO]
temperature_tolerance = 10.0_8
total_gen = 0
trigger_on = .false.
ufs = .false.
urr_ptables_on = .true.
verbosity = 7
weight_cutoff = 0.25_8
weight_survive = ONE
write_all_tracks = .false.
write_initial_source = .false.
! Deallocate arrays
call free_memory()
! Release compound datatypes
call h5tclose_f(hdf5_bank_t, err)
! Close FORTRAN interface.
call h5close_f(err)
#ifdef MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, err)
#endif
end subroutine openmc_finalize
!===============================================================================
! OPENMC_FIND determines the ID or a cell or material at a given point in space
!===============================================================================
function openmc_find(xyz, rtype, id, instance) result(err) bind(C)
real(C_DOUBLE), intent(in) :: xyz(3) ! Cartesian point
integer(C_INT), intent(in), value :: rtype ! 1 for cell, 2 for material
integer(C_INT32_T), intent(out) :: id
integer(C_INT32_T), intent(out) :: instance
integer(C_INT) :: err
logical :: found
type(Particle) :: p
call p % initialize()
p % coord(1) % xyz(:) = xyz
p % coord(1) % uvw(:) = [ZERO, ZERO, ONE]
call find_cell(p, found)
id = -1
instance = -1
err = E_UNASSIGNED
if (found) then
if (rtype == 1) then
id = cells(p % coord(p % n_coord) % cell) % id
elseif (rtype == 2) then
if (p % material == MATERIAL_VOID) then
id = 0
else
id = materials(p % material) % id
end if
end if
instance = p % cell_instance - 1
err = 0
else
err = E_CELL_NOT_FOUND
end if
end function openmc_find
!===============================================================================
! OPENMC_GET_CELL returns the index in the cells array of a cell with a given ID
!===============================================================================
function openmc_get_cell(id, index) result(err) bind(C)
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(cells)) then
if (cell_dict % has_key(id)) then
index = cell_dict % get_key(id)
err = 0
else
err = E_CELL_INVALID_ID
end if
else
err = E_CELL_NOT_ALLOCATED
end if
end function openmc_get_cell
!===============================================================================
! OPENMC_GET_MATERIAL returns the index in the materials array of a material
! with a given ID
!===============================================================================
function openmc_get_material(id, index) result(err) bind(C)
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(materials)) then
if (material_dict % has_key(id)) then
index = material_dict % get_key(id)
err = 0
else
err = E_MATERIAL_INVALID_ID
end if
else
err = E_MATERIAL_NOT_ALLOCATED
end if
end function openmc_get_material
!===============================================================================
! OPENMC_GET_NUCLIDE returns the index in the nuclides array of a nuclide
! with a given name
!===============================================================================
function openmc_get_nuclide(name, index) result(err) bind(C)
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT), intent(out) :: index
integer(C_INT) :: err
character(:), allocatable :: name_
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
if (allocated(nuclides)) then
if (nuclide_dict % has_key(to_lower(name_))) then
index = nuclide_dict % get_key(to_lower(name_))
err = 0
else
err = E_NUCLIDE_NOT_LOADED
end if
else
err = E_NUCLIDE_NOT_ALLOCATED
end if
end function openmc_get_nuclide
!===============================================================================
! OPENMC_GET_TALLY returns the index in the tallies array of a tally
! with a given ID
!===============================================================================
function openmc_get_tally(id, index) result(err) bind(C)
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(tallies)) then
if (tally_dict % has_key(id)) then
index = tally_dict % get_key(id)
err = 0
else
err = E_TALLY_INVALID_ID
end if
else
err = E_TALLY_NOT_ALLOCATED
end if
end function openmc_get_tally
!===============================================================================
! OPENMC_HARD_RESET reset tallies and timers as well as the pseudorandom
! generator state
!===============================================================================
subroutine openmc_hard_reset() bind(C)
! Reset all tallies and timers
call openmc_reset()
! Reset total generations and keff guess
keff = ONE
total_gen = 0
! Reset the random number generator state
seed = 1_8
call initialize_prng()
end subroutine openmc_hard_reset
!===============================================================================
! OPENMC_LOAD_NUCLIDE loads a nuclide from the cross section library
!===============================================================================
function openmc_load_nuclide(name) result(err) bind(C)
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT) :: err
integer :: i_library
integer :: n
integer(HID_T) :: file_id
integer(HID_T) :: group_id
character(:), allocatable :: name_
real(8) :: minmax(2) = [ZERO, INFINITY]
type(VectorReal) :: temperature
type(Nuclide), allocatable :: new_nuclides(:)
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
err = 0
if (.not. nuclide_dict % has_key(to_lower(name_))) then
if (library_dict % has_key(to_lower(name_))) then
! allocate extra space in nuclides array
n = n_nuclides_total
allocate(new_nuclides(n + 1))
new_nuclides(1:n) = nuclides(:)
call move_alloc(FROM=new_nuclides, TO=nuclides)
n = n + 1
i_library = library_dict % get_key(to_lower(name_))
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
call check_data_version(file_id)
! Read nuclide data from HDF5
group_id = open_group(file_id, name_)
call nuclides(n) % from_hdf5(group_id, temperature, &
temperature_method, temperature_tolerance, minmax, &
master)
call close_group(group_id)
call file_close(file_id)
! Add entry to nuclide dictionary
call nuclide_dict % add_key(to_lower(name_), n)
n_nuclides_total = n
! Assign resonant scattering data
if (res_scat_on) call assign_0K_elastic_scattering(nuclides(n))
! Initialize nuclide grid
call nuclides(n) % init_grid(energy_min_neutron, &
energy_max_neutron, n_log_bins)
else
err = E_NUCLIDE_NOT_IN_LIBRARY
end if
end if
end function openmc_load_nuclide
!===============================================================================
! OPENMC_MATERIAL_ADD_NUCLIDE
!===============================================================================
function openmc_material_add_nuclide(index, name, density) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
character(kind=C_CHAR) :: name(*)
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
integer :: j, k, n
real(8) :: awr
integer, allocatable :: new_nuclide(:)
real(8), allocatable :: new_density(:)
character(:), allocatable :: name_
name_ = to_f_string(name)
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
! Check if nuclide is already in material
do j = 1, size(m % nuclide)
k = m % nuclide(j)
if (nuclides(k) % name == name_) then
awr = nuclides(k) % awr
m % density = m % density + density - m % atom_density(j)
m % density_gpcc = m % density_gpcc + (density - &
m % atom_density(j)) * awr * MASS_NEUTRON / N_AVOGADRO
m % atom_density(j) = density
err = 0
end if
end do
! If nuclide wasn't found, extend nuclide/density arrays
if (err /= 0) then
! If nuclide hasn't been loaded, load it now
err = openmc_load_nuclide(name)
if (err == 0) then
! Extend arrays
n = size(m % nuclide)
allocate(new_nuclide(n + 1))
new_nuclide(1:n) = m % nuclide
call move_alloc(FROM=new_nuclide, TO=m % nuclide)
allocate(new_density(n + 1))
new_density(1:n) = m % atom_density
call move_alloc(FROM=new_density, TO=m % atom_density)
! Append new nuclide/density
k = nuclide_dict % get_key(to_lower(name_))
m % nuclide(n + 1) = k
m % atom_density(n + 1) = density
m % density = m % density + density
m % density_gpcc = m % density_gpcc + &
density * nuclides(k) % awr * MASS_NEUTRON / N_AVOGADRO
m % n_nuclides = n + 1
end if
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_add_nuclide
!===============================================================================
! OPENMC_MATERIAL_GET_DENSITIES returns an array of nuclide densities in a
! material
!===============================================================================
function openmc_material_get_densities(index, nuclides, densities, n) &
result(err) bind(C)
integer(C_INT32_T), value :: index
type(C_PTR), intent(out) :: nuclides
type(C_PTR), intent(out) :: densities
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
if (allocated(m % atom_density)) then
nuclides = C_LOC(m % nuclide(1))
densities = C_LOC(m % atom_density(1))
n = size(m % atom_density)
err = 0
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_get_densities
!===============================================================================
! OPENMC_MATERIAL_GET_ID returns the ID of a material
!===============================================================================
function openmc_material_get_id(index, id) result(err) bind(C)
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(materials)) then
id = materials(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_get_id
!===============================================================================
! OPENMC_MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm
!===============================================================================
function openmc_material_set_density(index, density) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
err = m % set_density(density, nuclides)
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_set_density
!===============================================================================
! OPENMC_MATERIAL_SET_DENSITIES sets the densities for a list of nuclides in a
! material. If the nuclides don't already exist in the material, they will be
! added
!===============================================================================
function openmc_material_set_densities(index, n, name, density) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
type(C_PTR), intent(in) :: name(n)
real(C_DOUBLE), intent(in) :: density(n)
integer(C_INT) :: err
integer :: i
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: name_
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
! If nuclide/density arrays are not correct size, reallocate
if (n /= size(m % nuclide)) then
deallocate(m % nuclide, m % atom_density, m % p0)
allocate(m % nuclide(n), m % atom_density(n), m % p0(n))
end if
do i = 1, n
! Convert C string to Fortran string
call c_f_pointer(name(i), string, [10])
name_ = to_lower(to_f_string(string))
if (.not. nuclide_dict % has_key(name_)) then
err = openmc_load_nuclide(string)
if (err < 0) return
end if
m % nuclide(i) = nuclide_dict % get_key(name_)
m % atom_density(i) = density(i)
end do
m % n_nuclides = n
! Set isotropic flags to flags
m % p0(:) = .false.
! Set total density to the sum of the vector
err = m % set_density(sum(density), nuclides)
! Assign S(a,b) tables
call m % assign_sab_tables(nuclides, sab_tables)
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_set_densities
!===============================================================================
! OPENMC_NUCLIDE_NAME returns the name of a nuclide with a given index
!===============================================================================
function openmc_nuclide_name(index, name) result(err) bind(C)
integer(C_INT), value, intent(in) :: index
type(c_ptr), intent(out) :: name
integer(C_INT) :: err
character(C_CHAR), pointer :: name_
err = E_UNASSIGNED
if (allocated(nuclides)) then
if (index >= 1 .and. index <= size(nuclides)) then
name_ => nuclides(index) % name(1:1)
name = C_LOC(name_)
err = 0
else
err = E_OUT_OF_BOUNDS
end if
else
err = E_NUCLIDE_NOT_ALLOCATED
end if
end function openmc_nuclide_name
!===============================================================================
! OPENMC_RESET resets tallies and timers
!===============================================================================
subroutine openmc_reset() bind(C)
integer :: i
if (allocated(tallies)) then
do i = 1, size(tallies)
tallies(i) % n_realizations = 0
if (allocated(tallies(i) % results)) then
tallies(i) % results(:, :, :) = ZERO
end if
end do
end if
! Reset global tallies
n_realizations = 0
if (allocated(global_tallies)) then
global_tallies(:, :) = ZERO
end if
k_col_abs = ZERO
k_col_tra = ZERO
k_abs_tra = ZERO
k_sum(:) = ZERO
! Turn off tally flags
tallies_on = .false.
active_batches = .false.
! Clear active tally lists
call active_analog_tallies % clear()
call active_tracklength_tallies % clear()
call active_current_tallies % clear()
call active_collision_tallies % clear()
call active_tallies % clear()
! Reset timers
call time_total % reset()
call time_total % reset()
call time_initialize % reset()
call time_read_xs % reset()
call time_unionize % reset()
call time_bank % reset()
call time_bank_sample % reset()
call time_bank_sendrecv % reset()
call time_tallies % reset()
call time_inactive % reset()
call time_active % reset()
call time_transport % reset()
call time_finalize % reset()
end subroutine openmc_reset
!===============================================================================
! OPENMC_TALLY_GET_ID returns the ID of a tally
!===============================================================================
function openmc_tally_get_id(index, id) result(err) bind(C)
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(tallies)) then
id = tallies(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_tally_get_id
!===============================================================================
! OPENMC_TALLY_NUCLIDES returns the list of nuclides assigned to a tally
!===============================================================================
function openmc_tally_get_nuclides(index, nuclides, n) result(err) bind(C)
integer(C_INT32_T), value :: index
type(C_PTR), intent(out) :: nuclides
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(tallies)) then
associate (t => tallies(index))
if (allocated(t % nuclide_bins)) then
nuclides = C_LOC(t % nuclide_bins(1))
n = size(t % nuclide_bins)
err = 0
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_tally_get_nuclides
!===============================================================================
! OPENMC_TALLY_RESULTS returns a pointer to a tally results array along with its
! shape. This allows a user to obtain in-memory tally results from Python
! directly.
!===============================================================================
function openmc_tally_results(index, ptr, shape_) result(err) bind(C)
integer(C_INT32_T), intent(in), value :: index
type(C_PTR), intent(out) :: ptr
integer(C_INT), intent(out) :: shape_(3)
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(tallies)) then
if (allocated(tallies(index) % results)) then
ptr = C_LOC(tallies(index) % results(1,1,1))
shape_(:) = shape(tallies(index) % results)
err = 0
end if
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_tally_results
!===============================================================================
! OPENMC_TALLY_SET_NUCLIDES sets the nuclides in the tally which results should
! be scored for
!===============================================================================
function openmc_tally_set_nuclides(index, n, nuclides) result(err) bind(C)
integer(C_INT32_T), value :: index
integer(C_INT), value :: n
type(C_PTR), intent(in) :: nuclides(n)
integer(C_INT) :: err
integer :: i
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: nuclide_
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(tallies)) then
associate (t => tallies(index))
if (allocated(t % nuclide_bins)) deallocate(t % nuclide_bins)
allocate(t % nuclide_bins(n))
t % n_nuclide_bins = n
do i = 1, n
! Convert C string to Fortran string
call c_f_pointer(nuclides(i), string, [10])
nuclide_ = to_lower(to_f_string(string))
select case (nuclide_)
case ('total')
t % nuclide_bins(i) = -1
case default
if (nuclide_dict % has_key(nuclide_)) then
t % nuclide_bins(i) = nuclide_dict % get_key(nuclide_)
else
err = E_NUCLIDE_NOT_LOADED
return
end if
end select
end do
! Recalculate total number of scoring bins
t % total_score_bins = t % n_score_bins * t % n_nuclide_bins
! (Re)allocate results array
if (allocated(t % results)) deallocate(t % results)
allocate(t % results(3, t % total_score_bins, t % total_filter_bins))
t % results(:,:,:) = ZERO
err = 0
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_tally_set_nuclides
!===============================================================================
! TO_F_STRING takes a null-terminated array of C chars and turns it into a
! deferred-length character string. Yay Fortran 2003!
!===============================================================================
function to_f_string(c_string) result(f_string)
character(kind=C_CHAR), intent(in) :: c_string(*)
character(:), allocatable :: f_string
integer :: i, n
! Determine length of original string
n = 0
do while (c_string(n + 1) /= C_NULL_CHAR)
n = n + 1
end do
! Copy C string character by character
allocate(character(len=n) :: f_string)
do i = 1, n
f_string(i:i) = c_string(i)
end do
end function to_f_string
end module openmc_api

View file

@ -108,6 +108,9 @@ contains
real(8) :: hxyz(3) ! cell dimensions of current ijk cell
real(8) :: vol ! volume of cell
real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy
#ifdef MPI
integer :: mpi_err ! MPI error code
#endif
! Get maximum of spatial and group indices
nx = cmfd % indices(1)
@ -232,6 +235,9 @@ contains
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
type(RegularMesh), pointer :: m ! point to mesh
#ifdef MPI
integer :: mpi_err
#endif
! Associate pointer
m => meshes(n_user_meshes + 1)

View file

@ -2,7 +2,6 @@ module cross_section
use algorithm, only: binary_search
use constants
use energy_grid, only: log_spacing
use error, only: fatal_error
use global
use list_header, only: ListElemInt
@ -36,6 +35,7 @@ contains
integer :: i_grid ! index into logarithmic mapping array or material
! union grid
real(8) :: atom_density ! atom density of a nuclide
real(8) :: sab_frac ! fraction of atoms affected by S(a,b)
logical :: check_sab ! should we check for S(a,b) table?
! Set all material macroscopic cross sections to zero
@ -60,21 +60,25 @@ contains
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! ========================================================================
! ======================================================================
! CHECK FOR S(A,B) TABLE
i_sab = 0
sab_frac = ZERO
! Check if this nuclide matches one of the S(a,b) tables specified -- this
! relies on i_sab_nuclides being in sorted order
! Check if this nuclide matches one of the S(a,b) tables specified.
! This relies on i_sab_nuclides being in sorted order
if (check_sab) then
if (i == mat % i_sab_nuclides(j)) then
! Get index in sab_tables
i_sab = mat % i_sab_tables(j)
sab_frac = mat % sab_fracs(j)
! If particle energy is greater than the highest energy for the S(a,b)
! table, don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) i_sab = 0
! If particle energy is greater than the highest energy for the
! S(a,b) table, then don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
i_sab = 0
end if
! Increment position in i_sab_nuclides
j = j + 1
@ -84,7 +88,7 @@ contains
end if
end if
! ========================================================================
! ======================================================================
! CALCULATE MICROSCOPIC CROSS SECTION
! Determine microscopic cross sections for this nuclide
@ -92,13 +96,14 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
.or. i_sab /= micro_xs(i_nuclide) % index_sab &
.or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, &
p % sqrtkT, sab_frac)
end if
! ========================================================================
! ======================================================================
! ADD TO MACROSCOPIC CROSS SECTION
! Copy atom density of nuclide in material
@ -133,13 +138,15 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT, &
sab_frac)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent
real(8), intent(in) :: sqrtkT ! square root of kT, material dependent
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
logical :: use_mp ! true if XS can be calculated with windowed multipole
integer :: i_temp ! index for temperature
@ -243,8 +250,10 @@ contains
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % thermal = ZERO
micro_xs(i_nuclide) % thermal_elastic = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
@ -271,20 +280,21 @@ contains
end if
! Initialize sab treatment to false
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % elastic_sab = ZERO
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % sab_frac = ZERO
! Initialize URR probability table treatment to false
micro_xs(i_nuclide) % use_ptable = .false.
micro_xs(i_nuclide) % use_ptable = .false.
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
! If there is S(a,b) data for this nuclide, we need to set the sab_scatter
! and sab_elastic cross sections and correct the total and elastic cross
! sections.
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
if (i_sab > 0) then
call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac)
end if
! if the particle is in the unresolved resonance range and there are
! If the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. nuc % urr_present .and. .not. use_mp) then
@ -295,7 +305,6 @@ contains
end if
micro_xs(i_nuclide) % last_E = E
micro_xs(i_nuclide) % last_index_sab = i_sab
micro_xs(i_nuclide) % last_sqrtkT = sqrtkT
end associate
@ -303,16 +312,16 @@ contains
!===============================================================================
! CALCULATE_SAB_XS determines the elastic and inelastic scattering
! cross-sections in the thermal energy range. These cross sections replace
! whatever data were taken from the normal Nuclide table.
! cross-sections in the thermal energy range. These cross sections replace a
! fraction of whatever data were taken from the normal Nuclide table.
!===============================================================================
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
real(8), intent(in) :: sqrtkT ! temperature
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
integer :: i_grid ! index on S(a,b) energy grid
integer :: i_temp ! temperature index
@ -341,7 +350,8 @@ contains
! Randomly sample between temperature i and i+1
f = (kT - sab_tables(i_sab) % kTs(i_temp)) / &
(sab_tables(i_sab) % kTs(i_temp + 1) - sab_tables(i_sab) % kTs(i_temp))
(sab_tables(i_sab) % kTs(i_temp + 1) &
- sab_tables(i_sab) % kTs(i_temp))
if (f > prn()) i_temp = i_temp + 1
end if
@ -402,16 +412,20 @@ contains
end if
end associate
! Store the S(a,b) cross sections.
micro_xs(i_nuclide) % thermal = sab_frac * (elastic + inelastic)
micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic
! Correct total and elastic cross sections
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % elastic + inelastic + elastic
micro_xs(i_nuclide) % elastic = inelastic + elastic
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total &
+ micro_xs(i_nuclide) % thermal &
- sab_frac * micro_xs(i_nuclide) % elastic
micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % thermal &
+ (ONE - sab_frac) * micro_xs(i_nuclide) % elastic
! Store S(a,b) elastic cross section for sampling later
micro_xs(i_nuclide) % elastic_sab = elastic
! Save temperature index
! Save temperature index and thermal fraction
micro_xs(i_nuclide) % index_temp_sab = i_temp
micro_xs(i_nuclide) % sab_frac = sab_frac
end subroutine calculate_sab_xs

View file

@ -1,5 +1,6 @@
module eigenvalue
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use constants, only: ZERO
@ -39,6 +40,7 @@ contains
& temp_sites(:) ! local array of extra sites on each node
#ifdef MPI
integer :: mpi_err ! MPI error code
integer(8) :: n ! number of sites to send/recv
integer :: neighbor ! processor to send/recv data from
#ifdef MPIF08
@ -95,8 +97,7 @@ contains
! skip ahead in the sequence using the starting index in the 'global'
! fission bank for each processor.
call set_particle_seed(int((current_batch - 1)*gen_per_batch + &
current_gen,8))
call set_particle_seed(int(total_gen + overall_generation(), 8))
call advance_prn_seed(start)
! Determine how many fission sites we need to sample from the source bank
@ -371,20 +372,28 @@ contains
subroutine calculate_generation_keff()
integer :: i ! overall generation
#ifdef MPI
integer :: mpi_err ! MPI error code
#endif
! Get keff for this generation by subtracting off the starting value
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation
! Determine overall generation
i = overall_generation()
#ifdef MPI
! Combine values across all processors
call MPI_ALLREDUCE(keff_generation, k_generation(overall_gen), 1, &
MPI_REAL8, MPI_SUM, mpi_intracomm, mpi_err)
call MPI_ALLREDUCE(keff_generation, k_generation(i), 1, MPI_REAL8, &
MPI_SUM, mpi_intracomm, mpi_err)
#else
k_generation(overall_gen) = keff_generation
k_generation(i) = keff_generation
#endif
! Normalize single batch estimate of k
! TODO: This should be normalized by total_weight, not by n_particles
k_generation(overall_gen) = k_generation(overall_gen) / n_particles
k_generation(i) = k_generation(i) / n_particles
end subroutine calculate_generation_keff
@ -396,22 +405,24 @@ contains
subroutine calculate_average_keff()
integer :: i ! overall generation within simulation
integer :: n ! number of active generations
real(8) :: alpha ! significance level for CI
real(8) :: t_value ! t-value for confidence intervals
! Determine number of active generations
n = overall_gen - n_inactive*gen_per_batch
! Determine overall generation and number of active generations
i = overall_generation()
n = i - n_inactive*gen_per_batch
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next
! generation
keff = k_generation(overall_gen)
keff = k_generation(i)
else
! Sample mean of keff
k_sum(1) = k_sum(1) + k_generation(overall_gen)
k_sum(2) = k_sum(2) + k_generation(overall_gen)**2
k_sum(1) = k_sum(1) + k_generation(i)
k_sum(2) = k_sum(2) + k_generation(i)**2
! Determine mean
keff = k_sum(1) / n
@ -433,8 +444,8 @@ contains
end subroutine calculate_average_keff
!===============================================================================
! CALCULATE_COMBINED_KEFF calculates a minimum variance estimate of k-effective
! based on a linear combination of the collision, absorption, and tracklength
! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on
! a linear combination of the collision, absorption, and tracklength
! estimates. The theory behind this can be found in M. Halperin, "Almost
! linearly-optimum combination of unbiased estimates," J. Am. Stat. Assoc., 56,
! 36-43 (1961), doi:10.1080/01621459.1961.10482088. The implementation here
@ -442,7 +453,9 @@ contains
! of keff confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
!===============================================================================
subroutine calculate_combined_keff()
function openmc_get_keff(k_combined) result(err) bind(C)
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
integer :: l ! loop index
integer :: i, j, k ! indices referring to collision, absorption, or track
@ -453,9 +466,14 @@ contains
real(8) :: g ! sum of weighting factors
real(8) :: S(3) ! sums used for variance calculation
k_combined = ZERO
! Make sure we have at least four realizations. Notice that at the end,
! there is a N-3 term in a denominator.
if (n_realizations <= 3) return
if (n_realizations <= 3) then
err = -1
return
end if
! Initialize variables
n = real(n_realizations, 8)
@ -517,7 +535,6 @@ contains
! Initialize variables
g = ZERO
S = ZERO
k_combined = ZERO
do l = 1, 3
! Permutations of estimates
@ -584,8 +601,9 @@ contains
k_combined(2) = sqrt(k_combined(2))
end if
err = 0
end subroutine calculate_combined_keff
end function openmc_get_keff
!===============================================================================
! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and
@ -600,6 +618,7 @@ contains
logical :: sites_outside ! were there sites outside the ufs mesh?
#ifdef MPI
integer :: n ! total number of ufs mesh cells
integer :: mpi_err ! MPI error code
#endif
if (current_batch == 1 .and. current_gen == 1) then

View file

@ -1,67 +0,0 @@
module energy_grid
use global
use list_header, only: ListReal
use output, only: write_message
implicit none
integer :: grid_method ! how to treat the energy grid
integer :: n_log_bins ! number of bins for logarithmic grid
real(8) :: log_spacing ! spacing on logarithmic grid
contains
!===============================================================================
! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding
! indices on a nuclide energy grid
!===============================================================================
subroutine logarithmic_grid()
integer :: i, j, k ! Loop indices
integer :: t ! temperature index
integer :: M ! Number of equally log-spaced bins
real(8) :: E_max ! Maximum energy in MeV
real(8) :: E_min ! Minimum energy in MeV
real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid
! Set minimum/maximum energies
E_max = energy_max_neutron
E_min = energy_min_neutron
! Determine equal-logarithmic energy spacing
M = n_log_bins
log_spacing = log(E_max/E_min)/M
! Create equally log-spaced energy grid
allocate(umesh(0:M))
umesh(:) = [(i*log_spacing, i=0, M)]
do i = 1, n_nuclides_total
associate (nuc => nuclides(i))
do t = 1, size(nuc % grid)
! Allocate logarithmic mapping for nuclide
allocate(nuc % grid(t) % grid_index(0:M))
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M
do while (log(nuc % grid(t) % energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(nuc % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == size(nuc % grid(t) % energy)) exit
j = j + 1
end do
nuc % grid(t) % grid_index(k) = j
end do
end do
end associate
end do
deallocate(umesh)
end subroutine logarithmic_grid
end module energy_grid

View file

@ -1,41 +0,0 @@
module finalize
use hdf5, only: h5tclose_f, h5close_f
use global
use hdf5_interface, only: hdf5_bank_t
use message_passing
implicit none
contains
!===============================================================================
! FINALIZE_RUN does all post-simulation tasks such as calculating tally
! statistics and writing out tallies
!===============================================================================
subroutine openmc_finalize()
integer :: hdf5_err
! Deallocate arrays
call free_memory()
! Release compound datatypes
call h5tclose_f(hdf5_bank_t, hdf5_err)
! Close FORTRAN interface.
call h5close_f(hdf5_err)
#ifdef MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, mpi_err)
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
end subroutine openmc_finalize
end module finalize

View file

@ -198,11 +198,9 @@ contains
integer :: distribcell_index
integer :: i_xyz(3) ! indices in lattice
integer :: n ! number of cells to search
integer :: index_cell ! index in cells array
integer :: i_cell ! index in cells array
integer :: i_universe ! index in universes array
logical :: use_search_cells ! use cells provided as argument
type(Cell), pointer :: c ! pointer to cell
class(Lattice), pointer :: lat ! pointer to lattice
type(Universe), pointer :: univ ! universe to search in
do j = p % n_coord + 1, MAX_COORD
call p % coord(j) % reset()
@ -210,172 +208,174 @@ contains
j = p % n_coord
! set size of list to search
i_universe = p % coord(j) % universe
if (present(search_cells)) then
use_search_cells = .true.
n = size(search_cells)
else
use_search_cells = .false.
univ => universes(p % coord(j) % universe)
n = size(univ % cells)
n = size(universes(i_universe) % cells)
end if
found = .false.
CELL_LOOP: do i = 1, n
! select cells based on whether we are searching a universe or a provided
! list of cells (this would be for lists of neighbor cells)
if (use_search_cells) then
index_cell = search_cells(i)
i_cell = search_cells(i)
! check to make sure search cell is in same universe
if (cells(index_cell) % universe /= p % coord(j) % universe) cycle
if (cells(i_cell) % universe /= i_universe) cycle
else
index_cell = univ % cells(i)
i_cell = universes(i_universe) % cells(i)
end if
! get pointer to cell
c => cells(index_cell)
! Move on to the next cell if the particle is not inside this cell
if (.not. cell_contains(c, p)) cycle
if (cell_contains(cells(i_cell), p)) then
! Set cell on this level
p % coord(j) % cell = i_cell
! Set cell on this level
p % coord(j) % cell = index_cell
! Show cell information on trace
if (verbosity >= 10 .or. trace) then
call write_message(" Entering cell " // trim(to_str(&
cells(i_cell) % id)))
end if
! Show cell information on trace
if (verbosity >= 10 .or. trace) then
call write_message(" Entering cell " // trim(to_str(c % id)))
found = .true.
exit
end if
CELL_TYPE: if (c % type == FILL_MATERIAL) then
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
! Save previous material and temperature
p % last_material = p % material
p % last_sqrtkT = p % sqrtkT
! Get distributed offset
if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then
! Distributed instances of this cell have different
! materials/temperatures. Determine which instance this is for
! assigning the matching material/temperature.
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])) then
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z)
end if
end if
end do
end if
! Save the material
if (size(c % material) > 1) then
p % material = c % material(offset + 1)
else
p % material = c % material(1)
end if
! Save the temperature
if (size(c % sqrtkT) > 1) then
p % sqrtkT = c % sqrtkT(offset + 1)
else
p % sqrtkT = c % sqrtkT(1)
end if
elseif (c % type == FILL_UNIVERSE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
! Store lower level coordinates
p % coord(j + 1) % xyz = p % coord(j) % xyz
p % coord(j + 1) % uvw = p % coord(j) % uvw
! Move particle to next level and set universe
j = j + 1
p % n_coord = j
p % coord(j) % universe = c % fill
! Apply translation
if (allocated(c % translation)) then
p % coord(j) % xyz = p % coord(j) % xyz - c % translation
end if
! Apply rotation
if (allocated(c % rotation_matrix)) then
p % coord(j) % xyz = matmul(c % rotation_matrix, p % coord(j) % xyz)
p % coord(j) % uvw = matmul(c % rotation_matrix, p % coord(j) % uvw)
p % coord(j) % rotated = .true.
end if
call find_cell(p, found)
j = p % n_coord
if (.not. found) exit
elseif (c % type == FILL_LATTICE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
! Set current lattice
lat => lattices(c % fill) % obj
! Determine lattice indices
i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw)
! Store lower level coordinates
p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz)
p % coord(j + 1) % uvw = p % coord(j) % uvw
! set particle lattice indices
p % coord(j + 1) % lattice = c % fill
p % coord(j + 1) % lattice_x = i_xyz(1)
p % coord(j + 1) % lattice_y = i_xyz(2)
p % coord(j + 1) % lattice_z = i_xyz(3)
! Set the next lowest coordinate level.
if (lat % are_valid_indices(i_xyz)) then
! Particle is inside the lattice.
p % coord(j + 1) % universe = &
lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
else
! Particle is outside the lattice.
if (lat % outer == NO_OUTER_UNIVERSE) then
call handle_lost_particle(p, "Particle " // trim(to_str(p %id)) &
// " is outside lattice " // trim(to_str(lat % id)) &
// " but the lattice has no defined outer universe.")
return
else
p % coord(j + 1) % universe = lat % outer
end if
end if
! Move particle to next level and search for the lower cells.
j = j + 1
p % n_coord = j
call find_cell(p, found)
j = p % n_coord
if (.not. found) exit
end if CELL_TYPE
! Found cell so we can return
found = .true.
return
end do CELL_LOOP
found = .false.
if (found) then
associate(c => cells(i_cell))
CELL_TYPE: if (c % type == FILL_MATERIAL) then
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
! Save previous material and temperature
p % last_material = p % material
p % last_sqrtkT = p % sqrtkT
! Get distributed offset
if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then
! Distributed instances of this cell have different
! materials/temperatures. Determine which instance this is for
! assigning the matching material/temperature.
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])) then
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z)
end if
end if
end do
! Keep track of which instance of the cell the particle is in
p % cell_instance = offset + 1
else
p % cell_instance = 1
end if
! Save the material
if (size(c % material) > 1) then
p % material = c % material(offset + 1)
else
p % material = c % material(1)
end if
! Save the temperature
if (size(c % sqrtkT) > 1) then
p % sqrtkT = c % sqrtkT(offset + 1)
else
p % sqrtkT = c % sqrtkT(1)
end if
elseif (c % type == FILL_UNIVERSE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
! Store lower level coordinates
p % coord(j + 1) % xyz = p % coord(j) % xyz
p % coord(j + 1) % uvw = p % coord(j) % uvw
! Move particle to next level and set universe
j = j + 1
p % n_coord = j
p % coord(j) % universe = c % fill
! Apply translation
if (allocated(c % translation)) then
p % coord(j) % xyz = p % coord(j) % xyz - c % translation
end if
! Apply rotation
if (allocated(c % rotation_matrix)) then
p % coord(j) % xyz = matmul(c % rotation_matrix, p % coord(j) % xyz)
p % coord(j) % uvw = matmul(c % rotation_matrix, p % coord(j) % uvw)
p % coord(j) % rotated = .true.
end if
call find_cell(p, found)
j = p % n_coord
elseif (c % type == FILL_LATTICE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
associate (lat => lattices(c % fill) % obj)
! Determine lattice indices
i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw)
! Store lower level coordinates
p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz)
p % coord(j + 1) % uvw = p % coord(j) % uvw
! set particle lattice indices
p % coord(j + 1) % lattice = c % fill
p % coord(j + 1) % lattice_x = i_xyz(1)
p % coord(j + 1) % lattice_y = i_xyz(2)
p % coord(j + 1) % lattice_z = i_xyz(3)
! Set the next lowest coordinate level.
if (lat % are_valid_indices(i_xyz)) then
! Particle is inside the lattice.
p % coord(j + 1) % universe = &
lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
else
! Particle is outside the lattice.
if (lat % outer == NO_OUTER_UNIVERSE) then
call handle_lost_particle(p, "Particle " // trim(to_str(p %id)) &
// " is outside lattice " // trim(to_str(lat % id)) &
// " but the lattice has no defined outer universe.")
return
else
p % coord(j + 1) % universe = lat % outer
end if
end if
end associate
! Move particle to next level and search for the lower cells.
j = j + 1
p % n_coord = j
call find_cell(p, found)
j = p % n_coord
end if CELL_TYPE
end associate
end if
end subroutine find_cell

View file

@ -43,11 +43,11 @@ module global
type(VolumeCalculation), allocatable :: volume_calcs(:)
! Size of main arrays
integer :: n_cells ! # of cells
integer(C_INT32_T), bind(C) :: n_cells ! # of cells
integer :: n_universes ! # of universes
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
integer(C_INT32_T), bind(C) :: n_materials ! # of materials
integer :: n_plots ! # of plots
! These dictionaries provide a fast lookup mechanism -- the key is the
@ -73,7 +73,8 @@ module global
! ============================================================================
! CROSS SECTION RELATED VARIABLES NEEDED REGARDLESS OF CE OR MG
integer :: n_nuclides_total ! Number of nuclide cross section tables
! Number of nuclide cross section tables
integer(C_INT), bind(C, name='n_nuclides') :: n_nuclides_total
! Cross section caches
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
@ -81,6 +82,10 @@ module global
! Dictionaries to look up cross sections and listings
type(DictCharInt) :: nuclide_dict
type(DictCharInt) :: library_dict
! Cross section libraries
type(Library), allocatable :: libraries(:)
! ============================================================================
! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES
@ -106,6 +111,10 @@ module global
logical :: temperature_multipole = .false.
real(8) :: temperature_tolerance = 10.0_8
real(8) :: temperature_default = 293.6_8
real(8) :: temperature_range(2) = [ZERO, ZERO]
integer :: n_log_bins ! number of bins for logarithmic grid
real(8) :: log_spacing ! spacing on logarithmic grid
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
@ -132,7 +141,7 @@ module global
integer :: max_order
! Whether or not to convert Legendres to tabulars
logical :: legendre_to_tabular = .True.
logical :: legendre_to_tabular = .true.
! Number of points to use in the Legendre to tabular conversion
integer :: legendre_to_tabular_points = 33
@ -185,7 +194,7 @@ module global
integer :: n_user_meshes = 0 ! # of structured user meshes
integer :: n_filters = 0 ! # of filters
integer :: n_user_filters = 0 ! # of user filters
integer :: n_tallies = 0 ! # of tallies
integer(C_INT32_T), bind(C) :: n_tallies = 0 ! # of tallies
integer :: n_user_tallies = 0 ! # of user tallies
! Tally derivatives
@ -214,9 +223,9 @@ module global
integer :: n_inactive ! # of inactive batches
integer :: n_active ! # of active batches
integer :: gen_per_batch = 1 ! # of generations per batch
integer :: current_batch = 0 ! current batch
integer :: current_gen = 0 ! current generation within a batch
integer :: overall_gen = 0 ! overall generation in the run
integer :: current_batch ! current batch
integer :: current_gen ! current generation within a batch
integer :: total_gen = 0 ! total number of generations simulated
! ============================================================================
! TALLY PRECISION TRIGGER VARIABLES
@ -249,7 +258,6 @@ module global
real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
real(8) :: k_combined(2) ! combined best estimate of k-effective
! Shannon entropy
logical :: entropy_on = .false.
@ -460,6 +468,7 @@ contains
if (allocated(surfaces)) deallocate(surfaces)
if (allocated(materials)) deallocate(materials)
if (allocated(plots)) deallocate(plots)
if (allocated(volume_calcs)) deallocate(volume_calcs)
! Deallocate geometry debugging information
if (allocated(overlap_check_cnt)) deallocate(overlap_check_cnt)
@ -472,6 +481,7 @@ contains
end do
deallocate(nuclides)
end if
if (allocated(libraries)) deallocate(libraries)
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
@ -480,7 +490,6 @@ contains
if (allocated(macro_xs)) deallocate(macro_xs)
if (allocated(sab_tables)) deallocate(sab_tables)
if (allocated(micro_xs)) deallocate(micro_xs)
! Deallocate external source
if (allocated(external_source)) deallocate(external_source)
@ -495,21 +504,24 @@ contains
if (allocated(meshes)) deallocate(meshes)
if (allocated(filters)) deallocate(filters)
if (allocated(tallies)) deallocate(tallies)
if (allocated(filter_matches)) deallocate(filter_matches)
! Deallocate fission and source bank and entropy
!$omp parallel
if (allocated(fission_bank)) deallocate(fission_bank)
if (allocated(filter_matches)) deallocate(filter_matches)
if (allocated(tally_derivs)) deallocate(tally_derivs)
!$omp end parallel
#ifdef _OPENMP
if (allocated(master_fission_bank)) deallocate(master_fission_bank)
#endif
if (allocated(source_bank)) deallocate(source_bank)
if (allocated(entropy_p)) deallocate(entropy_p)
! Deallocate array of work indices
if (allocated(work_index)) deallocate(work_index)
if (allocated(energy_bins)) deallocate(energy_bins)
if (allocated(energy_bin_avg)) deallocate(energy_bin_avg)
! Deallocate CMFD
call deallocate_cmfd(cmfd)
@ -536,6 +548,7 @@ contains
call plot_dict % clear()
call nuclide_dict % clear()
call sab_dict % clear()
call library_dict % clear()
! Clear statepoint and sourcepoint batch set
call statepoint_batch % clear()
@ -563,4 +576,13 @@ contains
end subroutine free_memory
!===============================================================================
! OVERALL_GENERATION determines the overall generation number
!===============================================================================
pure function overall_generation() result(gen)
integer :: gen
gen = gen_per_batch*(current_batch - 1) + current_gen
end function overall_generation
end module global

View file

@ -11,7 +11,6 @@ module initialize
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
@ -46,11 +45,28 @@ contains
! setting up timers, etc.
!===============================================================================
subroutine openmc_init(intracomm)
#ifdef MPIF08
type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator
#else
subroutine openmc_init(intracomm) bind(C)
integer, intent(in), optional :: intracomm ! MPI intracommunicator
! Copy the communicator to a new variable. This is done to avoid changing
! the signature of this subroutine. If MPI is being used but no communicator
! was passed, assume MPI_COMM_WORLD.
#ifdef MPI
#ifdef MPIF08
type(MPI_Comm), intent(in) :: comm ! MPI intracommunicator
if (present(intracomm)) then
comm % MPI_VAL = intracomm
else
comm = MPI_COMM_WORLD
end if
#else
integer :: comm
if (present(intracomm)) then
comm = intracomm
else
comm = MPI_COMM_WORLD
end if
#endif
#endif
! Start total and initialization timer
@ -59,7 +75,7 @@ contains
#ifdef MPI
! Setup MPI
call initialize_mpi(intracomm)
call initialize_mpi(comm)
#endif
! Initialize HDF5 interface
@ -101,12 +117,6 @@ contains
end if
if (run_mode /= MODE_PLOTTING) then
! Construct information needed for nuclear data
if (run_CE) then
! Construct log energy grid for cross-sections
call logarithmic_grid()
end if
! Allocate and setup tally stride, filter_matches, and tally maps
call configure_tallies()
@ -162,6 +172,7 @@ contains
integer, intent(in) :: intracomm ! MPI intracommunicator
#endif
integer :: mpi_err ! MPI error code
integer :: bank_blocks(5) ! Count for each datatype
#ifdef MPIF08
type(MPI_Datatype) :: bank_types(5)

View file

@ -7,7 +7,6 @@ module input_xml
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, &
get_temperatures, root_universe
@ -951,6 +950,9 @@ contains
if (check_for_node(root, "temperature_multipole")) then
call get_node_value(root, "temperature_multipole", temperature_multipole)
end if
if (check_for_node(root, "temperature_range")) then
call get_node_array(root, "temperature_range", temperature_range)
end if
! Check for tabular_legendre options
if (check_for_node(root, "tabular_legendre")) then
@ -2068,8 +2070,6 @@ contains
subroutine read_materials()
integer :: i, j
type(DictCharInt) :: library_dict
type(Library), allocatable :: libraries(:)
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
real(8), allocatable :: material_temps(:)
@ -2158,9 +2158,9 @@ contains
! Now that the cross_sections.xml or mgxs.h5 has been located, read it in
if (run_CE) then
call read_ce_cross_sections_xml(libraries)
call read_ce_cross_sections_xml()
else
call read_mg_cross_sections_header(libraries)
call read_mg_cross_sections_header()
end if
! Creating dictionary that maps the name of the material to the entry
@ -2181,7 +2181,7 @@ contains
end if
! Parse data from materials.xml
call read_materials_xml(libraries, library_dict, material_temps)
call read_materials_xml(material_temps)
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures(material_temps)
@ -2194,17 +2194,12 @@ contains
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
call time_read_xs % start()
call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
call read_ce_cross_sections(nuc_temps, sab_temps)
call time_read_xs % stop()
end if
! Clear dictionary
call library_dict % clear()
end subroutine read_materials
subroutine read_materials_xml(libraries, library_dict, material_temps)
type(Library), intent(in) :: libraries(:)
type(DictCharInt), intent(inout) :: library_dict
subroutine read_materials_xml(material_temps)
real(8), allocatable, intent(out) :: material_temps(:)
integer :: i ! loop index for materials
@ -2565,6 +2560,7 @@ contains
! table is indeed applied to multiple nuclides.
allocate(mat % sab_names(n_sab))
allocate(mat % i_sab_tables(n_sab))
allocate(mat % sab_fracs(n_sab))
do j = 1, n_sab
! Get pointer to S(a,b) table
@ -2578,6 +2574,13 @@ contains
name = trim(name)
mat % sab_names(j) = name
! Read the fraction of nuclei affected by this S(a,b) table
if (check_for_node(node_sab, "fraction")) then
call get_node_value(node_sab, "fraction", mat % sab_fracs(j))
else
mat % sab_fracs(j) = ONE
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. library_dict % has_key(to_lower(name))) then
call fatal_error("Could not find S(a,b) table " // trim(name) &
@ -3372,6 +3375,15 @@ contains
! =======================================================================
! READ DATA FOR FILTERS
! Check if user is using old XML format and throw an error if so
if (check_for_node(node_tal, "filter")) then
call fatal_error("Tally filters must be specified independently of &
&tallies in a <filter> element. The <tally> element itself should &
&have a list of filters that apply, e.g., <filters>1 2</filters> &
&where 1 and 2 are the IDs of filters specified outside of &
&<tally>.")
end if
! Determine number of filters
if (check_for_node(node_tal, "filters")) then
n_filter = node_word_count(node_tal, "filters")
@ -4813,9 +4825,7 @@ contains
! file contains a listing of the CE and MG cross sections that may be used.
!===============================================================================
subroutine read_ce_cross_sections_xml(libraries)
type(Library), allocatable, intent(out) :: libraries(:)
subroutine read_ce_cross_sections_xml()
integer :: i ! loop index
integer :: n
integer :: n_libraries
@ -4918,9 +4928,7 @@ contains
end subroutine read_ce_cross_sections_xml
subroutine read_mg_cross_sections_header(libraries)
type(Library), allocatable, intent(out) :: libraries(:)
subroutine read_mg_cross_sections_header()
integer :: i ! loop index
integer :: n_libraries
logical :: file_exists ! does mgxs.h5 exist?
@ -5161,103 +5169,7 @@ contains
end subroutine normalize_ao
!===============================================================================
! ASSIGN_SAB_TABLES assigns S(alpha,beta) tables to specific nuclides within
! materials so the code knows when to apply bound thermal scattering data
!===============================================================================
subroutine assign_sab_tables()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: k ! index over S(a,b) tables in material
integer :: m ! position for sorting
integer :: temp_nuclide ! temporary value for sorting
integer :: temp_table ! temporary value for sorting
logical :: found
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
if (.not. allocated(materials(i) % i_sab_tables)) cycle
associate (mat => materials(i))
ASSIGN_SAB: do k = 1, size(mat % i_sab_tables)
! In order to know which nuclide the S(a,b) table applies to, we need
! to search through the list of nuclides for one which has a matching
! name
found = .false.
associate (sab => sab_tables(mat % i_sab_tables(k)))
FIND_NUCLIDE: do j = 1, size(mat % nuclide)
if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then
call i_sab_tables % push_back(mat % i_sab_tables(k))
call i_sab_nuclides % push_back(j)
found = .true.
end if
end do FIND_NUCLIDE
end associate
! Check to make sure S(a,b) table matched a nuclide
if (.not. found) then
call fatal_error("S(a,b) table " // trim(mat % &
sab_names(k)) // " did not match any nuclide on material " &
// trim(to_str(mat % id)))
end if
end do ASSIGN_SAB
! Update i_sab_tables and i_sab_nuclides
deallocate(mat % i_sab_tables)
m = i_sab_tables % size()
allocate(mat % i_sab_tables(m))
allocate(mat % i_sab_nuclides(m))
mat % i_sab_tables(:) = i_sab_tables % data(1:m)
mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m)
! Clear entries in vectors for next material
call i_sab_tables % clear()
call i_sab_nuclides % clear()
! If there are multiple S(a,b) tables, we need to make sure that the
! entries in i_sab_nuclides are sorted or else they won't be applied
! correctly in the cross_section module. The algorithm here is a simple
! insertion sort -- don't need anything fancy!
if (size(mat % i_sab_tables) > 1) then
SORT_SAB: do k = 2, size(mat % i_sab_tables)
! Save value to move
m = k
temp_nuclide = mat % i_sab_nuclides(k)
temp_table = mat % i_sab_tables(k)
MOVE_OVER: do
! Check if insertion value is greater than (m-1)th value
if (temp_nuclide >= mat % i_sab_nuclides(m-1)) exit
! Move values over until hitting one that's not larger
mat % i_sab_nuclides(m) = mat % i_sab_nuclides(m-1)
mat % i_sab_tables(m) = mat % i_sab_tables(m-1)
m = m - 1
! Exit if we've reached the beginning of the list
if (m == 1) exit
end do MOVE_OVER
! Put the original value into its new position
mat % i_sab_nuclides(m) = temp_nuclide
mat % i_sab_tables(m) = temp_table
end do SORT_SAB
end if
! Deallocate temporary arrays for names of nuclides and S(a,b) tables
if (allocated(mat % names)) deallocate(mat % names)
end associate
end do
end subroutine assign_sab_tables
subroutine read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
type(Library), intent(in) :: libraries(:)
type(DictCharInt), intent(inout) :: library_dict
subroutine read_ce_cross_sections(nuc_temps, sab_temps)
type(VectorReal), intent(in) :: nuc_temps(:)
type(VectorReal), intent(in) :: sab_temps(:)
@ -5273,9 +5185,6 @@ contains
allocate(nuclides(n_nuclides_total))
allocate(sab_tables(n_sab_tables))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! Read cross sections
do i = 1, size(materials)
@ -5296,7 +5205,8 @@ contains
! Read nuclide data from HDF5
group_id = open_group(file_id, name)
call nuclides(i_nuclide) % from_hdf5(group_id, nuc_temps(i_nuclide), &
temperature_method, temperature_tolerance, master)
temperature_method, temperature_tolerance, temperature_range, &
master)
call close_group(group_id)
call file_close(file_id)
@ -5326,6 +5236,13 @@ contains
end do
end do
! Set up logarithmic grid for nuclides
do i = 1, size(nuclides)
call nuclides(i) % init_grid(energy_min_neutron, &
energy_max_neutron, n_log_bins)
end do
log_spacing = log(energy_max_neutron/energy_min_neutron) / n_log_bins
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
if (.not. allocated(materials(i) % sab_names)) cycle
@ -5348,7 +5265,7 @@ contains
! Read S(a,b) data from HDF5
group_id = open_group(file_id, name)
call sab_tables(i_sab) % from_hdf5(group_id, sab_temps(i_sab), &
temperature_method, temperature_tolerance)
temperature_method, temperature_tolerance, temperature_range)
call close_group(group_id)
call file_close(file_id)
@ -5356,10 +5273,10 @@ contains
call already_read % add(name)
end if
end do
end do
! Associate S(a,b) tables with specific nuclides
call assign_sab_tables()
! Associate S(a,b) tables with specific nuclides
call materials(i) % assign_sab_tables(nuclides, sab_tables)
end do
! Show which nuclide results in lowest energy for neutron transport
do i = 1, size(nuclides)

View file

@ -1,20 +1,25 @@
program main
use constants
use finalize, only: openmc_finalize
use global
use initialize, only: openmc_init
use message_passing
use particle_restart, only: run_particle_restart
use plot, only: run_plot
use simulation, only: run_simulation
use volume_calc, only: run_volume_calculations
use openmc_api, only: openmc_init, openmc_finalize, openmc_run, &
openmc_plot_geometry, openmc_calculate_volumes
use particle_restart, only: run_particle_restart
implicit none
#ifdef MPI
integer :: mpi_err ! MPI error code
#endif
! Initialize run -- when run with MPI, pass communicator
#ifdef MPI
#ifdef MPIF08
call openmc_init(MPI_COMM_WORLD % MPI_VAL)
#else
call openmc_init(MPI_COMM_WORLD)
#endif
#else
call openmc_init()
#endif
@ -22,16 +27,22 @@ program main
! start problem based on mode
select case (run_mode)
case (MODE_FIXEDSOURCE, MODE_EIGENVALUE)
call run_simulation()
call openmc_run()
case (MODE_PLOTTING)
call run_plot()
call openmc_plot_geometry()
case (MODE_PARTICLE)
if (master) call run_particle_restart()
case (MODE_VOLUME)
call run_volume_calculations()
call openmc_calculate_volumes()
end select
! finalize run
call openmc_finalize()
#ifdef MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
end program main

View file

@ -1,5 +1,12 @@
module material_header
use constants
use error, only: fatal_error
use nuclide_header, only: Nuclide
use sab_header, only: SAlphaBeta
use stl_vector, only: VectorReal, VectorInt
use string, only: to_str
implicit none
!===============================================================================
@ -22,10 +29,11 @@ module material_header
! Unionized energy grid information
integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers
! S(a,b) data references
! S(a,b) data
integer :: n_sab = 0 ! number of S(a,b) tables
integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide
integer, allocatable :: i_sab_tables(:) ! index in sab_tables
real(8), allocatable :: sab_fracs(:) ! how often to use S(a,b)
! Temporary names read during initialization
character(20), allocatable :: names(:) ! isotope names
@ -38,6 +46,163 @@ module material_header
! enforce isotropic scattering in lab
logical, allocatable :: p0(:)
contains
procedure :: set_density => material_set_density
procedure :: assign_sab_tables => material_assign_sab_tables
end type Material
contains
!===============================================================================
! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm.
!===============================================================================
function material_set_density(m, density, nuclides) result(err)
class(Material), intent(inout) :: m
real(8), intent(in) :: density
type(Nuclide), intent(in) :: nuclides(:)
integer :: err
integer :: i
real(8) :: sum_percent
real(8) :: awr
err = -1
if (allocated(m % atom_density)) then
! Set total density based on value provided
m % density = density
! Determine normalized atom percents
sum_percent = sum(m % atom_density)
m % atom_density(:) = m % atom_density / sum_percent
! Recalculate nuclide atom densities based on given density
m % atom_density(:) = density * m % atom_density
! Calculate density in g/cm^3.
m % density_gpcc = ZERO
do i = 1, m % n_nuclides
awr = nuclides(m % nuclide(i)) % awr
m % density_gpcc = m % density_gpcc &
+ m % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
end do
err = 0
end if
end function material_set_density
!===============================================================================
! ASSIGN_SAB_TABLES assigns S(alpha,beta) tables to specific nuclides within
! materials so the code knows when to apply bound thermal scattering data
!===============================================================================
subroutine material_assign_sab_tables(this, nuclides, sab_tables)
class(Material), intent(inout) :: this
type(Nuclide), intent(in) :: nuclides(:)
type(SAlphaBeta), intent(in) :: sab_tables(:)
integer :: j ! index over nuclides in material
integer :: k ! index over S(a,b) tables in material
integer :: m ! position for sorting
integer :: temp_nuclide ! temporary value for sorting
integer :: temp_table ! temporary value for sorting
real(8) :: temp_frac ! temporary value for sorting
logical :: found
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
type(VectorReal) :: sab_fracs
if (.not. allocated(this % i_sab_tables)) return
ASSIGN_SAB: do k = 1, size(this % i_sab_tables)
! In order to know which nuclide the S(a,b) table applies to, we need
! to search through the list of nuclides for one which has a matching
! name
found = .false.
associate (sab => sab_tables(this % i_sab_tables(k)))
FIND_NUCLIDE: do j = 1, size(this % nuclide)
if (any(sab % nuclides == nuclides(this % nuclide(j)) % name)) then
call i_sab_tables % push_back(this % i_sab_tables(k))
call i_sab_nuclides % push_back(j)
call sab_fracs % push_back(this % sab_fracs(k))
found = .true.
end if
end do FIND_NUCLIDE
end associate
! Check to make sure S(a,b) table matched a nuclide
if (.not. found) then
call fatal_error("S(a,b) table " // trim(this % &
sab_names(k)) // " did not match any nuclide on material " &
// trim(to_str(this % id)))
end if
end do ASSIGN_SAB
! Make sure each nuclide only appears in one table.
do j = 1, i_sab_nuclides % size()
do k = j+1, i_sab_nuclides % size()
if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then
call fatal_error(trim( &
nuclides(this % nuclide(i_sab_nuclides % data(j))) % name) &
// " in material " // trim(to_str(this % id)) // " was found &
&in multiple S(a,b) tables. Each nuclide can only appear in &
&one S(a,b) table per material.")
end if
end do
end do
! Update i_sab_tables and i_sab_nuclides
deallocate(this % i_sab_tables)
deallocate(this % sab_fracs)
if (allocated(this % i_sab_nuclides)) deallocate(this % i_sab_nuclides)
m = i_sab_tables % size()
allocate(this % i_sab_tables(m))
allocate(this % i_sab_nuclides(m))
allocate(this % sab_fracs(m))
this % i_sab_tables(:) = i_sab_tables % data(1:m)
this % i_sab_nuclides(:) = i_sab_nuclides % data(1:m)
this % sab_fracs(:) = sab_fracs % data(1:m)
! Clear entries in vectors for next material
call i_sab_tables % clear()
call i_sab_nuclides % clear()
call sab_fracs % clear()
! If there are multiple S(a,b) tables, we need to make sure that the
! entries in i_sab_nuclides are sorted or else they won't be applied
! correctly in the cross_section module. The algorithm here is a simple
! insertion sort -- don't need anything fancy!
if (size(this % i_sab_tables) > 1) then
SORT_SAB: do k = 2, size(this % i_sab_tables)
! Save value to move
m = k
temp_nuclide = this % i_sab_nuclides(k)
temp_table = this % i_sab_tables(k)
temp_frac = this % i_sab_tables(k)
MOVE_OVER: do
! Check if insertion value is greater than (m-1)th value
if (temp_nuclide >= this % i_sab_nuclides(m-1)) exit
! Move values over until hitting one that's not larger
this % i_sab_nuclides(m) = this % i_sab_nuclides(m-1)
this % i_sab_tables(m) = this % i_sab_tables(m-1)
this % sab_fracs(m) = this % sab_fracs(m-1)
m = m - 1
! Exit if we've reached the beginning of the list
if (m == 1) exit
end do MOVE_OVER
! Put the original value into its new position
this % i_sab_nuclides(m) = temp_nuclide
this % i_sab_tables(m) = temp_table
this % sab_fracs(m) = temp_frac
end do SORT_SAB
end if
! Deallocate temporary arrays for names of nuclides and S(a,b) tables
if (allocated(this % names)) deallocate(this % names)
end subroutine material_assign_sab_tables
end module material_header

View file

@ -137,21 +137,22 @@ contains
real(8), intent(in), optional :: energies(:) ! energy grid to search
integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc)
logical, intent(inout), optional :: sites_outside ! were there sites outside mesh?
real(8), allocatable :: cnt_(:,:,:,:)
integer :: i ! loop index for local fission sites
integer :: n_sites ! size of bank array
integer :: ijk(3) ! indices on mesh
integer :: n_groups ! number of groups in energies
integer :: n ! number of energy groups / size
integer :: e_bin ! energy_bin
#ifdef MPI
integer :: mpi_err ! MPI error code
#endif
logical :: in_mesh ! was single site outside mesh?
logical :: outside ! was any site outside mesh?
#ifdef MPI
integer :: n ! total size of count variable
real(8) :: dummy ! temporary receive buffer for non-root reductions
#endif
! initialize variables
cnt = ZERO
allocate(cnt_(size(cnt,1), size(cnt,2), size(cnt,3), size(cnt,4)))
cnt_ = ZERO
outside = .false.
! Set size of bank
@ -163,9 +164,9 @@ contains
! Determine number of energies in group structure
if (present(energies)) then
n_groups = size(energies) - 1
n = size(energies) - 1
else
n_groups = 1
n = 1
end if
! loop over fission sites and count how many are in each mesh box
@ -183,42 +184,33 @@ contains
if (present(energies)) then
if (bank_array(i) % E < energies(1)) then
e_bin = 1
elseif (bank_array(i) % E > energies(n_groups+1)) then
e_bin = n_groups
elseif (bank_array(i) % E > energies(n + 1)) then
e_bin = n
else
e_bin = binary_search(energies, n_groups + 1, bank_array(i) % E)
e_bin = binary_search(energies, n + 1, bank_array(i) % E)
end if
else
e_bin = 1
end if
! add to appropriate mesh box
cnt(e_bin,ijk(1),ijk(2),ijk(3)) = cnt(e_bin,ijk(1),ijk(2),ijk(3)) + &
cnt_(e_bin,ijk(1),ijk(2),ijk(3)) = cnt_(e_bin,ijk(1),ijk(2),ijk(3)) + &
bank_array(i) % wgt
end do FISSION_SITES
#ifdef MPI
! determine total number of mesh cells
n = n_groups * size(cnt,2) * size(cnt,3) * size(cnt,4)
! collect values from all processors
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, cnt, n, MPI_REAL8, MPI_SUM, 0, &
mpi_intracomm, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(cnt, dummy, n, MPI_REAL8, MPI_SUM, 0, &
mpi_intracomm, mpi_err)
end if
n = size(cnt_)
call MPI_REDUCE(cnt_, cnt, n, MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
! Check if there were sites outside the mesh for any processor
if (present(sites_outside)) then
call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, &
mpi_intracomm, mpi_err)
end if
#else
sites_outside = outside
cnt = cnt_
#endif
end subroutine count_bank_sites

View file

@ -16,7 +16,6 @@ module message_passing
integer :: rank = 0 ! rank of process
logical :: master = .true. ! master process?
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
integer :: mpi_err ! MPI error code
#ifdef MPIF08
type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank
type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator

View file

@ -73,9 +73,6 @@ contains
! allocate arrays for MGXS storage and cross section cache
allocate(nuclides_MG(n_nuclides_total))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! ==========================================================================
! READ ALL MGXS CROSS SECTION TABLES

View file

@ -97,38 +97,41 @@ module nuclide_header
contains
procedure :: clear => nuclide_clear
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
procedure :: nu => nuclide_nu
procedure, private :: create_derived => nuclide_create_derived
end type Nuclide
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a
! particular nuclide at the current energy
! NUCLIDEMICROXS contains cached microscopic cross sections for a particular
! nuclide at the current energy
!===============================================================================
type NuclideMicroXS
integer :: index_grid ! index on nuclide energy grid
integer :: index_temp ! temperature index for nuclide
real(8) :: last_E = ZERO ! last evaluated energy
real(8) :: interp_factor ! interpolation factor on nuc. energy grid
real(8) :: total ! microscropic total xs
real(8) :: elastic ! microscopic elastic scattering xs
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
! Microscopic cross sections in barns
real(8) :: total
real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is
! averaged over bound and non-bound nuclei
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
real(8) :: thermal ! Bound thermal elastic & inelastic scattering
real(8) :: thermal_elastic ! Bound thermal elastic scattering
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide
integer :: index_temp_sab ! temperature index for sab_tables
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
! Indicies and factors needed to compute cross sections from the data tables
integer :: index_grid ! Index on nuclide energy grid
integer :: index_temp ! Temperature index for nuclide
real(8) :: interp_factor ! Interpolation factor on nuc. energy grid
integer :: index_sab = NONE ! Index in sab_tables
integer :: index_temp_sab ! Temperature index for sab_tables
real(8) :: sab_frac ! Fraction of atoms affected by S(a,b)
logical :: use_ptable ! In URR range with probability tables?
! Information for URR probability table use
logical :: use_ptable ! in URR range with probability tables?
! Information for Doppler broadening
! Energy and temperature last used to evaluate these cross sections. If
! these values have changed, then the cross sections must be re-evaluated.
real(8) :: last_E = ZERO ! Last evaluated energy
real(8) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann
! constant * temperature (eV))
! constant * temperature (eV))
end type NuclideMicroXS
!===============================================================================
@ -168,12 +171,13 @@ module nuclide_header
end subroutine nuclide_clear
subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance, &
master)
minmax, master)
class(Nuclide), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
type(VectorReal), intent(in) :: temperature ! list of desired temperatures
type(VectorReal), intent(in) :: temperature ! list of desired temperatures
integer, intent(inout) :: method
real(8), intent(in) :: tolerance
real(8), intent(in) :: minmax(2) ! range of temperatures
logical, intent(in) :: master ! if this is the master proc
integer :: i
@ -233,7 +237,18 @@ module nuclide_header
method = TEMPERATURE_NEAREST
end if
! Determine actual temperatures to read
! Determine actual temperatures to read -- start by checking whether a
! temperature range was given, in which case all temperatures in the range
! are loaded irrespective of what temperatures actually appear in the model
if (minmax(2) > ZERO) then
do i = 1, size(temps_available)
temp_actual = temps_available(i)
if (minmax(1) <= temp_actual .and. temp_actual <= minmax(2)) then
call temps_to_read % push_back(nint(temp_actual))
end if
end do
end if
select case (method)
case (TEMPERATURE_NEAREST)
! Find nearest temperatures
@ -671,4 +686,42 @@ module nuclide_header
end function nuclide_nu
subroutine nuclide_init_grid(this, E_min, E_max, M)
class(Nuclide), intent(inout) :: this
real(8), intent(in) :: E_min ! Minimum energy in MeV
real(8), intent(in) :: E_max ! Maximum energy in MeV
integer, intent(in) :: M ! Number of equally log-spaced bins
integer :: i, j, k ! Loop indices
integer :: t ! temperature index
real(8) :: spacing
real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid
! Determine equal-logarithmic energy spacing
spacing = log(E_max/E_min)/M
! Create equally log-spaced energy grid
allocate(umesh(0:M))
umesh(:) = [(i*spacing, i=0, M)]
do t = 1, size(this % grid)
! Allocate logarithmic mapping for nuclide
allocate(this % grid(t) % grid_index(0:M))
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M
do while (log(this % grid(t) % energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(this % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == size(this % grid(t) % energy)) exit
j = j + 1
end do
this % grid(t) % grid_index(k) = j
end do
end do
end subroutine nuclide_init_grid
end module nuclide_header

View file

@ -1,8 +1,10 @@
module output
use, intrinsic :: ISO_C_BINDING
use, intrinsic :: ISO_FORTRAN_ENV
use constants
use eigenvalue, only: openmc_get_keff
use endf, only: reaction_name
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
@ -366,17 +368,23 @@ contains
subroutine print_generation()
integer :: i ! overall generation
integer :: n ! number of active generations
! Determine overall generation and number of active generations
i = overall_generation()
n = i - n_inactive*gen_per_batch
! write out information about batch and generation
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(current_gen))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_generation(overall_gen)
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(overall_gen)
entropy(i)
if (overall_gen - n_inactive*gen_per_batch > 1) then
if (n > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
end if
@ -393,18 +401,25 @@ contains
subroutine print_batch_keff()
integer :: i ! overall generation
integer :: n ! number of active generations
! Determine overall generation and number of active generations
i = overall_generation()
n = i - n_inactive*gen_per_batch
! write out information batch and option independent output
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_generation(overall_gen)
k_generation(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(current_batch*gen_per_batch)
entropy(i)
! write out accumulated k-effective if after first active batch
if (overall_gen - n_inactive*gen_per_batch > 1) then
if (n > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
else
@ -590,50 +605,58 @@ contains
subroutine print_results()
integer :: n ! number of realizations
real(8) :: alpha ! significance level for CI
real(8) :: t_value ! t-value for confidence intervals
real(8) :: t_n1 ! t-value with N-1 degrees of freedom
real(8) :: t_n3 ! t-value with N-3 degrees of freedom
real(8) :: x(2) ! mean and standard deviation
real(C_DOUBLE) :: k_combined(2)
integer(C_INT) :: err
! display header block for results
call header("Results", 4)
n = n_realizations
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1)
! Adjust sum_sq
global_tallies(RESULT_SUM_SQ,:) = t_value * global_tallies(RESULT_SUM_SQ,:)
! Adjust combined estimator
if (n_realizations > 3) then
t_value = t_percentile(ONE - alpha/TWO, n_realizations - 3)
k_combined(2) = t_value * k_combined(2)
end if
t_n1 = t_percentile(ONE - alpha/TWO, n - 1)
t_n3 = t_percentile(ONE - alpha/TWO, n - 3)
else
t_n1 = ONE
t_n3 = ONE
end if
! write global tallies
if (n_realizations > 1) then
if (run_mode == MODE_EIGENVALUE) then
write(ou,102) "k-effective (Collision)", global_tallies(RESULT_SUM, &
K_COLLISION), global_tallies(RESULT_SUM_SQ, K_COLLISION)
write(ou,102) "k-effective (Track-length)", global_tallies(RESULT_SUM, &
K_TRACKLENGTH), global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH)
write(ou,102) "k-effective (Absorption)", global_tallies(RESULT_SUM, &
K_ABSORPTION), global_tallies(RESULT_SUM_SQ, K_ABSORPTION)
if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined
end if
write(ou,102) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE), &
global_tallies(RESULT_SUM_SQ, LEAKAGE)
if (n > 1) then
associate (r => global_tallies(RESULT_SUM:RESULT_SUM_SQ, :))
if (run_mode == MODE_EIGENVALUE) then
x(:) = mean_stdev(r(:, K_COLLISION), n)
write(ou,102) "k-effective (Collision)", x(1), t_n1 * x(2)
x(:) = mean_stdev(r(:, K_TRACKLENGTH), n)
write(ou,102) "k-effective (Track-length)", x(1), t_n1 * x(2)
x(:) = mean_stdev(r(:, K_ABSORPTION), n)
write(ou,102) "k-effective (Absorption)", x(1), t_n1 * x(2)
if (n > 3) then
err = openmc_get_keff(k_combined)
write(ou,102) "Combined k-effective", k_combined(1), &
t_n3 * k_combined(2)
end if
end if
x(:) = mean_stdev(global_tallies(:, LEAKAGE), n)
write(ou,102) "Leakage Fraction", x(1), t_n1 * x(2)
end associate
else
if (master) call warning("Could not compute uncertainties -- only one &
&active batch simulated!")
if (run_mode == MODE_EIGENVALUE) then
write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION)
write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH)
write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION)
write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION) / n
write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION) / n
end if
write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE)
write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE) / n
end if
write(ou,*)
@ -698,8 +721,10 @@ contains
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer :: unit_tally ! tallies.out file unit
integer :: nr ! number of realizations
real(8) :: t_value ! t-values for confidence intervals
real(8) :: alpha ! significance level for CI
real(8) :: x(2) ! mean and standard deviation
character(MAX_FILE_LEN) :: filename ! name of output file
character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function
character(36) :: score_name ! names of scoring function
@ -745,20 +770,20 @@ contains
if (confidence_intervals) then
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1)
else
t_value = ONE
end if
TALLY_LOOP: do i = 1, n_tallies
t => tallies(i)
nr = t % n_realizations
if (confidence_intervals) then
! Calculate t-value for confidence intervals
if (confidence_intervals) then
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, t % n_realizations - 1)
end if
! Multiply uncertainty by t-value
t % results(RESULT_SUM_SQ,:,:) = t_value * t % results(RESULT_SUM_SQ,:,:)
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, nr - 1)
else
t_value = ONE
end if
! Write header block
@ -890,24 +915,27 @@ contains
do l = 1, t % n_user_score_bins
k = k + 1
score_index = score_index + 1
associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :))
select case(t % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
score_name = 'P' // trim(to_str(t % moment_order(k))) // " " // &
score_names(abs(t % score_bins(k)))
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index)))
repeat(" ", indent), score_name, to_str(x(1)), &
trim(to_str(t_value * x(2)))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
do n_order = 0, t % moment_order(k)
score_index = score_index + 1
score_name = 'P' // trim(to_str(n_order)) // " " //&
score_names(abs(t % score_bins(k)))
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index)))
to_str(x(1)), trim(to_str(t_value * x(2)))
end do
k = k + t % moment_order(k)
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
@ -919,11 +947,10 @@ contains
score_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order)) // " " &
// score_names(abs(t % score_bins(k)))
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,&
filter_index)))
to_str(x(1)), trim(to_str(t_value * x(2)))
end do
end do
k = k + (t % moment_order(k) + 1)**2 - 1
@ -933,11 +960,12 @@ contains
else
score_name = score_names(abs(t % score_bins(k)))
end if
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index)))
to_str(x(1)), trim(to_str(t_value * x(2)))
end select
end associate
end do
indent = indent - 2
@ -975,11 +1003,15 @@ contains
integer :: stride_surf ! stride for surface filter
integer :: n ! number of incoming energy bins
integer :: filter_index ! index in results array for filters
integer :: nr ! number of realizations
real(8) :: x(2) ! mean and standard deviation
logical :: print_ebin ! should incoming energy bin be displayed?
logical :: energy_filters ! energy filters present
character(MAX_LINE_LEN) :: string
type(RegularMesh), pointer :: m
nr = t % n_realizations
! Get pointer to mesh
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
select type(filt => filters(i_filter_mesh) % obj)
@ -1050,104 +1082,101 @@ contains
% bins % data(1) - 1) * t % stride(j)
end do
! Left Surface
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Left", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(OUT_LEFT - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(OUT_LEFT - 1) * stride_surf)))
associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :))
! Left Surface
x(:) = mean_stdev(r(:, 1, filter_index + (OUT_LEFT - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Left", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(IN_LEFT - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(IN_LEFT - 1) * stride_surf)))
"Outgoing Current on Left", to_str(x(1)), trim(to_str(x(2)))
x(:) = mean_stdev(r(:, 1, filter_index + (IN_LEFT - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Left", to_str(x(1)), trim(to_str(x(2)))
! Right Surface
x(:) = mean_stdev(r(:, 1, filter_index + (OUT_RIGHT - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Right", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(OUT_RIGHT - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(OUT_RIGHT - 1) * stride_surf)))
"Outgoing Current on Right", to_str(x(1)), trim(to_str(x(2)))
x(:) = mean_stdev(r(:, 1, filter_index + (IN_RIGHT - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Right", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(IN_RIGHT - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(IN_RIGHT - 1) * stride_surf)))
"Incoming Current on Right", to_str(x(1)), trim(to_str(x(2)))
if (n_dim >= 2) then
! Back Surface
x(:) = mean_stdev(r(:, 1, filter_index + (OUT_BACK - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Back", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(OUT_BACK - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(OUT_BACK - 1) * stride_surf)))
"Outgoing Current on Back", to_str(x(1)), trim(to_str(x(2)))
x(:) = mean_stdev(r(:, 1, filter_index + (IN_BACK - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Back", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(IN_BACK - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(IN_BACK - 1) * stride_surf)))
"Incoming Current on Back", to_str(x(1)), trim(to_str(x(2)))
! Front Surface
x(:) = mean_stdev(r(:, 1, filter_index + (OUT_FRONT - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Front", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(OUT_FRONT - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(OUT_FRONT - 1) * stride_surf)))
"Outgoing Current on Front", to_str(x(1)), trim(to_str(x(2)))
x(:) = mean_stdev(r(:, 1, filter_index + (IN_FRONT - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Front", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(IN_FRONT - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(IN_FRONT - 1) * stride_surf)))
"Incoming Current on Front", to_str(x(1)), trim(to_str(x(2)))
end if
if (n_dim == 3) then
! Bottom Surface
x(:) = mean_stdev(r(:, 1, filter_index + (OUT_BOTTOM - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Bottom", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(OUT_BOTTOM - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(OUT_BOTTOM - 1) * stride_surf)))
"Outgoing Current on Bottom", to_str(x(1)), trim(to_str(x(2)))
x(:) = mean_stdev(r(:, 1, filter_index + (IN_BOTTOM - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Bottom", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(IN_BOTTOM - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(IN_BOTTOM - 1) * stride_surf)))
"Incoming Current on Bottom", to_str(x(1)), trim(to_str(x(2)))
! Top Surface
x(:) = mean_stdev(r(:, 1, filter_index + (OUT_TOP - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Top", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(OUT_TOP - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(OUT_TOP - 1) * stride_surf)))
"Outgoing Current on Top", to_str(x(1)), trim(to_str(x(2)))
x(:) = mean_stdev(r(:, 1, filter_index + (IN_TOP - 1) * &
stride_surf), nr)
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Top", &
to_str(t % results(RESULT_SUM,1,filter_index + &
(IN_TOP - 1) * stride_surf)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index + &
(IN_TOP - 1) * stride_surf)))
"Incoming Current on Top", to_str(x(1)), trim(to_str(x(2)))
end if
end associate
end do
end do
end subroutine write_surface_current
!===============================================================================
! MEAN_STDEV computes the sample mean and standard deviation of the mean of a
! single tally score
!===============================================================================
pure function mean_stdev(result_, n) result(x)
real(8), intent(in) :: result_(2) ! sum and sum-of-squares
integer, intent(in) :: n ! number of realizations
real(8) :: x(2) ! mean and standard deviation
x(1) = result_(1) / n
if (n > 1) then
x(2) = sqrt((result_(2) / n - x(1)*x(1))/(n - 1))
else
x(2) = ZERO
end if
end function mean_stdev
end module output

View file

@ -46,6 +46,7 @@ module particle_header
! Particle coordinates
integer :: n_coord ! number of current coordinates
integer :: cell_instance ! offset for distributed properties
type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels
! Particle coordinates before crossing a surface

View file

@ -32,6 +32,8 @@ contains
! Set verbosity high
verbosity = 10
allocate(micro_xs(n_nuclides_total))
! Initialize the particle to be tracked
call p % initialize()
@ -44,8 +46,7 @@ contains
! Compute random number seed
select case (previous_run_mode)
case (MODE_EIGENVALUE)
particle_seed = ((current_batch - 1)*gen_per_batch + &
current_gen - 1)*n_particles + p % id
particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id
case (MODE_FIXEDSOURCE)
particle_seed = p % id
end select
@ -58,6 +59,8 @@ contains
! Write output if particle made it
call print_particle(p)
deallocate(micro_xs)
end subroutine run_particle_restart
!===============================================================================

View file

@ -324,6 +324,7 @@ contains
real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering
real(8) :: phi ! azimuthal angle for iso-in-lab scattering
real(8) :: kT ! temperature in eV
logical :: sampled ! whether or not a reaction type has been sampled
type(Nuclide), pointer :: nuc
! copy incoming direction
@ -337,36 +338,43 @@ contains
! For tallying purposes, this routine might be called directly. In that
! case, we need to sample a reaction via the cutoff variable
prob = ZERO
cutoff = prn() * (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption)
sampled = .false.
prob = prob + micro_xs(i_nuclide) % elastic
prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal
if (prob > cutoff) then
! =======================================================================
! ELASTIC SCATTERING
if (micro_xs(i_nuclide) % index_sab /= NONE) then
! S(a,b) scattering
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, &
p % E, p % coord(1) % uvw, p % mu)
! NON-S(A,B) ELASTIC SCATTERING
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, &
p % E, p % coord(1) % uvw, p % mu, p % wgt)
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
p % event_MT = ELASTIC
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, &
p % coord(1) % uvw, p % mu, p % wgt)
else
p % event_MT = ELASTIC
sampled = .true.
end if
prob = micro_xs(i_nuclide) % elastic
if (prob > cutoff .and. .not. sampled) then
! =======================================================================
! S(A,B) SCATTERING
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, &
p % coord(1) % uvw, p % mu)
p % event_MT = ELASTIC
sampled = .true.
end if
if (.not. sampled) then
! =======================================================================
! INELASTIC SCATTERING
@ -388,7 +396,7 @@ contains
if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF &
.or. rx % MT == N_2NF .or. rx % MT == N_3NF) cycle
! some materials have gas production cross sections with MT > 200 that
! Some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle
@ -406,24 +414,24 @@ contains
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, nuc%reactions(i), p)
p % event_MT = nuc%reactions(i)%MT
p % event_MT = nuc % reactions(i) % MT
end if
! Set event component
p % event = EVENT_SCATTER
! sample new outgoing angle for isotropic in lab scattering
! Sample new outgoing angle for isotropic in lab scattering
if (materials(p % material) % p0(i_nuc_mat)) then
! sample isotropic-in-lab outgoing direction
! Sample isotropic-in-lab outgoing direction
uvw_new(1) = TWO * prn() - ONE
phi = TWO * PI * prn()
uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1))
uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1))
p % mu = dot_product(uvw_old, uvw_new)
! change direction of particle
! Change direction of particle
p % coord(1) % uvw = uvw_new
end if
@ -559,8 +567,8 @@ contains
associate (sab => sab_tables(i_sab) % data(i_temp))
! Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs(i_nuclide) % elastic_sab / &
micro_xs(i_nuclide) % elastic) then
if (prn() < micro_xs(i_nuclide) % thermal_elastic / &
micro_xs(i_nuclide) % thermal) then
! elastic scattering
! Get index and interpolation factor for elastic grid
@ -818,7 +826,6 @@ contains
real(8) :: xs_low ! 0K xs at lowest practical relative energy
real(8) :: xs_up ! 0K xs at highest practical relative energy
real(8) :: m ! slope for interpolation
real(8) :: xi ! pseudorandom number on [0,1)
real(8) :: R ! rejection criterion for DBRC / target speed
real(8) :: cdf_low ! xs cdf at lowest practical relative energy
real(8) :: cdf_up ! xs cdf at highest practical relative energy

View file

@ -1,5 +1,9 @@
module plot
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use error, only: fatal_error
use geometry, only: find_cell, check_cell_overlap
@ -14,12 +18,10 @@ module plot
use progress_header, only: ProgressBar
use string, only: to_str
use hdf5
implicit none
private
public :: run_plot
public :: openmc_plot_geometry
integer, parameter :: RED = 1
integer, parameter :: GREEN = 2
@ -31,7 +33,7 @@ contains
! RUN_PLOT controls the logic for making one or many plots
!===============================================================================
subroutine run_plot()
subroutine openmc_plot_geometry() bind(C)
integer :: i ! loop index for plots
@ -51,7 +53,7 @@ contains
end associate
end do
end subroutine run_plot
end subroutine openmc_plot_geometry
!===============================================================================
! POSITION_RGB computes the red/green/blue values for a given plot with the

View file

@ -28,7 +28,8 @@ element materials {
}* &
element sab {
(element name { xsd:string } | attribute name { xsd:string })
(element name { xsd:string } | attribute name { xsd:string }) &
(element fraction { xsd:double } | attribute fraction { xsd:double })?
}*
}+ &

View file

@ -119,14 +119,26 @@
</zeroOrMore>
<zeroOrMore>
<element name="sab">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="fraction">
<data type="double"/>
</element>
<attribute name="fraction">
<data type="double"/>
</attribute>
</choice>
</optional>
</interleave>
</element>
</zeroOrMore>
</interleave>

View file

@ -116,6 +116,8 @@ element settings {
element temperature_multipole { xsd:boolean }? &
element temperature_range { list { xsd:double, xsd:double } }? &
element temperature_tolerance { xsd:double }? &
element threads { xsd:positiveInteger }? &

View file

@ -517,6 +517,14 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="temperature_range">
<list>
<data type="double"/>
<data type="double"/>
</list>
</element>
</optional>
<optional>
<element name="temperature_tolerance">
<data type="double"/>

View file

@ -80,12 +80,14 @@ module sab_header
contains
subroutine salphabeta_from_hdf5(this, group_id, temperature, method, tolerance)
subroutine salphabeta_from_hdf5(this, group_id, temperature, method, &
tolerance, minmax)
class(SAlphaBeta), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
type(VectorReal), intent(in) :: temperature ! list of temperatures
integer, intent(in) :: method
real(8), intent(in) :: tolerance
real(8), intent(in) :: minmax(2)
integer :: i, j
integer :: t
@ -143,6 +145,18 @@ contains
end do
call sort(temps_available)
! Determine actual temperatures to read -- start by checking whether a
! temperature range was given, in which case all temperatures in the range
! are loaded irrespective of what temperatures actually appear in the model
if (minmax(2) > ZERO) then
do i = 1, size(temps_available)
temp_actual = temps_available(i)
if (minmax(1) <= temp_actual .and. temp_actual <= minmax(2)) then
call temps_to_read % push_back(nint(temp_actual))
end if
end do
end if
select case (method)
case (TEMPERATURE_NEAREST)
! Determine actual temperatures to read

View file

@ -1,10 +1,12 @@
module simulation
use, intrinsic :: ISO_C_BINDING
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
calculate_combined_keff, calculate_generation_keff, &
shannon_entropy, synchronize_bank, keff_generation
calculate_generation_keff, shannon_entropy, &
synchronize_bank, keff_generation, k_sum
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
@ -18,40 +20,28 @@ module simulation
use source, only: initialize_source, sample_external_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
tally_statistics
use tally, only: accumulate_tallies, setup_active_usertallies
use trigger, only: check_triggers
use tracking, only: transport
use volume_calc, only: run_volume_calculations
implicit none
private
public :: run_simulation
public :: openmc_run
contains
!===============================================================================
! RUN_SIMULATION encompasses all the main logic where iterations are performed
! OPENMC_RUN encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a fixed source or k-eigenvalue
! calculation.
!===============================================================================
subroutine run_simulation()
subroutine openmc_run() bind(C)
type(Particle) :: p
integer(8) :: i_work
if (.not. restart_run) call initialize_source()
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", 3)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", 3)
if (verbosity >= 7) call print_columns()
end if
end if
call initialize_simulation()
! Turn on inactive timer
call time_inactive % start()
@ -115,7 +105,7 @@ contains
! Clear particle
call p % clear()
end subroutine run_simulation
end subroutine openmc_run
!===============================================================================
! INITIALIZE_HISTORY
@ -137,7 +127,7 @@ contains
p % id = work_index(rank) + index_source
! set random number seed
particle_seed = (overall_gen - 1)*n_particles + p % id
particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! set particle trace
@ -202,9 +192,6 @@ contains
subroutine initialize_generation()
! set overall generation number
overall_gen = gen_per_batch*(current_batch - 1) + current_gen
if (run_mode == MODE_EIGENVALUE) then
! Reset number of fission bank sites
n_bank = 0
@ -269,13 +256,20 @@ contains
call calculate_average_keff()
! Write generation output
if (master .and. current_gen /= gen_per_batch .and. verbosity >= 7) &
call print_generation()
if (master .and. verbosity >= 7) then
if (current_gen /= gen_per_batch) then
call print_generation()
else
call print_batch_keff()
end if
end if
elseif (run_mode == MODE_FIXEDSOURCE) then
! For fixed-source mode, we need to sample the external source
if (path_source == '') then
do i = 1, work
call set_particle_seed(overall_gen*n_particles + work_index(rank) + i)
call set_particle_seed((total_gen + overall_generation()) * &
n_particles + work_index(rank) + i)
call sample_external_source(source_bank(i))
end do
end if
@ -291,9 +285,13 @@ contains
subroutine finalize_batch()
! Collect tallies
#ifdef MPI
integer :: mpi_err ! MPI error code
#endif
! Reduce tallies onto master process and accumulate
call time_tallies % start()
call synchronize_tallies()
call accumulate_tallies()
call time_tallies % stop()
! Reset global tally results
@ -305,12 +303,6 @@ contains
if (run_mode == MODE_EIGENVALUE) then
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Display output
if (master .and. verbosity >= 7) call print_batch_keff()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
end if
! Check_triggers
@ -335,13 +327,6 @@ contains
call write_source_point()
end if
if (master .and. current_batch == n_max_batches .and. &
run_mode == MODE_EIGENVALUE) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
end if
end subroutine finalize_batch
!===============================================================================
@ -358,7 +343,6 @@ contains
if (run_mode == MODE_EIGENVALUE) then
do current_gen = 1, gen_per_batch
overall_gen = overall_gen + 1
call calculate_average_keff()
! print out batch keff
@ -379,24 +363,88 @@ contains
end subroutine replay_batch_history
!===============================================================================
! INITIALIZE_SIMULATION
!===============================================================================
subroutine initialize_simulation()
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
if (.not. restart_run) call initialize_source()
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", 3)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", 3)
if (verbosity >= 7) call print_columns()
end if
end if
end subroutine initialize_simulation
!===============================================================================
! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
! execution time and results
!===============================================================================
subroutine finalize_simulation
subroutine finalize_simulation()
#ifdef MPI
integer :: i ! loop index for tallies
integer :: n ! size of arrays
integer :: mpi_err ! MPI error code
integer(8) :: temp
real(8) :: tempr(3) ! temporary array for communication
#endif
!$omp parallel
deallocate(micro_xs)
!$omp end parallel
! Increment total number of generations
total_gen = total_gen + n_batches*gen_per_batch
! Start finalization timer
call time_finalize%start()
call time_finalize % start()
! Calculate statistics for tallies and write to tallies.out
if (master) then
if (n_realizations > 1) call tally_statistics()
#ifdef MPI
! Broadcast tally results so that each process has access to results
if (allocated(tallies)) then
do i = 1, size(tallies)
n = size(tallies(i) % results)
call MPI_BCAST(tallies(i) % results, n, MPI_DOUBLE, 0, &
mpi_intracomm, mpi_err)
end do
end if
if (output_tallies) then
if (master) call write_tallies()
! Also broadcast global tally results
n = size(global_tallies)
call MPI_BCAST(global_tallies, n, MPI_DOUBLE, 0, mpi_intracomm, mpi_err)
! These guys are needed so that non-master processes can calculate the
! combined estimate of k-effective
tempr(1) = k_col_abs
tempr(2) = k_col_tra
tempr(3) = k_abs_tra
call MPI_BCAST(tempr, 3, MPI_REAL8, 0, mpi_intracomm, mpi_err)
k_col_abs = tempr(1)
k_col_tra = tempr(2)
k_abs_tra = tempr(3)
if (check_overlaps) then
call MPI_REDUCE(overlap_check_cnt, temp, n_cells, MPI_INTEGER8, &
MPI_SUM, 0, mpi_intracomm, mpi_err)
overlap_check_cnt = temp
end if
if (check_overlaps) call reduce_overlap_count()
#endif
! Write tally results to tallies.out
if (output_tallies .and. master) call write_tallies()
! Stop timers and show timing statistics
call time_finalize%stop()
@ -409,22 +457,4 @@ contains
end subroutine finalize_simulation
!===============================================================================
! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master
!===============================================================================
subroutine reduce_overlap_count()
#ifdef MPI
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, &
MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err)
else
call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, &
MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err)
end if
#endif
end subroutine reduce_overlap_count
end module simulation

View file

@ -73,7 +73,7 @@ contains
src => source_bank(i)
! initialize random number seed
id = work_index(rank) + i
id = total_gen*n_particles + work_index(rank) + i
call set_particle_seed(id)
! sample external source distribution

View file

@ -11,11 +11,12 @@ module state_point
! intervals, using the <state_point ... /> tag.
!===============================================================================
use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use eigenvalue, only: openmc_get_keff
use endf, only: reaction_name
use error, only: fatal_error, warning
use global
@ -46,6 +47,8 @@ contains
integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
mesh_group, filters_group, filter_group, derivs_group, &
deriv_group, runtime_group
integer(C_INT) :: err
real(C_DOUBLE) :: k_combined(2)
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(MAX_FILE_LEN) :: filename
type(TallyObject), pointer :: tally
@ -117,6 +120,7 @@ contains
call write_dataset(file_id, "k_col_abs", k_col_abs)
call write_dataset(file_id, "k_col_tra", k_col_tra)
call write_dataset(file_id, "k_abs_tra", k_abs_tra)
err = openmc_get_keff(k_combined)
call write_dataset(file_id, "k_combined", k_combined)
! Write out CMFD info
@ -518,7 +522,8 @@ contains
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results
real(8), target :: global_temp(3,N_GLOBAL_TALLIES)
#ifdef MPI
real(8) :: dummy ! temporary receive buffer for non-root reduces
integer :: mpi_err ! MPI error code
real(8) :: dummy ! temporary receive buffer for non-root reduces
#endif
type(TallyObject) :: dummy_tally
@ -535,18 +540,15 @@ contains
tallies_group = open_group(file_id, "tallies")
end if
! Copy global tallies into temporary array for reducing
n_bins = 3 * N_GLOBAL_TALLIES
global_temp(:,:) = global_tallies(:,:)
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
#ifdef MPI
call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
! Reduce global tallies
n_bins = size(global_tallies)
call MPI_REDUCE(global_tallies, global_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
#endif
if (master) then
! Transfer values to value on master
if (current_batch == n_max_batches .or. satisfy_triggers) then
global_tallies(:,:) = global_temp(:,:)
@ -554,12 +556,6 @@ contains
! Write out global tallies sum and sum_sq
call write_dataset(file_id, "global_tallies", global_temp)
else
! Receive buffer not significant at other processors
#ifdef MPI
call MPI_REDUCE(global_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
#endif
end if
if (tallies_on) then
@ -647,7 +643,6 @@ contains
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
real(8) :: real_array(3)
logical :: source_present
character(MAX_WORD_LEN) :: word
@ -727,7 +722,6 @@ contains
call read_dataset(k_col_abs, file_id, "k_col_abs")
call read_dataset(k_col_tra, file_id, "k_col_tra")
call read_dataset(k_abs_tra, file_id, "k_abs_tra")
call read_dataset(real_array(1:2), file_id, "k_combined")
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))
@ -848,6 +842,7 @@ contains
#else
integer :: i
#ifdef MPI
integer :: mpi_err ! MPI error code
type(Bank), allocatable, target :: temp_source(:)
#endif
#endif

View file

@ -22,10 +22,6 @@ module tally
implicit none
integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array
!$omp threadprivate(position)
procedure(score_general_), pointer :: score_general => null()
procedure(score_analog_tally_), pointer :: score_analog_tally => null()
@ -2410,9 +2406,6 @@ contains
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
! Reset tally map positioning
position = 0
end subroutine score_analog_tally_ce
subroutine score_analog_tally_mg(p)
@ -2558,9 +2551,6 @@ contains
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
! Reset tally map positioning
position = 0
end subroutine score_analog_tally_mg
!===============================================================================
@ -2956,9 +2946,6 @@ contains
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
! Reset tally map positioning
position = 0
end subroutine score_tracklength_tally
!===============================================================================
@ -3132,9 +3119,6 @@ contains
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
! Reset tally map positioning
position = 0
end subroutine score_collision_tally
!===============================================================================
@ -4334,11 +4318,11 @@ contains
end subroutine zero_flux_derivs
!===============================================================================
! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history
! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history
! within the batch to a new random variable
!===============================================================================
subroutine synchronize_tallies()
subroutine accumulate_tallies()
integer :: i
real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff
@ -4388,7 +4372,7 @@ contains
end do
end if
end subroutine synchronize_tallies
end subroutine accumulate_tallies
!===============================================================================
! REDUCE_TALLY_RESULTS collects all the results from tallies onto one processor
@ -4398,73 +4382,55 @@ contains
subroutine reduce_tally_results()
integer :: i
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:) ! contiguous array of results
real(8) :: global_temp(N_GLOBAL_TALLIES)
real(8) :: dummy ! temporary receive buffer for non-root reduces
type(TallyObject), pointer :: t
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: n_bins ! total number of bins
integer :: mpi_err ! MPI error code
real(C_DOUBLE), allocatable :: tally_temp(:,:) ! contiguous array of results
real(C_DOUBLE), allocatable :: tally_temp2(:,:) ! reduced contiguous results
real(C_DOUBLE) :: temp(N_GLOBAL_TALLIES), temp2(N_GLOBAL_TALLIES)
do i = 1, active_tallies % size()
t => tallies(active_tallies % data(i))
associate (t => tallies(active_tallies % data(i)))
m = t % total_score_bins
n = t % total_filter_bins
n_bins = m*n
m = size(t % results, 2)
n = size(t % results, 3)
n_bins = m*n
allocate(tally_temp(m,n))
allocate(tally_temp(m,n), tally_temp2(m,n))
tally_temp = t % results(RESULT_VALUE,:,:)
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into
! a receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
! Reduce contiguous set of tally results
tally_temp = t % results(RESULT_VALUE,:,:)
call MPI_REDUCE(tally_temp, tally_temp2, n_bins, MPI_DOUBLE, &
MPI_SUM, 0, mpi_intracomm, mpi_err)
! Transfer values to value on master
t % results(RESULT_VALUE,:,:) = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, &
MPI_SUM, 0, mpi_intracomm, mpi_err)
if (master) then
! Transfer values to value on master
t % results(RESULT_VALUE,:,:) = tally_temp2
else
! Reset value on other processors
t % results(RESULT_VALUE,:,:) = ZERO
end if
! Reset value on other processors
t % results(RESULT_VALUE,:,:) = ZERO
end if
deallocate(tally_temp)
deallocate(tally_temp, tally_temp2)
end associate
end do
! Copy global tallies into array to be reduced
global_temp = global_tallies(RESULT_VALUE, :)
! Reduce global tallies onto master
temp = global_tallies(RESULT_VALUE, :)
call MPI_REDUCE(temp, temp2, N_GLOBAL_TALLIES, MPI_DOUBLE, MPI_SUM, &
0, mpi_intracomm, mpi_err)
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
! Transfer values back to global_tallies on master
global_tallies(RESULT_VALUE, :) = global_temp
global_tallies(RESULT_VALUE, :) = temp2
else
! Receive buffer not significant at other processors
call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
! Reset value on other processors
global_tallies(RESULT_VALUE, :) = ZERO
end if
! We also need to determine the total starting weight of particles from the
! last realization
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
end if
temp(1) = total_weight
call MPI_REDUCE(temp, total_weight, 1, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
end subroutine reduce_tally_results
#endif
@ -4502,45 +4468,6 @@ contains
end subroutine accumulate_tally
!===============================================================================
! TALLY_STATISTICS computes the mean and standard deviation of the mean of each
! tally and stores them in the RESULT_SUM and RESULT_SUM_SQ positions
!===============================================================================
subroutine tally_statistics()
integer :: i ! index in tallies array
integer :: j, k ! score/filter indices
integer :: n ! number of realizations
! Calculate sample mean and standard deviation of the mean -- note that we
! have used Bessel's correction so that the estimator of the variance of the
! sample mean is unbiased.
do i = 1, n_tallies
n = tallies(i) % n_realizations
associate (r => tallies(i) % results)
do k = 1, size(r, 3)
do j = 1, size(r, 2)
r(RESULT_SUM, j, k) = r(RESULT_SUM, j, k) / n
r(RESULT_SUM_SQ, j, k) = sqrt((r(RESULT_SUM_SQ, j, k)/n - &
r(RESULT_SUM, j, k) * r(RESULT_SUM, j, k))/(n - 1))
end do
end do
end associate
end do
! Calculate statistics for global tallies
n = n_realizations
associate (r => global_tallies)
do j = 1, size(r, 2)
r(RESULT_SUM, j) = r(RESULT_SUM, j) / n
r(RESULT_SUM_SQ, j) = sqrt((r(RESULT_SUM_SQ, j)/n - &
r(RESULT_SUM, j) * r(RESULT_SUM, j))/(n - 1))
end do
end associate
end subroutine tally_statistics
!===============================================================================
! SETUP_ACTIVE_USERTALLIES
!===============================================================================

View file

@ -1,10 +1,13 @@
module trigger
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
use message_passing
#endif
use constants
use eigenvalue, only: openmc_get_keff
use global
use string, only: to_str
use output, only: warning, write_message
@ -101,10 +104,12 @@ contains
integer :: score_index ! scoring bin index
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer(C_INT) :: err
real(8) :: uncertainty ! trigger uncertainty
real(8) :: std_dev = ZERO ! trigger standard deviation
real(8) :: rel_err = ZERO ! trigger relative error
real(8) :: ratio ! ratio of the uncertainty/trigger threshold
real(C_DOUBLE) :: k_combined(2)
type(TallyObject), pointer :: t ! tally pointer
type(TriggerObject), pointer :: trigger ! tally trigger
@ -119,6 +124,7 @@ contains
! Check eigenvalue trigger
if (run_mode == MODE_EIGENVALUE) then
if (keff_trigger % trigger_type /= 0) then
err = openmc_get_keff(k_combined)
select case (keff_trigger % trigger_type)
case(VARIANCE)
uncertainty = k_combined(2) ** 2

View file

@ -1,5 +1,7 @@
module volume_calc
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
#ifdef _OPENMP
use omp_lib
@ -21,16 +23,16 @@ module volume_calc
implicit none
private
public :: run_volume_calculations
public :: openmc_calculate_volumes
contains
!===============================================================================
! RUN_VOLUME_CALCULATIONS runs each of the stochastic volume calculations that
! OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations that
! the user has specified and writes results to HDF5 files
!===============================================================================
subroutine run_volume_calculations()
subroutine openmc_calculate_volumes() bind(C)
integer :: i, j
integer :: n
real(8), allocatable :: volume(:,:) ! volume mean/stdev in each domain
@ -93,7 +95,7 @@ contains
call write_message("Elapsed time: " // trim(to_str(time_volume % &
get_value())) // " s", 6)
end if
end subroutine run_volume_calculations
end subroutine openmc_calculate_volumes
!===============================================================================
! GET_VOLUME stochastically determines the volume of a set of domains along with
@ -130,6 +132,7 @@ contains
integer :: min_samples ! minimum number of samples per process
integer :: remainder ! leftover samples from uneven divide
#ifdef MPI
integer :: mpi_err ! MPI error code
integer :: m ! index over materials
integer :: n ! number of materials
integer, allocatable :: data(:) ! array used to send number of hits

View file

@ -188,6 +188,10 @@ class Test(object):
build_str += "-Dopenmp=OFF "
if self.coverage:
build_str += "-Dcoverage=ON "
if self.phdf5:
build_str += "-DHDF5_PREFER_PARALLEL=ON "
else:
build_str += "-DHDF5_PREFER_PARALLEL=OFF "
self.build_opts = build_str
return self.build_opts
@ -214,6 +218,9 @@ class Test(object):
# Runs cmake when in non-script mode
def run_cmake(self):
build_opts = self.build_opts.split()
self.cmake += build_opts
os.environ['FC'] = self.fc
os.environ['CC'] = self.cc
os.environ['CXX'] = self.cxx
@ -221,10 +228,10 @@ class Test(object):
os.environ['MPI_DIR'] = MPI_DIR
if self.phdf5:
os.environ['HDF5_ROOT'] = PHDF5_DIR
self.cmake.append('-DHDF5_PREFER_PARALLEL=ON')
else:
os.environ['HDF5_ROOT'] = HDF5_DIR
build_opts = self.build_opts.split()
self.cmake += build_opts
self.cmake.append('-DHDF5_PREFER_PARALLEL=OFF')
rc = call(self.cmake)
if rc != 0:
self.success = False

View file

@ -1,19 +1,19 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell fill="202" id="1" region="10000 -10001 10002 -10003 10004 -10005" universe="0" />
<cell id="27" material="1" region="-1" universe="3" />
<cell id="28" material="2" region="1 -2" universe="3" />
<cell id="29" material="4" region="2" universe="3" />
<cell id="30" material="4" region="-3" universe="4" />
<cell id="31" material="2" region="3 -4" universe="4" />
<cell id="32" material="4" region="4" universe="4" />
<cell id="70" material="4" region="35 -36" universe="7" />
<cell fill="101" id="80" region="35 -36" universe="8" />
<lattice id="101" name="Fuel assembly (upper half)">
<pitch>1.26 1.26</pitch>
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<universes>
<cell fill="202" id="1" region="9 -10 11 -12 13 -14" universe="0" />
<cell id="27" material="1" region="-1" universe="3" />
<cell id="28" material="2" region="1 -2" universe="3" />
<cell id="29" material="4" region="2" universe="3" />
<cell id="30" material="4" region="-3" universe="4" />
<cell id="31" material="2" region="3 -4" universe="4" />
<cell id="32" material="4" region="4" universe="4" />
<cell id="70" material="4" region="35 -36" universe="7" />
<cell fill="101" id="80" region="35 -36" universe="8" />
<lattice id="101" name="Fuel assembly (upper half)">
<pitch>1.26 1.26</pitch>
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
@ -31,197 +31,197 @@
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 </universes>
</lattice>
<lattice id="202" name="3x3 Core Lattice">
<pitch>21.42 21.42</pitch>
<dimension>3 3</dimension>
<lower_left>-32.13 -32.13</lower_left>
<universes>
</lattice>
<lattice id="202" name="3x3 Core Lattice">
<pitch>21.42 21.42</pitch>
<dimension>3 3</dimension>
<lower_left>-32.13 -32.13</lower_left>
<universes>
8 7 7
8 8 8
7 7 7 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.41" id="1" type="z-cylinder" />
<surface boundary="reflective" coeffs="-32.13" id="10000" type="x-plane" />
<surface boundary="reflective" coeffs="32.13" id="10001" type="x-plane" />
<surface boundary="reflective" coeffs="-32.13" id="10002" type="y-plane" />
<surface boundary="reflective" coeffs="32.13" id="10003" type="y-plane" />
<surface boundary="reflective" coeffs="0" id="10004" type="z-plane" />
<surface boundary="reflective" coeffs="32.13" id="10005" type="z-plane" />
<surface coeffs="0.0 0.0 0.475" id="2" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.56" id="3" type="z-cylinder" />
<surface coeffs="0.0" id="35" type="z-plane" />
<surface coeffs="183.0" id="36" type="z-plane" />
<surface coeffs="0.0 0.0 0.62" id="4" type="z-cylinder" />
</lattice>
<surface coeffs="0.0 0.0 0.41" id="1" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.475" id="2" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.56" id="3" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.62" id="4" type="z-cylinder" />
<surface boundary="reflective" coeffs="-32.13" id="9" type="x-plane" />
<surface boundary="reflective" coeffs="32.13" id="10" type="x-plane" />
<surface boundary="reflective" coeffs="-32.13" id="11" type="y-plane" />
<surface boundary="reflective" coeffs="32.13" id="12" type="y-plane" />
<surface boundary="reflective" coeffs="0" id="13" type="z-plane" />
<surface boundary="reflective" coeffs="32.13" id="14" type="z-plane" />
<surface coeffs="0.0" id="35" type="z-plane" />
<surface coeffs="183.0" id="36" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="1" name="UOX fuel">
<density units="g/cm3" value="10.062" />
<nuclide ao="4.9476e-06" name="U234" />
<nuclide ao="0.00048218" name="U235" />
<nuclide ao="0.021504" name="U238" />
<nuclide ao="1.0801e-08" name="Xe135" />
<nuclide ao="0.045737" name="O16" />
</material>
<material id="10" name="Top nozzle region">
<density units="g/cm3" value="1.746" />
<nuclide name="H1" wo="0.035887" />
<nuclide name="O16" wo="0.2847761" />
<nuclide name="B10" wo="0.000115699" />
<nuclide name="B11" wo="0.000527075" />
<nuclide name="Fe54" wo="0.02644016154" />
<nuclide name="Fe56" wo="0.43037146399" />
<nuclide name="Fe57" wo="0.0101152584" />
<nuclide name="Fe58" wo="0.00137211607" />
<nuclide name="Ni58" wo="0.04104621835" />
<nuclide name="Mn55" wo="0.0135739" />
<nuclide name="Cr52" wo="0.107931450781" />
<sab name="c_H_in_H2O" />
</material>
<material id="11" name="Top of fuel assemblies">
<density units="g/cm3" value="3.044" />
<nuclide name="H1" wo="0.0162913" />
<nuclide name="O16" wo="0.1292776" />
<nuclide name="B10" wo="5.25228e-05" />
<nuclide name="B11" wo="0.000239272" />
<nuclide name="Zr90" wo="0.43313403903" />
<nuclide name="Zr91" wo="0.09549277374" />
<nuclide name="Zr92" wo="0.14759527104" />
<nuclide name="Zr94" wo="0.15280552077" />
<nuclide name="Zr96" wo="0.02511169542" />
<sab name="c_H_in_H2O" />
</material>
<material id="12" name="Bottom of fuel assemblies">
<density units="g/cm3" value="1.762" />
<nuclide name="H1" wo="0.0292856" />
<nuclide name="O16" wo="0.2323919" />
<nuclide name="B10" wo="9.44159e-05" />
<nuclide name="B11" wo="0.00043012" />
<nuclide name="Zr90" wo="0.3741373658" />
<nuclide name="Zr91" wo="0.0824858164" />
<nuclide name="Zr92" wo="0.1274914944" />
<nuclide name="Zr94" wo="0.1319920622" />
<nuclide name="Zr96" wo="0.0216912612" />
<sab name="c_H_in_H2O" />
</material>
<material id="2" name="Zircaloy">
<density units="g/cm3" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="3" name="Cold borated water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<nuclide ao="0.000649" name="B10" />
<nuclide ao="0.002689" name="B11" />
<sab name="c_H_in_H2O" />
</material>
<material id="4" name="Hot borated water">
<density units="atom/b-cm" value="0.06614" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<nuclide ao="0.000649" name="B10" />
<nuclide ao="0.002689" name="B11" />
<sab name="c_H_in_H2O" />
</material>
<material id="5" name="Reactor pressure vessel steel">
<density units="g/cm3" value="7.9" />
<nuclide name="Fe54" wo="0.05437098" />
<nuclide name="Fe56" wo="0.88500663" />
<nuclide name="Fe57" wo="0.0208008" />
<nuclide name="Fe58" wo="0.00282159" />
<nuclide name="Ni58" wo="0.0067198" />
<nuclide name="Ni60" wo="0.0026776" />
<nuclide name="Mn55" wo="0.01" />
<nuclide name="Cr52" wo="0.002092475" />
<nuclide name="C0" wo="0.0025" />
<nuclide name="Cu63" wo="0.0013696" />
</material>
<material id="6" name="Lower radial reflector">
<density units="g/cm3" value="4.32" />
<nuclide name="H1" wo="0.0095661" />
<nuclide name="O16" wo="0.0759107" />
<nuclide name="B10" wo="3.08409e-05" />
<nuclide name="B11" wo="0.000140499" />
<nuclide name="Fe54" wo="0.035620772088" />
<nuclide name="Fe56" wo="0.579805982228" />
<nuclide name="Fe57" wo="0.01362750048" />
<nuclide name="Fe58" wo="0.001848545204" />
<nuclide name="Ni58" wo="0.055298376566" />
<nuclide name="Mn55" wo="0.018287" />
<nuclide name="Cr52" wo="0.145407678031" />
<sab name="c_H_in_H2O" />
</material>
<material id="7" name="Upper radial reflector / Top plate region">
<density units="g/cm3" value="4.28" />
<nuclide name="H1" wo="0.0086117" />
<nuclide name="O16" wo="0.0683369" />
<nuclide name="B10" wo="2.77638e-05" />
<nuclide name="B11" wo="0.000126481" />
<nuclide name="Fe54" wo="0.035953677186" />
<nuclide name="Fe56" wo="0.585224740891" />
<nuclide name="Fe57" wo="0.01375486056" />
<nuclide name="Fe58" wo="0.001865821363" />
<nuclide name="Ni58" wo="0.055815129186" />
<nuclide name="Mn55" wo="0.0184579" />
<nuclide name="Cr52" wo="0.146766614995" />
<sab name="c_H_in_H2O" />
</material>
<material id="8" name="Bottom plate region">
<density units="g/cm3" value="7.184" />
<nuclide name="H1" wo="0.0011505" />
<nuclide name="O16" wo="0.0091296" />
<nuclide name="B10" wo="3.70915e-06" />
<nuclide name="B11" wo="1.68974e-05" />
<nuclide name="Fe54" wo="0.03855611055" />
<nuclide name="Fe56" wo="0.627585036425" />
<nuclide name="Fe57" wo="0.014750478" />
<nuclide name="Fe58" wo="0.002000875025" />
<nuclide name="Ni58" wo="0.059855207342" />
<nuclide name="Mn55" wo="0.019794" />
<nuclide name="Cr52" wo="0.157390026871" />
<sab name="c_H_in_H2O" />
</material>
<material id="9" name="Bottom nozzle region">
<density units="g/cm3" value="2.53" />
<nuclide name="H1" wo="0.0245014" />
<nuclide name="O16" wo="0.1944274" />
<nuclide name="B10" wo="7.89917e-05" />
<nuclide name="B11" wo="0.000359854" />
<nuclide name="Fe54" wo="0.030411411144" />
<nuclide name="Fe56" wo="0.495012237964" />
<nuclide name="Fe57" wo="0.01163454624" />
<nuclide name="Fe58" wo="0.001578204652" />
<nuclide name="Ni58" wo="0.047211231662" />
<nuclide name="Mn55" wo="0.0156126" />
<nuclide name="Cr52" wo="0.124142524198" />
<sab name="c_H_in_H2O" />
</material>
<material id="1" name="UOX fuel">
<density units="g/cm3" value="10.062" />
<nuclide ao="4.9476e-06" name="U234" />
<nuclide ao="0.00048218" name="U235" />
<nuclide ao="0.021504" name="U238" />
<nuclide ao="1.0801e-08" name="Xe135" />
<nuclide ao="0.045737" name="O16" />
</material>
<material id="2" name="Zircaloy">
<density units="g/cm3" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="3" name="Cold borated water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<nuclide ao="0.000649" name="B10" />
<nuclide ao="0.002689" name="B11" />
<sab name="c_H_in_H2O" />
</material>
<material id="4" name="Hot borated water">
<density units="atom/b-cm" value="0.06614" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<nuclide ao="0.000649" name="B10" />
<nuclide ao="0.002689" name="B11" />
<sab name="c_H_in_H2O" />
</material>
<material id="5" name="Reactor pressure vessel steel">
<density units="g/cm3" value="7.9" />
<nuclide name="Fe54" wo="0.05437098" />
<nuclide name="Fe56" wo="0.88500663" />
<nuclide name="Fe57" wo="0.0208008" />
<nuclide name="Fe58" wo="0.00282159" />
<nuclide name="Ni58" wo="0.0067198" />
<nuclide name="Ni60" wo="0.0026776" />
<nuclide name="Mn55" wo="0.01" />
<nuclide name="Cr52" wo="0.002092475" />
<nuclide name="C0" wo="0.0025" />
<nuclide name="Cu63" wo="0.0013696" />
</material>
<material id="6" name="Lower radial reflector">
<density units="g/cm3" value="4.32" />
<nuclide name="H1" wo="0.0095661" />
<nuclide name="O16" wo="0.0759107" />
<nuclide name="B10" wo="3.08409e-05" />
<nuclide name="B11" wo="0.000140499" />
<nuclide name="Fe54" wo="0.035620772088" />
<nuclide name="Fe56" wo="0.579805982228" />
<nuclide name="Fe57" wo="0.01362750048" />
<nuclide name="Fe58" wo="0.001848545204" />
<nuclide name="Ni58" wo="0.055298376566" />
<nuclide name="Mn55" wo="0.018287" />
<nuclide name="Cr52" wo="0.145407678031" />
<sab name="c_H_in_H2O" />
</material>
<material id="7" name="Upper radial reflector / Top plate region">
<density units="g/cm3" value="4.28" />
<nuclide name="H1" wo="0.0086117" />
<nuclide name="O16" wo="0.0683369" />
<nuclide name="B10" wo="2.77638e-05" />
<nuclide name="B11" wo="0.000126481" />
<nuclide name="Fe54" wo="0.035953677186" />
<nuclide name="Fe56" wo="0.585224740891" />
<nuclide name="Fe57" wo="0.01375486056" />
<nuclide name="Fe58" wo="0.001865821363" />
<nuclide name="Ni58" wo="0.055815129186" />
<nuclide name="Mn55" wo="0.0184579" />
<nuclide name="Cr52" wo="0.146766614995" />
<sab name="c_H_in_H2O" />
</material>
<material id="8" name="Bottom plate region">
<density units="g/cm3" value="7.184" />
<nuclide name="H1" wo="0.0011505" />
<nuclide name="O16" wo="0.0091296" />
<nuclide name="B10" wo="3.70915e-06" />
<nuclide name="B11" wo="1.68974e-05" />
<nuclide name="Fe54" wo="0.03855611055" />
<nuclide name="Fe56" wo="0.627585036425" />
<nuclide name="Fe57" wo="0.014750478" />
<nuclide name="Fe58" wo="0.002000875025" />
<nuclide name="Ni58" wo="0.059855207342" />
<nuclide name="Mn55" wo="0.019794" />
<nuclide name="Cr52" wo="0.157390026871" />
<sab name="c_H_in_H2O" />
</material>
<material id="9" name="Bottom nozzle region">
<density units="g/cm3" value="2.53" />
<nuclide name="H1" wo="0.0245014" />
<nuclide name="O16" wo="0.1944274" />
<nuclide name="B10" wo="7.89917e-05" />
<nuclide name="B11" wo="0.000359854" />
<nuclide name="Fe54" wo="0.030411411144" />
<nuclide name="Fe56" wo="0.495012237964" />
<nuclide name="Fe57" wo="0.01163454624" />
<nuclide name="Fe58" wo="0.001578204652" />
<nuclide name="Ni58" wo="0.047211231662" />
<nuclide name="Mn55" wo="0.0156126" />
<nuclide name="Cr52" wo="0.124142524198" />
<sab name="c_H_in_H2O" />
</material>
<material id="10" name="Top nozzle region">
<density units="g/cm3" value="1.746" />
<nuclide name="H1" wo="0.035887" />
<nuclide name="O16" wo="0.2847761" />
<nuclide name="B10" wo="0.000115699" />
<nuclide name="B11" wo="0.000527075" />
<nuclide name="Fe54" wo="0.02644016154" />
<nuclide name="Fe56" wo="0.43037146399" />
<nuclide name="Fe57" wo="0.0101152584" />
<nuclide name="Fe58" wo="0.00137211607" />
<nuclide name="Ni58" wo="0.04104621835" />
<nuclide name="Mn55" wo="0.0135739" />
<nuclide name="Cr52" wo="0.107931450781" />
<sab name="c_H_in_H2O" />
</material>
<material id="11" name="Top of fuel assemblies">
<density units="g/cm3" value="3.044" />
<nuclide name="H1" wo="0.0162913" />
<nuclide name="O16" wo="0.1292776" />
<nuclide name="B10" wo="5.25228e-05" />
<nuclide name="B11" wo="0.000239272" />
<nuclide name="Zr90" wo="0.43313403903" />
<nuclide name="Zr91" wo="0.09549277374" />
<nuclide name="Zr92" wo="0.14759527104" />
<nuclide name="Zr94" wo="0.15280552077" />
<nuclide name="Zr96" wo="0.02511169542" />
<sab name="c_H_in_H2O" />
</material>
<material id="12" name="Bottom of fuel assemblies">
<density units="g/cm3" value="1.762" />
<nuclide name="H1" wo="0.0292856" />
<nuclide name="O16" wo="0.2323919" />
<nuclide name="B10" wo="9.44159e-05" />
<nuclide name="B11" wo="0.00043012" />
<nuclide name="Zr90" wo="0.3741373658" />
<nuclide name="Zr91" wo="0.0824858164" />
<nuclide name="Zr92" wo="0.1274914944" />
<nuclide name="Zr94" wo="0.1319920622" />
<nuclide name="Zr96" wo="0.0216912612" />
<sab name="c_H_in_H2O" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="fission">
<parameters>-32 -32 0 32 32 32</parameters>
</space>
</source>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="fission">
<parameters>-32 -32 0 32 32 32</parameters>
</space>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<filter id="10000" type="distribcell">
<bins>27</bins>
</filter>
<tally id="27" name="distribcell tally">
<filters>10000</filters>
<scores>nu-fission</scores>
</tally>
<filter id="1" type="distribcell">
<bins>27</bins>
</filter>
<tally id="27" name="distribcell tally">
<filters>1</filters>
<scores>nu-fission</scores>
</tally>
</tallies>

View file

@ -1 +1 @@
f388c07481a7333bde8044d2e8805954ec41cc3e73d00f1b92b6acff03dad06e62a74e8b7595994e877f69a051ff34d8e498de69ac0efec6e396434007c10b0c
3b76b468b9c0e7df6508189b75748a1b7b4f2b37486396749e1a15e536526336f70b60bb207095c39ecbd46822e8c8705ea81184a3c8546e7da09bb905d74637

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