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

Merge pull request #903 from paulromano/c-api

Browse source 
The OpenMC C API
This commit is contained in:
Sterling Harper 2017-08-03 13:03:20 -04:00 • committed by GitHub
commit 5776eb89d4
51 changed files with 2938 additions and 873 deletions
Download patch file Download diff file Expand all files Collapse all files

108
CMakeLists.txt
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
#===============================================================================
@ -100,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
@ -109,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()
@ -214,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}")
@ -264,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
@ -281,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
@ -346,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})
@ -385,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

6
docs/source/conf.py
View file

@ -25,9 +25,9 @@ except ImportError:
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'scipy', 'scipy.sparse', 'scipy.interpolate', 'scipy.integrate',
'scipy.optimize', 'scipy.special', '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

1
docs/source/index.rst
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

11
docs/source/io_formats/settings.rst
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:
-----------------------------------

39
docs/source/pythonapi/capi.rst Normal file
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

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

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

21
examples/xml/basic/tallies.xml
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>

6
examples/xml/lattice/nested/tallies.xml
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>

6
examples/xml/lattice/simple/tallies.xml
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>

11
examples/xml/pincell/tallies.xml
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>

9
examples/xml/pincell_multigroup/tallies.xml
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>

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()

4
openmc/nuclide.py
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

23
openmc/settings.py
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)

17
setup.py
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',
@ -55,10 +68,6 @@ if have_setuptools:
'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

6
src/cmfd_execute.F90
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)

1
src/cross_section.F90
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

53
src/eigenvalue.F90
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

67
src/energy_grid.F90
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

41
src/finalize.F90
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

298
src/geometry.F90
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

46
src/global.F90
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
@ -184,7 +193,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
@ -213,9 +222,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
@ -248,7 +257,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.
@ -459,6 +467,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)
@ -471,6 +480,7 @@ contains
end do
deallocate(nuclides)
end if
if (allocated(libraries)) deallocate(libraries)
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
@ -479,7 +489,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)
@ -494,21 +503,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)
@ -534,6 +546,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()
@ -561,4 +574,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

35
src/initialize.F90
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)

171
src/input_xml.F90
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
@ -4807,9 +4802,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
@ -4912,9 +4905,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?
@ -5155,126 +5146,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
real(8) :: temp_frac ! temporary value for sorting
logical :: found
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
type(VectorReal) :: sab_fracs
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)
call sab_fracs % push_back(mat % 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(mat % &
sab_names(k)) // " did not match any nuclide on material " &
// trim(to_str(mat % 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(mat % nuclide(i_sab_nuclides % data(j))) % name) &
// " in material " // trim(to_str(mat % 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(mat % i_sab_tables)
deallocate(mat % sab_fracs)
m = i_sab_tables % size()
allocate(mat % i_sab_tables(m))
allocate(mat % i_sab_nuclides(m))
allocate(mat % sab_fracs(m))
mat % i_sab_tables(:) = i_sab_tables % data(1:m)
mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m)
mat % 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(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)
temp_frac = 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)
mat % sab_fracs(m) = mat % 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
mat % i_sab_nuclides(m) = temp_nuclide
mat % i_sab_tables(m) = temp_table
mat % 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(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(:)
@ -5290,9 +5162,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)
@ -5313,7 +5182,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)
@ -5343,6 +5213,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
@ -5365,7 +5242,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)
@ -5373,10 +5250,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)

29
src/main.F90
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

164
src/material_header.F90
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
!===============================================================================
@ -39,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

40
src/mesh.F90
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

1
src/message_passing.F90
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

3
src/mgxs_data.F90
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

57
src/nuclide_header.F90
View file

@ -97,6 +97,7 @@ 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
@ -170,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
@ -235,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
@ -673,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

255
src/output.F90
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
@ -744,20 +769,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
@ -889,24 +914,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, &
@ -918,11 +946,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
@ -932,11 +959,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
@ -974,11 +1002,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)
@ -1049,104 +1081,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

1
src/particle_header.F90
View file

@ -45,6 +45,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
! Energy Data

7
src/particle_restart.F90
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
!===============================================================================

12
src/plot.F90
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

2
src/relaxng/settings.rnc
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 }? &

8
src/relaxng/settings.rng
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"/>

16
src/sab_header.F90
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

168
src/simulation.F90
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

2
src/source.F90
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

31
src/state_point.F90
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

143
src/tally.F90
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
!===============================================================================
@ -4198,11 +4182,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
@ -4252,7 +4236,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
@ -4262,73 +4246,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
@ -4366,45 +4332,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
!===============================================================================

6
src/trigger.F90
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

11
src/volume_calc.F90
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