diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b07b3c482..c21b44cc16 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 + $ + ${CMAKE_CURRENT_SOURCE_DIR}/openmc/capi/_$ + COMMENT "Copying libopenmc to Python module directory") + #=============================================================================== # Install executable, scripts, manpage, license #=============================================================================== diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst new file mode 100644 index 0000000000..1e30eacf18 --- /dev/null +++ b/docs/source/capi/index.rst @@ -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 diff --git a/docs/source/conf.py b/docs/source/conf.py index 59928c5e15..95930fedba 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -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 diff --git a/docs/source/index.rst b/docs/source/index.rst index ce4a5f6f86..00b09db9e7 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -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 diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index e9bb5ef535..0aa058cdba 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -661,6 +661,17 @@ methods like "nearest" and "interpolation" in the resolved resonance range. *Default*: False +------------------------------- +```` Element +------------------------------- + +The ```` 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: ----------------------------------- diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst new file mode 100644 index 0000000000..6e4e860fad --- /dev/null +++ b/docs/source/pythonapi/capi.rst @@ -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 diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 5492f362df..bb99dd470f 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -47,5 +47,6 @@ Modules mgxs model data + capi examples openmoc diff --git a/examples/xml/basic/tallies.xml b/examples/xml/basic/tallies.xml index 26e61995da..a125e9ca2a 100644 --- a/examples/xml/basic/tallies.xml +++ b/examples/xml/basic/tallies.xml @@ -1,21 +1,30 @@ + + 100 + + + + 0 20.0e6 + + + + 0 20.0e6 + + - + 1 total scatter nu-scatter absorption fission nu-fission - - + 1 2 total scatter nu-scatter absorption fission nu-fission - - - + 1 2 3 scatter nu-scatter nu-fission diff --git a/examples/xml/lattice/nested/tallies.xml b/examples/xml/lattice/nested/tallies.xml index 89c0774f15..d342174a99 100644 --- a/examples/xml/lattice/nested/tallies.xml +++ b/examples/xml/lattice/nested/tallies.xml @@ -8,8 +8,12 @@ 1.0 1.0 + + 1 + + - + 1 total diff --git a/examples/xml/lattice/simple/tallies.xml b/examples/xml/lattice/simple/tallies.xml index 89c0774f15..d342174a99 100644 --- a/examples/xml/lattice/simple/tallies.xml +++ b/examples/xml/lattice/simple/tallies.xml @@ -8,8 +8,12 @@ 1.0 1.0 + + 1 + + - + 1 total diff --git a/examples/xml/pincell/tallies.xml b/examples/xml/pincell/tallies.xml index f740530cf7..642bdd7cd3 100644 --- a/examples/xml/pincell/tallies.xml +++ b/examples/xml/pincell/tallies.xml @@ -7,9 +7,16 @@ 0.62992 0.62992 1.e50 + + 1 + + + + 0. 4. 20.0e6 + + - - + 1 2 flux fission nu-fission diff --git a/examples/xml/pincell_multigroup/tallies.xml b/examples/xml/pincell_multigroup/tallies.xml index 1591a556b4..d84e129f2f 100644 --- a/examples/xml/pincell_multigroup/tallies.xml +++ b/examples/xml/pincell_multigroup/tallies.xml @@ -5,9 +5,14 @@ -0.63 -0.63 -1e+50 0.63 0.63 1e+50 + + 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + + + 1 + - - + 1 2 flux fission nu-fission diff --git a/openmc/capi/__init__.py b/openmc/capi/__init__.py new file mode 100644 index 0000000000..355304a310 --- /dev/null +++ b/openmc/capi/__init__.py @@ -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 * diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py new file mode 100644 index 0000000000..f519d2f6dc --- /dev/null +++ b/openmc/capi/cell.py @@ -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() diff --git a/openmc/capi/core.py b/openmc/capi/core.py new file mode 100644 index 0000000000..fe5416dd61 --- /dev/null +++ b/openmc/capi/core.py @@ -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() diff --git a/openmc/capi/error.py b/openmc/capi/error.py new file mode 100644 index 0000000000..2058ceb204 --- /dev/null +++ b/openmc/capi/error.py @@ -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)) diff --git a/openmc/capi/material.py b/openmc/capi/material.py new file mode 100644 index 0000000000..d07d7aa8a9 --- /dev/null +++ b/openmc/capi/material.py @@ -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() diff --git a/openmc/capi/nuclide.py b/openmc/capi/nuclide.py new file mode 100644 index 0000000000..4c6248d6be --- /dev/null +++ b/openmc/capi/nuclide.py @@ -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() diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py new file mode 100644 index 0000000000..611f1ce879 --- /dev/null +++ b/openmc/capi/tally.py @@ -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() diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 03062ee0e2..fc0d7ef07f 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -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 diff --git a/openmc/settings.py b/openmc/settings.py index b1ff4c0c65..3efdf1ab7a 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -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) diff --git a/setup.py b/setup.py index e7eb9281ca..67bc734035 100755 --- a/setup.py +++ b/setup.py @@ -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 diff --git a/src/api.F90 b/src/api.F90 new file mode 100644 index 0000000000..d45360d284 --- /dev/null +++ b/src/api.F90 @@ -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 diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 0d4a9f4d35..e87bff930c 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -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) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 68955f2dfc..83cb4ea257 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 86169d98d4..accf80b310 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -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 diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 deleted file mode 100644 index 47408943e3..0000000000 --- a/src/energy_grid.F90 +++ /dev/null @@ -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 diff --git a/src/finalize.F90 b/src/finalize.F90 deleted file mode 100644 index 0f1af23560..0000000000 --- a/src/finalize.F90 +++ /dev/null @@ -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 diff --git a/src/geometry.F90 b/src/geometry.F90 index f80bd7eb70..4444562321 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 1d45bcbe8d..2f47a9dce0 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index e9718fee10..1e29c638c8 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index aa41dc4b5f..2fc6ef9188 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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) diff --git a/src/main.F90 b/src/main.F90 index 1cd3a9e547..fadf9f3bbd 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -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 diff --git a/src/material_header.F90 b/src/material_header.F90 index 17fe8d2482..5d45c40d60 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/mesh.F90 b/src/mesh.F90 index b778cbdf06..c25e3f7a6a 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -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 diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 6c13b00e88..0f2b94d1af 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -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 diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index a43d45c0e1..ad20aca83b 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 5c5131abc7..bf1ea97b28 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index 10d4f0511c..9835225e06 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/particle_header.F90 b/src/particle_header.F90 index a7309393b2..f28e549761 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -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 diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 262cef8021..4438c48c00 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -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 !=============================================================================== diff --git a/src/plot.F90 b/src/plot.F90 index 6a3a191703..067d0c56fc 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -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 diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index a0f100cd85..23a0cfd4b7 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -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 }? & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 2a9a170f9a..e75dc7c15e 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -517,6 +517,14 @@ + + + + + + + + diff --git a/src/sab_header.F90 b/src/sab_header.F90 index 1d64711b8e..e0b79f2485 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -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 diff --git a/src/simulation.F90 b/src/simulation.F90 index 0cdb4c66d0..7427443bbe 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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 diff --git a/src/source.F90 b/src/source.F90 index 81931245a6..517715b563 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index e4c9a43ae7..e2a476e7e9 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -11,11 +11,12 @@ module state_point ! intervals, using the 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 diff --git a/src/tally.F90 b/src/tally.F90 index 0ee041cbe4..cf62412f03 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 !=============================================================================== diff --git a/src/trigger.F90 b/src/trigger.F90 index 9eb39e9cc2..dc0a10b153 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -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 diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index d63955017b..f8cf4085a5 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -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