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

Merging with the upstream, especially including hdf5 C++ changes

Browse source
This commit is contained in:
Adam G Nelson 2018-05-06 07:16:21 -04:00
commit 2000b4f998
88 changed files with 2695 additions and 3020 deletions
Download patch file Download diff file Expand all files Collapse all files

36
CMakeLists.txt
View file

@ -5,14 +5,10 @@ project(openmc Fortran C CXX)
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include)
# Set module path
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules)
# Make sure Fortran module directory is included when building
include_directories(${CMAKE_BINARY_DIR}/include)
#===============================================================================
# Architecture specific definitions
#===============================================================================
@ -73,7 +69,7 @@ if(NOT DEFINED HDF5_PREFER_PARALLEL)
endif()
endif()
find_package(HDF5 COMPONENTS Fortran_HL)
find_package(HDF5 COMPONENTS HL)
if(NOT HDF5_FOUND)
message(FATAL_ERROR "Could not find HDF5")
endif()
@ -436,20 +432,24 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/trigger_header.F90
)
set(LIBOPENMC_CXX_SRC
src/error.h
src/hdf5_interface.h
src/initialize.cpp
src/finalize.cpp
src/hdf5_interface.cpp
src/math_functions.h
src/math_functions.cpp
src/message_passing.cpp
src/plot.cpp
src/random_lcg.cpp
src/random_lcg.h
src/simulation.cpp
src/state_point.cpp
src/surface.cpp
src/surface.h
src/xml_interface.h
src/pugixml/pugixml.cpp
src/pugixml/pugixml.hpp)
src/xml_interface.cpp
src/pugixml/pugixml.cpp)
add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc)
add_executable(${program} src/main.F90)
set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc
PUBLIC_HEADER include/openmc.h)
add_executable(${program} src/main.cpp)
#===============================================================================
# Add compiler/linker flags
@ -458,11 +458,11 @@ add_executable(${program} src/main.F90)
set_property(TARGET ${program} libopenmc pugixml_fortran
PROPERTY LINKER_LANGUAGE Fortran)
target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS})
target_include_directories(libopenmc PUBLIC include ${HDF5_INCLUDE_DIRS})
# The executable and the faddeeva package use only one language. They can be
# set via target_compile_options which accepts a list.
target_compile_options(${program} PUBLIC ${f90flags})
target_compile_options(${program} PUBLIC ${cxxflags})
target_compile_options(faddeeva PRIVATE ${cflags})
# The libopenmc library has both F90 and C++ so the compile flags must be set
@ -503,7 +503,9 @@ add_custom_command(TARGET libopenmc POST_BUILD
install(TARGETS ${program} libopenmc
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include
)
install(DIRECTORY src/relaxng DESTINATION share/openmc)
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)

56
docs/source/capi/index.rst
View file

@ -46,10 +46,13 @@ Type Definitions
Functions
---------
.. c:function:: void openmc_calculate_volumes()
.. c:function:: int openmc_calculate_volumes()
Run a stochastic volume calculation
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
Get the fill for a cell
@ -192,11 +195,14 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: void openmc_finalize()
.. c:function:: int openmc_finalize()
Finalize a simulation
.. c:function:: void openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance)
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance)
Determine the ID of the cell/material containing a given point
@ -207,6 +213,8 @@ Functions
occurs, the ID is -1.
:param int32_t* instance: If a cell is repeated in the geometry, the instance
of the cell that was found and zero otherwise.
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_cell_index(int32_t id, int32_t* index)
@ -266,17 +274,24 @@ Functions
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_hard_reset()
.. c:function:: int openmc_hard_reset()
Reset tallies, timers, and pseudo-random number generator state
.. c:function:: void openmc_init(const int* intracomm)
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_init(int argc, char** argv, const void* intracomm)
Initialize OpenMC
:param int argc: Number of command-line arguments (including command)
:param char** argv: Command-line arguments
:param intracomm: MPI intracommunicator. If MPI is not being used, a null
pointer should be passed.
:type intracomm: const int*
:type intracomm: const void*
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_load_nuclide(char name[])
@ -394,26 +409,41 @@ Functions
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_plot_geometry()
.. c:function:: int openmc_plot_geometry()
Run plotting mode.
.. c:function:: void openmc_reset()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_reset()
Resets all tally scores
.. c:function:: void openmc_run()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_run()
Run a simulation
.. c:function:: void openmc_simulation_finalize()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_simulation_finalize()
Finalize a simulation.
.. c:function:: void openmc_simulation_init()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_simulation_init()
Initialize a simulation. Must be called after openmc_init().
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_source_bank(struct Bank** ptr, int64_t* n)
Return a pointer to the source bank array.
@ -433,13 +463,15 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: void openmc_statepoint_write(const char filename[])
.. c:function:: int openmc_statepoint_write(const char filename[])
Write a statepoint file
:param filename: Name of file to create. If a null pointer is passed, a
filename is assigned automatically.
:type filename: const char[]
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_tally_get_id(int32_t index, int32_t* id)

16
docs/source/io_formats/source.rst
View file

@ -8,13 +8,13 @@ Normally, source data is stored in a state point file. However, it is possible
to request that the source be written separately, in which case the format used
is that documented here.
**/filetype** (*char[]*)
**/**
String indicating the type of file.
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
**/source_bank** (Compound type)
Source bank information for each particle. The compound type has fields
``wgt``, ``xyz``, ``uvw``, ``E``, and ``delayed_group``, which
represent the weight, position, direction, energy, energy group, and
delayed_group of the source particle, respectively.
:Datasets:
- **source_bank** (Compound type) -- Source bank information for each
particle. The compound type has fields ``wgt``, ``xyz``, ``uvw``,
``E``, and ``delayed_group``, which represent the weight, position,
direction, energy, energy group, and delayed_group of the source
particle, respectively.

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

@ -141,16 +141,10 @@ Prerequisites
recommend that your HDF5 installation be built with parallel I/O
features. An example of configuring HDF5_ is listed below::
FC=mpifort ./configure --enable-fortran --enable-parallel
FC=mpifort ./configure --enable-parallel
You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial.
.. important::
If you are building HDF5 version 1.8.x or earlier, you must include
``--enable-fortran2003`` as well when configuring HDF5 or else OpenMC
will not be able to compile.
.. admonition:: Optional
:class: note

45
include/openmc.h
View file

@ -16,7 +16,7 @@ extern "C" {
int delayed_group;
};
void openmc_calculate_volumes();
int openmc_calculate_volumes();
int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
int openmc_cell_get_id(int32_t index, int32_t* id);
int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices);
@ -34,7 +34,7 @@ extern "C" {
int openmc_filter_get_type(int32_t index, const char** type);
int openmc_filter_set_id(int32_t index, int32_t id);
int openmc_filter_set_type(int32_t index, const char* type);
void openmc_finalize();
int openmc_finalize();
int openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance);
int openmc_get_cell_index(int32_t id, int32_t* index);
int openmc_get_filter_index(int32_t id, int32_t* index);
@ -45,8 +45,9 @@ extern "C" {
int openmc_get_nuclide_index(const char name[], int* index);
int64_t openmc_get_seed();
int openmc_get_tally_index(int32_t id, int32_t* index);
void openmc_hard_reset();
void openmc_init(const int* intracomm);
int openmc_hard_reset();
int openmc_init(int argc, char* argv[], const void* intracomm);
int openmc_init_f(const int* intracomm);
int openmc_legendre_filter_get_order(int32_t index, int* order);
int openmc_legendre_filter_set_order(int32_t index, int order);
int openmc_load_nuclide(char name[]);
@ -70,12 +71,13 @@ extern "C" {
int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, char** name);
void openmc_plot_geometry();
void openmc_reset();
int openmc_particle_restart();
int openmc_plot_geometry();
int openmc_reset();
int openmc_run();
void openmc_set_seed(int64_t new_seed);
void openmc_simulation_finalize();
void openmc_simulation_init();
int openmc_simulation_finalize();
int openmc_simulation_init();
int openmc_source_bank(struct Bank** ptr, int64_t* n);
int openmc_source_set_strength(int32_t index, double strength);
int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
@ -87,7 +89,7 @@ extern "C" {
int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
void openmc_statepoint_write(const char filename[]);
int openmc_statepoint_write(const char filename[]);
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_id(int32_t index, int32_t* id);
int openmc_tally_get_filters(int32_t index, int32_t** indices, int* n);
@ -122,8 +124,8 @@ extern "C" {
// Global variables
extern char openmc_err_msg[256];
extern double keff;
extern double keff_std;
extern double openmc_keff;
extern double openmc_keff_std;
extern int32_t n_batches;
extern int32_t n_cells;
extern int32_t n_filters;
@ -139,9 +141,24 @@ extern "C" {
extern int32_t n_surfaces;
extern int32_t n_tallies;
extern int32_t n_universes;
extern int run_mode;
extern bool simulation_initialized;
extern int verbosity;
extern int openmc_run_mode;
extern bool openmc_simulation_initialized;
extern int openmc_verbosity;
// Variables that are shared by necessity (can be removed from public header
// later)
extern bool openmc_master;
extern int openmc_n_procs;
extern int openmc_n_threads;
extern int openmc_rank;
extern int64_t openmc_work;
// Run modes
constexpr int RUN_MODE_FIXEDSOURCE {1};
constexpr int RUN_MODE_EIGENVALUE {2};
constexpr int RUN_MODE_PLOTTING {3};
constexpr int RUN_MODE_PARTICLE {4};
constexpr int RUN_MODE_VOLUME {5};
#ifdef __cplusplus
}

75
openmc/capi/core.py
View file

@ -1,6 +1,6 @@
from contextlib import contextmanager
from ctypes import (CDLL, c_int, c_int32, c_int64, c_double, c_char_p,
POINTER, Structure)
from ctypes import (CDLL, c_int, c_int32, c_int64, c_double, c_char_p, c_char,
POINTER, Structure, c_void_p, create_string_buffer)
from warnings import warn
import numpy as np
@ -20,32 +20,41 @@ class _Bank(Structure):
('delayed_group', c_int)]
_dll.openmc_calculate_volumes.restype = None
_dll.openmc_finalize.restype = None
_dll.openmc_calculate_volumes.restype = c_int
_dll.openmc_calculate_volumes.errcheck = _error_handler
_dll.openmc_finalize.restype = c_int
_dll.openmc_finalize.errcheck = _error_handler
_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_hard_reset.restype = c_int
_dll.openmc_hard_reset.errcheck = _error_handler
_dll.openmc_init.argtypes = [c_int, POINTER(POINTER(c_char)), c_void_p]
_dll.openmc_init.restype = c_int
_dll.openmc_init.errcheck = _error_handler
_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_next_batch.argtypes = [POINTER(c_int)]
_dll.openmc_next_batch.restype = c_int
_dll.openmc_next_batch.errcheck = _error_handler
_dll.openmc_plot_geometry.restype = None
_dll.openmc_plot_geometry.restype = c_int
_dll.openmc_plot_geometry.restype = _error_handler
_dll.openmc_run.restype = c_int
_dll.openmc_run.errcheck = _error_handler
_dll.openmc_reset.restype = None
_dll.openmc_reset.restype = c_int
_dll.openmc_reset.errcheck = _error_handler
_dll.openmc_source_bank.argtypes = [POINTER(POINTER(_Bank)), POINTER(c_int64)]
_dll.openmc_source_bank.restype = c_int
_dll.openmc_source_bank.errcheck = _error_handler
_dll.openmc_simulation_init.restype = None
_dll.openmc_simulation_finalize.restype = None
_dll.openmc_simulation_init.restype = c_int
_dll.openmc_simulation_init.errcheck = _error_handler
_dll.openmc_simulation_finalize.restype = c_int
_dll.openmc_simulation_finalize.errcheck = _error_handler
_dll.openmc_statepoint_write.argtypes = [POINTER(c_char_p)]
_dll.openmc_statepoint_write.restype = None
_dll.openmc_statepoint_write.restype = c_int
_dll.openmc_statepoint_write.errcheck = _error_handler
def calculate_volumes():
@ -106,25 +115,42 @@ def hard_reset():
_dll.openmc_hard_reset()
def init(intracomm=None):
def init(args=None, intracomm=None):
"""Initialize OpenMC
Parameters
----------
args : list of str
Command-line arguments
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))
if args is not None:
args = ['openmc'] + list(args)
argc = len(args)
# Create the argv array. Note that it is actually expected to be of
# length argc + 1 with the final item being a null pointer.
argv = (POINTER(c_char) * (argc + 1))()
for i, arg in enumerate(args):
argv[i] = create_string_buffer(arg.encode())
else:
_dll.openmc_init(None)
argc = 0
argv = None
if intracomm is not None:
# If an mpi4py communicator was passed, convert it to void* to be passed
# to openmc_init
try:
from mpi4py import MPI
except ImportError:
intracomm = None
else:
address = MPI._addressof(intracomm)
intracomm = c_void_p(address)
_dll.openmc_init(argc, argv, intracomm)
def iter_batches():
@ -178,8 +204,9 @@ def keff():
return tuple(k)
else:
# Otherwise, return the tracklength estimator
mean = c_double.in_dll(_dll, 'keff').value
std_dev = c_double.in_dll(_dll, 'keff_std').value if n > 1 else np.inf
mean = c_double.in_dll(_dll, 'openmc_keff').value
std_dev = c_double.in_dll(_dll, 'openmc_keff_std').value \
if n > 1 else np.inf
return (mean, std_dev)

6
openmc/capi/settings.py
View file

@ -20,11 +20,11 @@ class _Settings(object):
generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch')
inactive = _DLLGlobal(c_int32, 'n_inactive')
particles = _DLLGlobal(c_int64, 'n_particles')
verbosity = _DLLGlobal(c_int, 'verbosity')
verbosity = _DLLGlobal(c_int, 'openmc_verbosity')
@property
def run_mode(self):
i = c_int.in_dll(_dll, 'run_mode').value
i = c_int.in_dll(_dll, 'openmc_run_mode').value
try:
return _RUN_MODES[i]
except KeyError:
@ -32,7 +32,7 @@ class _Settings(object):
@run_mode.setter
def run_mode(self, mode):
current_idx = c_int.in_dll(_dll, 'run_mode')
current_idx = c_int.in_dll(_dll, 'openmc_run_mode')
for idx, mode_value in _RUN_MODES.items():
if mode_value == mode:
current_idx.value = idx

2
openmc/deplete/operator.py
View file

@ -268,7 +268,7 @@ class Operator(TransportOperator):
# Initialize OpenMC library
comm.barrier()
openmc.capi.init(comm)
openmc.capi.init(intracomm=comm)
# Generate tallies in memory
self._generate_tallies()

2
openmc/statepoint.py
View file

@ -150,6 +150,8 @@ class StatePoint(object):
def __exit__(self, *exc):
self._f.close()
if self._summary is not None:
self._summary._f.close()
@property
def cmfd_on(self):

4
src/angle_distribution.F90
View file

@ -1,12 +1,10 @@
module angle_distribution
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, ONE, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer, Tabular
use hdf5_interface, only: read_attribute, get_shape, read_dataset, &
open_dataset, close_dataset
open_dataset, close_dataset, HID_T, HSIZE_T
use random_lcg, only: prn
implicit none

2
src/angleenergy_header.F90
View file

@ -1,6 +1,6 @@
module angleenergy_header
use hdf5, only: HID_T
use hdf5_interface, only: HID_T
!===============================================================================
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy

42
src/api.F90
View file

@ -2,8 +2,6 @@ module openmc_api
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, h5tclose_f, h5close_f
use bank_header, only: openmc_source_bank
use constants, only: K_BOLTZMANN
use eigenvalue, only: k_sum, openmc_get_keff
@ -16,7 +14,7 @@ module openmc_api
use mesh_header
use message_passing
use nuclide_header
use initialize, only: openmc_init
use initialize, only: openmc_init_f
use particle_header, only: Particle
use plot, only: openmc_plot_geometry
use random_lcg, only: openmc_get_seed, openmc_set_seed
@ -65,7 +63,7 @@ module openmc_api
public :: openmc_get_tally_index
public :: openmc_global_tallies
public :: openmc_hard_reset
public :: openmc_init
public :: openmc_init_f
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
@ -81,7 +79,6 @@ module openmc_api
public :: openmc_nuclide_name
public :: openmc_plot_geometry
public :: openmc_reset
public :: openmc_run
public :: openmc_set_seed
public :: openmc_simulation_finalize
public :: openmc_simulation_init
@ -117,12 +114,16 @@ contains
! variables
!===============================================================================
subroutine openmc_finalize() bind(C)
function openmc_finalize() result(err) bind(C)
integer(C_INT) :: err
integer :: err
interface
subroutine openmc_free_bank() bind(C)
end subroutine openmc_free_bank
end interface
! Clear results
call openmc_reset()
err = openmc_reset()
! Reset global variables
assume_separate = .false.
@ -182,18 +183,14 @@ contains
! Deallocate arrays
call free_memory()
! Release compound datatypes
call h5tclose_f(hdf5_bank_t, err)
! Close FORTRAN interface.
call h5close_f(err)
err = 0
#ifdef OPENMC_MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, err)
call openmc_free_bank()
#endif
end subroutine openmc_finalize
end function openmc_finalize
!===============================================================================
! OPENMC_FIND determines the ID or a cell or material at a given point in space
@ -244,9 +241,11 @@ contains
! generator state
!===============================================================================
subroutine openmc_hard_reset() bind(C)
function openmc_hard_reset() result(err) bind(C)
integer(C_INT) :: err
! Reset all tallies and timers
call openmc_reset()
err = openmc_reset()
! Reset total generations and keff guess
keff = ONE
@ -254,13 +253,15 @@ contains
! Reset the random number generator state
call openmc_set_seed(DEFAULT_SEED)
end subroutine openmc_hard_reset
end function openmc_hard_reset
!===============================================================================
! OPENMC_RESET resets tallies and timers
!===============================================================================
subroutine openmc_reset() bind(C)
function openmc_reset() result(err) bind(C)
integer(C_INT) :: err
integer :: i
if (allocated(tallies)) then
@ -308,7 +309,8 @@ contains
call time_transport % reset()
call time_finalize % reset()
end subroutine openmc_reset
err = 0
end function openmc_reset
!===============================================================================
! FREE_MEMORY deallocates and clears all global allocatable arrays in the

2
src/endf_header.F90
View file

@ -1,7 +1,5 @@
module endf_header
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, &
LOG_LINEAR, LOG_LOG

2
src/energy_distribution.F90
View file

@ -1,7 +1,5 @@
module energy_distribution
use hdf5
use algorithm, only: binary_search
use constants, only: ZERO, ONE, HALF, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tabulated1D

8
src/error.F90
View file

@ -111,6 +111,14 @@ contains
end subroutine warning
subroutine warning_from_c(message, message_len) bind(C)
integer(C_INT), intent(in), value :: message_len
character(kind=C_CHAR), intent(in) :: message(message_len)
character(message_len+1) :: message_out
write(message_out, *) message
call warning(message_out)
end subroutine
!===============================================================================
! FATAL_ERROR alerts the user that an error has been encountered and displays a
! message about the particular problem. Errors are considered 'fatal' and hence

23
src/error.h
View file

@ -9,25 +9,38 @@
namespace openmc {
extern "C" void fatal_error_from_c(const char *message, int message_len);
extern "C" void fatal_error_from_c(const char* message, int message_len);
extern "C" void warning_from_c(const char* message, int message_len);
inline
void fatal_error(const char *message)
{
fatal_error_from_c(message, strlen(message));
}
inline
void fatal_error(const std::string &message)
{
fatal_error_from_c(message.c_str(), message.length());
}
inline
void fatal_error(const std::stringstream &message)
{
std::string out {message.str()};
fatal_error_from_c(out.c_str(), out.length());
fatal_error(message.str());
}
inline
void warning(const std::string& message)
{
warning_from_c(message.c_str(), message.length());
}
inline
void warning(const std::stringstream& message)
{
warning(message.str());
}
} // namespace openmc

10
src/finalize.cpp Normal file
View file

@ -0,0 +1,10 @@
#include "finalize.h"
#include "message_passing.h"
void openmc_free_bank()
{
#ifdef OPENMC_MPI
MPI_Type_free(&openmc::mpi::bank);
#endif
}

6
src/finalize.h Normal file
View file

@ -0,0 +1,6 @@
#ifndef FINALIZE_H
#define FINALIZE_H
extern "C" void openmc_free_bank();
#endif // FINALIZE_H

3
src/geometry.F90
View file

@ -355,9 +355,10 @@ contains
else
! Particle is outside the lattice.
if (lat % outer == NO_OUTER_UNIVERSE) then
call p % mark_as_lost("Particle " // trim(to_str(p %id)) &
call warning("Particle " // trim(to_str(p %id)) &
// " is outside lattice " // trim(to_str(lat % id)) &
// " but the lattice has no defined outer universe.")
found = .false.
return
else
p % coord(j + 1) % universe = lat % outer

2803
src/hdf5_interface.F90
View file

File diff suppressed because it is too large Load diff

658
src/hdf5_interface.cpp Normal file
View file

@ -0,0 +1,658 @@
#include "hdf5_interface.h"
#include <array>
#include <cstring>
#include <sstream>
#include <string>
#include "hdf5.h"
#include "hdf5_hl.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#include "message_passing.h"
#endif
#include "error.h"
namespace openmc {
bool
attribute_exists(hid_t obj_id, const char* name)
{
htri_t out = H5Aexists_by_name(obj_id, ".", name, H5P_DEFAULT);
return out > 0;
}
size_t
attribute_typesize(hid_t obj_id, const char* name)
{
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
hid_t filetype = H5Aget_type(attr);
size_t n = H5Tget_size(filetype);
H5Tclose(filetype);
H5Aclose(attr);
return n;
}
void
get_shape(hid_t obj_id, hsize_t* dims)
{
auto type = H5Iget_type(obj_id);
hid_t dspace;
if (type == H5I_DATASET) {
dspace = H5Dget_space(obj_id);
} else if (type == H5I_ATTR) {
dspace = H5Aget_space(obj_id);
}
H5Sget_simple_extent_dims(dspace, dims, nullptr);
H5Sclose(dspace);
}
void
get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims)
{
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
hid_t dspace = H5Aget_space(attr);
H5Sget_simple_extent_dims(dspace, dims, nullptr);
H5Sclose(dspace);
H5Aclose(attr);
}
hid_t
create_group(hid_t parent_id, char const *name)
{
hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (out < 0) {
std::stringstream err_msg;
err_msg << "Failed to create HDF5 group \"" << name << "\"";
fatal_error(err_msg);
}
return out;
}
hid_t
create_group(hid_t parent_id, const std::string &name)
{
return create_group(parent_id, name.c_str());
}
void
close_dataset(hid_t dataset_id)
{
if (H5Dclose(dataset_id) < 0) fatal_error("Failed to close dataset");
}
void
close_group(hid_t group_id)
{
if (H5Gclose(group_id) < 0) fatal_error("Failed to close group");
}
int
dataset_ndims(hid_t dset)
{
hid_t dspace = H5Dget_space(dset);
int ndims = H5Sget_simple_extent_ndims(dspace);
H5Sclose(dspace);
return ndims;
}
size_t
dataset_typesize(hid_t dset)
{
hid_t filetype = H5Dget_type(dset);
size_t n = H5Tget_size(filetype);
H5Tclose(filetype);
return n;
}
hid_t
file_open(const char* filename, char mode, bool parallel)
{
bool create;
unsigned int flags;
switch (mode) {
case 'r':
case 'a':
create = false;
flags = (mode == 'r' ? H5F_ACC_RDONLY : H5F_ACC_RDWR);
break;
case 'w':
case 'x':
create = true;
flags = (mode == 'x' ? H5F_ACC_EXCL : H5F_ACC_TRUNC);
break;
default:
std::stringstream err_msg;
err_msg << "Invalid file mode: " << mode;
fatal_error(err_msg);
}
hid_t plist = H5P_DEFAULT;
#ifdef PHDF5
if (parallel) {
// Setup file access property list with parallel I/O access
plist = H5Pcreate(H5P_FILE_ACCESS);
H5Pset_fapl_mpio(plist, openmc::mpi::intracomm, MPI_INFO_NULL);
}
#endif
// Open the file collectively
hid_t file_id;
if (create) {
file_id = H5Fcreate(filename, flags, H5P_DEFAULT, plist);
} else {
file_id = H5Fopen(filename, flags, plist);
}
if (file_id < 0) {
std::stringstream msg;
msg << "Failed to open HDF5 file with mode '" << mode << "': " << filename;
fatal_error(msg);
}
#ifdef PHDF5
// Close the property list
if (parallel) H5Pclose(plist);
#endif
return file_id;
}
hid_t
file_open(const std::string& filename, char mode, bool parallel=false)
{
file_open(filename.c_str(), mode, parallel);
}
void file_close(hid_t file_id)
{
H5Fclose(file_id);
}
void
get_name(hid_t obj_id, char* name)
{
size_t size = 1 + H5Iget_name(obj_id, nullptr, 0);
H5Iget_name(obj_id, name, size);
}
int get_num_datasets(hid_t group_id)
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get number of groups
H5O_info_t oinfo;
int ndatasets = 0;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type == H5O_TYPE_DATASET) ndatasets += 1;
}
return ndatasets;
}
int get_num_groups(hid_t group_id)
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get number of groups
H5O_info_t oinfo;
int ngroups = 0;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type == H5O_TYPE_GROUP) ngroups += 1;
}
return ngroups;
}
void
get_datasets(hid_t group_id, char* name[])
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get names
H5O_info_t oinfo;
hsize_t count = 0;
size_t size;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type != H5O_TYPE_DATASET) continue;
// Get size of name
size = 1 + H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC,
i, nullptr, 0, H5P_DEFAULT);
// Read name
H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i,
name[count], size, H5P_DEFAULT);
count += 1;
}
}
void
get_groups(hid_t group_id, char* name[])
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get names
H5O_info_t oinfo;
hsize_t count = 0;
size_t size;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type != H5O_TYPE_GROUP) continue;
// Get size of name
size = 1 + H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC,
i, nullptr, 0, H5P_DEFAULT);
// Read name
H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i,
name[count], size, H5P_DEFAULT);
count += 1;
}
}
bool
object_exists(hid_t object_id, const char* name)
{
htri_t out = H5LTpath_valid(object_id, name, true);
if (out < 0) {
std::stringstream err_msg;
err_msg << "Failed to check if object \"" << name << "\" exists.";
fatal_error(err_msg);
}
return (out > 0);
}
hid_t
open_dataset(hid_t group_id, const char* name)
{
if (object_exists(group_id, name)) {
return H5Dopen(group_id, name, H5P_DEFAULT);
} else {
std::stringstream err_msg;
err_msg << "Group \"" << name << "\" does not exist";
fatal_error(err_msg);
}
}
hid_t
open_group(hid_t group_id, const char* name)
{
if (object_exists(group_id, name)) {
return H5Gopen(group_id, name, H5P_DEFAULT);
} else {
std::stringstream err_msg;
err_msg << "Group \"" << name << "\" does not exist";
fatal_error(err_msg);
}
}
void
read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer)
{
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
H5Aread(attr, mem_type_id, buffer);
H5Aclose(attr);
}
void
read_attr_double(hid_t obj_id, const char* name, double* buffer)
{
read_attr(obj_id, name, H5T_NATIVE_DOUBLE, buffer);
}
void
read_attr_int(hid_t obj_id, const char* name, int* buffer)
{
read_attr(obj_id, name, H5T_NATIVE_INT, buffer);
}
void
read_attr_string(hid_t obj_id, const char* name, size_t slen, char* buffer)
{
// Create datatype for a string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, slen + 1);
// Read data into buffer
read_attr(obj_id, name, datatype, buffer);
// Free resources
H5Tclose(datatype);
}
void
read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
void* buffer, bool indep)
{
hid_t dset = obj_id;
if (name) dset = open_dataset(obj_id, name);
if (using_mpio_device(dset)) {
#ifdef PHDF5
// Set up collective vs independent I/O
auto data_xfer_mode = indep ? H5FD_MPIO_INDEPENDENT : H5FD_MPIO_COLLECTIVE;
// Create dataset transfer property list
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, data_xfer_mode);
// Read data
H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, plist, buffer);
H5Pclose(plist);
#endif
} else {
H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
}
if (name) H5Dclose(dset);
}
void
read_double(hid_t obj_id, const char* name, double* buffer, bool indep)
{
read_dataset(obj_id, name, H5T_NATIVE_DOUBLE, buffer, indep);
}
void
read_int(hid_t obj_id, const char* name, int* buffer, bool indep)
{
read_dataset(obj_id, name, H5T_NATIVE_INT, buffer, indep);
}
void
read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep)
{
read_dataset(obj_id, name, H5T_NATIVE_LLONG, buffer, indep);
}
void
read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep)
{
// Create datatype for a string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, slen + 1);
// Read data into buffer
read_dataset(obj_id, name, datatype, buffer, indep);
// Free resources
H5Tclose(datatype);
}
void
read_complex(hid_t obj_id, const char* name, double _Complex* buffer, bool indep)
{
// Create compound datatype for complex numbers
struct complex_t {
double re;
double im;
};
complex_t tmp;
hid_t complex_id = H5Tcreate(H5T_COMPOUND, sizeof tmp);
H5Tinsert(complex_id, "r", HOFFSET(complex_t, re), H5T_NATIVE_DOUBLE);
H5Tinsert(complex_id, "i", HOFFSET(complex_t, im), H5T_NATIVE_DOUBLE);
// Read data
read_dataset(obj_id, name, complex_id, buffer, indep);
// Free resources
H5Tclose(complex_id);
}
void
read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
{
// Create dataspace for hyperslab in memory
hsize_t dims[] {n_filter, n_score, 3};
hsize_t start[] {0, 0, 1};
hsize_t count[] {n_filter, n_score, 2};
hid_t memspace = H5Screate_simple(3, dims, nullptr);
H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Create and write dataset
hid_t dset = H5Dopen(group_id, "results", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, memspace, H5S_ALL, H5P_DEFAULT, results);
// Free resources
H5Dclose(dset);
H5Sclose(memspace);
}
void
write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer)
{
// If array is given, create a simple dataspace. Otherwise, create a scalar
// datascape.
hid_t dspace;
if (ndim > 0) {
dspace = H5Screate_simple(ndim, dims, nullptr);
} else {
dspace = H5Screate(H5S_SCALAR);
}
// Create attribute and Write data
hid_t attr = H5Acreate(obj_id, name, mem_type_id, dspace,
H5P_DEFAULT, H5P_DEFAULT);
H5Awrite(attr, mem_type_id, buffer);
// Free resources
H5Aclose(attr);
H5Sclose(dspace);
}
void
write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
const double* buffer)
{
write_attr(obj_id, ndim, dims, name, H5T_NATIVE_DOUBLE, buffer);
}
void
write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
const int* buffer)
{
write_attr(obj_id, ndim, dims, name, H5T_NATIVE_INT, buffer);
}
void
write_attr_string(hid_t obj_id, const char* name, const char* buffer)
{
size_t n = strlen(buffer);
if (n > 0) {
// Set up appropriate datatype for a fixed-length string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, n);
write_attr(obj_id, 0, nullptr, name, datatype, buffer);
// Free resources
H5Tclose(datatype);
}
}
void
write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer, bool indep)
{
// If array is given, create a simple dataspace. Otherwise, create a scalar
// datascape.
hid_t dspace;
if (ndim > 0) {
dspace = H5Screate_simple(ndim, dims, nullptr);
} else {
dspace = H5Screate(H5S_SCALAR);
}
hid_t dset = H5Dcreate(group_id, name, mem_type_id, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (using_mpio_device(group_id)) {
#ifdef PHDF5
// Set up collective vs independent I/O
auto data_xfer_mode = indep ? H5FD_MPIO_INDEPENDENT : H5FD_MPIO_COLLECTIVE;
// Create dataset transfer property list
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, data_xfer_mode);
// Write data
H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, plist, buffer);
H5Pclose(plist);
#endif
} else {
H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
}
// Free resources
H5Dclose(dset);
H5Sclose(dspace);
}
void
write_double(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
const double* buffer, bool indep)
{
write_dataset(group_id, ndim, dims, name, H5T_NATIVE_DOUBLE, buffer, indep);
}
void
write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
const int* buffer, bool indep)
{
write_dataset(group_id, ndim, dims, name, H5T_NATIVE_INT, buffer, indep);
}
void
write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
const long long* buffer, bool indep)
{
write_dataset(group_id, ndim, dims, name, H5T_NATIVE_LLONG, buffer, indep);
}
void
write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
const char* name, const char* buffer, bool indep)
{
if (slen > 0) {
// Set up appropriate datatype for a fixed-length string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, slen);
write_dataset(group_id, ndim, dims, name, datatype, buffer, indep);
// Free resources
H5Tclose(datatype);
}
}
void
write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep)
{
write_string(group_id, 0, nullptr, buffer.length(), name, buffer.c_str(), indep);
}
void
write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results)
{
// Set dimensions of sum/sum_sq hyperslab to store
hsize_t count[] {n_filter, n_score, 2};
hid_t dspace = H5Screate_simple(3, count, nullptr);
// Set dimensions of results array
hsize_t dims[] {n_filter, n_score, 3};
hsize_t start[] {0, 0, 1};
hid_t memspace = H5Screate_simple(3, dims, nullptr);
H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Create and write dataset
hid_t dset = H5Dcreate(group_id, "results", H5T_NATIVE_DOUBLE, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dset, H5T_NATIVE_DOUBLE, memspace, H5S_ALL, H5P_DEFAULT, results);
// Free resources
H5Dclose(dset);
H5Sclose(memspace);
H5Sclose(dspace);
}
bool
using_mpio_device(hid_t obj_id)
{
// Determine file that this object is part of
hid_t file_id = H5Iget_file_id(obj_id);
// Get file access property list
hid_t fapl_id = H5Fget_access_plist(file_id);
// Get low-level driver identifier
hid_t driver = H5Pget_driver(fapl_id);
// Free resources
H5Pclose(fapl_id);
H5Fclose(file_id);
return driver == H5FD_MPIO;
}
} // namespace openmc

142
src/hdf5_interface.h
View file

@ -1,90 +1,96 @@
#ifndef HDF5_INTERFACE_H
#define HDF5_INTERFACE_H
#include <array>
#include <string.h>
#include "hdf5.h"
#include "hdf5_hl.h"
#include "error.h"
#include <array>
#include <string>
#include <sstream>
#include <complex.h>
namespace openmc {
extern "C" bool attribute_exists(hid_t obj_id, const char* name);
extern "C" size_t attribute_typesize(hid_t obj_id, const char* name);
extern "C" hid_t create_group(hid_t parent_id, const char* name);
hid_t create_group(hid_t parent_id, const std::string& name);
extern "C" void close_dataset(hid_t dataset_id);
extern "C" void close_group(hid_t group_id);
extern "C" int dataset_ndims(hid_t dset);
extern "C" size_t dataset_typesize(hid_t dset);
extern "C" hid_t file_open(const char* filename, char mode, bool parallel);
hid_t file_open(const std::string& filename, char mode, bool parallel);
extern "C" void file_close(hid_t file_id);
extern "C" void get_name(hid_t obj_id, char* name);
extern "C" int get_num_datasets(hid_t group_id);
extern "C" int get_num_groups(hid_t group_id);
extern "C" void get_datasets(hid_t group_id, char* name[]);
extern "C" void get_groups(hid_t group_id, char* name[]);
extern "C" void get_shape(hid_t obj_id, hsize_t* dims);
extern "C" void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims);
extern "C" bool object_exists(hid_t object_id, const char* name);
extern "C" hid_t open_dataset(hid_t group_id, const char* name);
extern "C" hid_t open_group(hid_t group_id, const char* name);
bool using_mpio_device(hid_t obj_id);
hid_t
create_group(hid_t parent_id, char const *name)
{
hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (out < 0) {
std::string err_msg{"Failed to create HDF5 group \""};
err_msg += name;
err_msg += "\"";
fatal_error(err_msg);
}
return out;
}
void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id,
const void* buffer);
extern "C" void read_attr_double(hid_t obj_id, const char* name, double* buffer);
extern "C" void read_attr_int(hid_t obj_id, const char* name, int* buffer);
extern "C" void read_attr_string(hid_t obj_id, const char* name, size_t slen,
char* buffer);
void read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
void* buffer, bool indep);
extern "C" void read_double(hid_t obj_id, const char* name, double* buffer,
bool indep);
extern "C" void read_int(hid_t obj_id, const char* name, int* buffer,
bool indep);
extern "C" void read_llong(hid_t obj_id, const char* name, long long* buffer,
bool indep);
extern "C" void read_string(hid_t obj_id, const char* name, size_t slen,
char* buffer, bool indep);
extern "C" void read_complex(hid_t obj_id, const char* name,
double _Complex* buffer, bool indep);
extern "C" void read_tally_results(hid_t group_id, hsize_t n_filter,
hsize_t n_score, double* results);
void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer);
extern "C" void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
const char* name, const double* buffer);
extern "C" void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
const char* name, const int* buffer);
extern "C" void write_attr_string(hid_t obj_id, const char* name, const char* buffer);
hid_t
create_group(hid_t parent_id, const std::string &name)
{
return create_group(parent_id, name.c_str());
}
void write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer, bool indep);
extern "C" void write_double(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const double* buffer, bool indep);
extern "C" void write_int(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const int* buffer, bool indep);
extern "C" void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const long long* buffer, bool indep);
extern "C" void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
const char* name, char const* buffer, bool indep);
void write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep);
void
close_group(hid_t group_id)
{
herr_t err = H5Gclose(group_id);
if (err < 0) {
fatal_error("Failed to close HDF5 group");
}
}
extern "C" void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
const double* results);
template<std::size_t array_len> void
write_double_1D(hid_t group_id, char const *name,
std::array<double, array_len> &buffer)
{
hsize_t dims[1]{array_len};
hid_t dataspace = H5Screate_simple(1, dims, NULL);
hid_t dataset = H5Dcreate(group_id, name, H5T_NATIVE_DOUBLE, dataspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&buffer[0]);
H5Sclose(dataspace);
H5Dclose(dataset);
}
void
write_string(hid_t group_id, char const *name, char const *buffer)
{
size_t buffer_len = strlen(buffer);
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, buffer_len);
hid_t dataspace = H5Screate(H5S_SCALAR);
hid_t dataset = H5Dcreate(group_id, name, datatype, dataspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
H5Tclose(datatype);
H5Sclose(dataspace);
H5Dclose(dataset);
}
void
write_string(hid_t group_id, char const *name, const std::string &buffer)
{
write_string(group_id, name, buffer.c_str());
hsize_t dims[1] {array_len};
write_dataset(group_id, 1, dims, name, H5T_NATIVE_DOUBLE,
buffer.data(), false);
}
} // namespace openmc

262
src/initialize.F90
View file

@ -1,37 +1,27 @@
module initialize
use, intrinsic :: ISO_C_BINDING, only: c_loc
use, intrinsic :: ISO_C_BINDING
use hdf5
#ifdef _OPENMP
use omp_lib
#endif
use bank_header, only: Bank
use constants
use set_header, only: SetInt
use error, only: fatal_error, warning, write_message
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
root_universe
use hdf5_interface, only: file_open, read_attribute, file_close, &
hdf5_bank_t, hdf5_integer8_t
use input_xml, only: read_input_xml
use material_header, only: Material
use message_passing
use mgxs_data, only: read_mgxs, create_macro_xs
use output, only: print_version, print_usage
use random_lcg, only: openmc_set_seed
use settings
#ifdef _OPENMP
use simulation_header, only: n_threads
#endif
use string, only: to_str, starts_with, ends_with, str_to_int
use tally_header, only: TallyObject
use tally_filter
use string, only: ends_with, to_f_string
use timer_header
implicit none
type(C_PTR), bind(C) :: openmc_path_input
type(C_PTR), bind(C) :: openmc_path_statepoint
type(C_PTR), bind(C) :: openmc_path_sourcepoint
type(C_PTR), bind(C) :: openmc_path_particle_restart
contains
!===============================================================================
@ -41,8 +31,9 @@ contains
! setting up timers, etc.
!===============================================================================
subroutine openmc_init(intracomm) bind(C)
function openmc_init_f(intracomm) result(err) bind(C)
integer, intent(in), optional :: intracomm ! MPI intracommunicator
integer(C_INT) :: err
#ifdef _OPENMP
character(MAX_WORD_LEN) :: envvar
@ -86,9 +77,6 @@ contains
end if
#endif
! Initialize HDF5 interface
call hdf5_initialize()
! Read command line arguments
call read_command_line()
@ -105,7 +93,8 @@ contains
! Stop initialization timer
call time_initialize%stop()
end subroutine openmc_init
err = 0
end function openmc_init_f
#ifdef OPENMC_MPI
!===============================================================================
@ -129,27 +118,11 @@ contains
integer :: bank_types(5) ! Datatypes
#endif
integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements
logical :: init_called
type(Bank) :: b
! Indicate that MPI is turned on
mpi_enabled = .true.
! Initialize MPI
call MPI_INITIALIZED(init_called, mpi_err)
if (.not. init_called) call MPI_INIT(mpi_err)
! Determine number of processors and rank of each processor
mpi_intracomm = intracomm
call MPI_COMM_SIZE(mpi_intracomm, n_procs, mpi_err)
call MPI_COMM_RANK(mpi_intracomm, rank, mpi_err)
! Determine master
if (rank == 0) then
master = .true.
else
master = .false.
end if
! ==========================================================================
! CREATE MPI_BANK TYPE
@ -174,203 +147,46 @@ contains
end subroutine initialize_mpi
#endif
!===============================================================================
! HDF5_INITIALIZE
!===============================================================================
subroutine hdf5_initialize()
type(Bank), target :: tmpb(2) ! temporary Bank
integer :: hdf5_err
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates
! Initialize FORTRAN interface.
call h5open_f(hdf5_err)
! Create compound type for xyz and uvw
call h5tarray_create_f(H5T_NATIVE_DOUBLE, 1, dims, coordinates_t, hdf5_err)
! Create the compound datatype for Bank
call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(2))), hdf5_bank_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "wgt", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%wgt)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "xyz", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%xyz)), coordinates_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "uvw", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%uvw)), coordinates_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "E", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%E)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "delayed_group", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%delayed_group)), H5T_NATIVE_INTEGER, hdf5_err)
! Determine type for integer(8)
hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND)
end subroutine hdf5_initialize
!===============================================================================
! READ_COMMAND_LINE reads all parameters from the command line
!===============================================================================
subroutine read_command_line()
! Arguments were already read on C++ side (initialize.cpp). Here we just
! convert the C-style strings to Fortran style
integer :: i ! loop index
integer :: argc ! number of command line arguments
integer :: last_flag ! index of last flag
character(MAX_WORD_LEN) :: filetype
integer(HID_T) :: file_id
character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments
character(kind=C_CHAR), pointer :: string(:)
interface
function is_null(ptr) result(x) bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: ptr
logical(C_BOOL) :: x
end function is_null
end interface
! Check number of command line arguments and allocate argv
argc = COMMAND_ARGUMENT_COUNT()
! Allocate and retrieve command arguments
allocate(argv(argc))
do i = 1, argc
call GET_COMMAND_ARGUMENT(i, argv(i))
end do
! Process command arguments
last_flag = 0
i = 1
do while (i <= argc)
! Check for flags
if (starts_with(argv(i), "-")) then
select case (argv(i))
case ('-p', '-plot', '--plot')
run_mode = MODE_PLOTTING
check_overlaps = .true.
case ('-n', '-particles', '--particles')
! Read number of particles per cycle
i = i + 1
n_particles = str_to_int(argv(i))
! Check that number specified was valid
if (n_particles == ERROR_INT) then
call fatal_error("Must specify integer after " // trim(argv(i-1)) &
&// " command-line flag.")
end if
case ('-r', '-restart', '--restart')
! Read path for state point/particle restart
i = i + 1
! Check what type of file this is
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_attribute(filetype, file_id, 'filetype')
call file_close(file_id)
! Set path and flag for type of run
select case (trim(filetype))
case ('statepoint')
path_state_point = argv(i)
restart_run = .true.
case ('particle restart')
path_particle_restart = argv(i)
particle_restart_run = .true.
case default
call fatal_error("Unrecognized file after restart flag: " // filetype // ".")
end select
! If its a restart run check for additional source file
if (restart_run .and. i + 1 <= argc) then
! Increment arg
i = i + 1
! Check if it has extension we can read
if (ends_with(argv(i), '.h5')) then
! Check file type is a source file
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_attribute(filetype, file_id, 'filetype')
call file_close(file_id)
if (filetype /= 'source') then
call fatal_error("Second file after restart flag must be a &
&source file")
end if
! It is a source file
path_source_point = argv(i)
else ! Different option is specified not a source file
! Source is in statepoint file
path_source_point = path_state_point
! Set argument back
i = i - 1
end if
else ! No command line arg after statepoint
! Source is assumed to be in statepoint file
path_source_point = path_state_point
end if
case ('-g', '-geometry-debug', '--geometry-debug')
check_overlaps = .true.
case ('-c', '--volume')
run_mode = MODE_VOLUME
case ('-s', '--threads')
! Read number of threads
i = i + 1
#ifdef _OPENMP
! Read and set number of OpenMP threads
n_threads = int(str_to_int(argv(i)), 4)
if (n_threads < 1) then
call fatal_error("Invalid number of threads specified on command &
&line.")
end if
call omp_set_num_threads(n_threads)
#else
if (master) call warning("Ignoring number of threads specified on &
&command line.")
#endif
case ('-?', '-h', '-help', '--help')
call print_usage()
stop
case ('-v', '-version', '--version')
call print_version()
stop
case ('-t', '-track', '--track')
write_all_tracks = .true.
case default
call fatal_error("Unknown command line option: " // argv(i))
end select
last_flag = i
end if
! Increment counter
i = i + 1
end do
! Determine directory where XML input files are
if (argc > 0 .and. last_flag < argc) then
path_input = argv(last_flag + 1)
if (.not. is_null(openmc_path_input)) then
call c_f_pointer(openmc_path_input, string, [255])
path_input = to_f_string(string)
else
path_input = ''
end if
! Add slash at end of directory if it isn't there
if (.not. ends_with(path_input, "/") .and. len_trim(path_input) > 0) then
path_input = trim(path_input) // "/"
if (.not. is_null(openmc_path_statepoint)) then
call c_f_pointer(openmc_path_statepoint, string, [255])
path_state_point = to_f_string(string)
end if
if (.not. is_null(openmc_path_sourcepoint)) then
call c_f_pointer(openmc_path_sourcepoint, string, [255])
path_source_point = to_f_string(string)
end if
if (.not. is_null(openmc_path_particle_restart)) then
call c_f_pointer(openmc_path_particle_restart, string, [255])
path_particle_restart = to_f_string(string)
end if
! Free memory from argv
deallocate(argv)
! TODO: Check that directory exists
! Add slash at end of directory if it isn't there
if (len_trim(path_input) > 0 .and. .not. ends_with(path_input, "/")) then
path_input = trim(path_input) // "/"
end if
end subroutine read_command_line
end module initialize

229
src/initialize.cpp Normal file
View file

@ -0,0 +1,229 @@
#include "initialize.h"
#include <cstddef>
#include <cstring>
#include <iostream>
#include <sstream>
#include <string>
#include "error.h"
#include "hdf5_interface.h"
#include "message_passing.h"
#include "openmc.h"
#ifdef _OPENMP
#include "omp.h"
#endif
// data/functions from Fortran side
extern "C" bool openmc_check_overlaps;
extern "C" bool openmc_write_all_tracks;
extern "C" bool openmc_particle_restart_run;
extern "C" bool openmc_restart_run;
extern "C" void print_usage();
extern "C" void print_version();
// Paths to various files
extern "C" {
char* openmc_path_input;
char* openmc_path_statepoint;
char* openmc_path_sourcepoint;
char* openmc_path_particle_restart;
bool is_null(void* ptr) {return !ptr;}
}
int openmc_init(int argc, char* argv[], const void* intracomm)
{
#ifdef OPENMC_MPI
// Check if intracomm was passed
MPI_Comm comm;
if (intracomm) {
comm = *static_cast<const MPI_Comm *>(intracomm);
} else {
comm = MPI_COMM_WORLD;
}
// Initialize MPI for C++
openmc::initialize_mpi(comm);
#endif
// Parse command-line arguments
int err = openmc::parse_command_line(argc, argv);
if (err) return err;
// Continue with rest of initialization
#ifdef OPENMC_MPI
MPI_Fint fcomm = MPI_Comm_c2f(comm);
openmc_init_f(&fcomm);
#else
openmc_init_f(nullptr);
#endif
return 0;
}
namespace openmc {
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm)
{
openmc::mpi::intracomm = intracomm;
// Initialize MPI
int flag;
MPI_Initialized(&flag);
if (!flag) MPI_Init(nullptr, nullptr);
// Determine number of processes and rank for each
MPI_Comm_size(intracomm, &openmc::mpi::n_procs);
MPI_Comm_rank(intracomm, &openmc::mpi::rank);
// Set variable for Fortran side
openmc_n_procs = openmc::mpi::n_procs;
openmc_rank = openmc::mpi::rank;
openmc_master = (openmc::mpi::rank == 0);
// Create bank datatype
Bank b;
MPI_Aint disp[] {
offsetof(Bank, wgt),
offsetof(Bank, xyz),
offsetof(Bank, uvw),
offsetof(Bank, E),
offsetof(Bank, delayed_group)
};
int blocks[] {1, 3, 3, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
MPI_Type_create_struct(5, blocks, disp, types, &openmc::mpi::bank);
MPI_Type_commit(&openmc::mpi::bank);
}
#endif // OPENMC_MPI
inline bool ends_with(std::string const& value, std::string const& ending)
{
if (ending.size() > value.size()) return false;
return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
}
int
parse_command_line(int argc, char* argv[])
{
char buffer[256]; // buffer for reading attribute
int last_flag = 0;
for (int i=1; i < argc; ++i) {
std::string arg {argv[i]};
if (arg[0] == '-') {
if (arg == "-p" || arg == "--plot") {
openmc_run_mode = RUN_MODE_PLOTTING;
openmc_check_overlaps = true;
} else if (arg == "-n" || arg == "--particles") {
i += 1;
n_particles = std::stoll(argv[i]);
} else if (arg == "-r" || arg == "--restart") {
i += 1;
// Check what type of file this is
hid_t file_id = file_open(argv[i], 'r', true);
size_t len = attribute_typesize(file_id, "filetype");
read_attr_string(file_id, "filetype", len, buffer);
file_close(file_id);
std::string filetype {buffer};
// Set path and flag for type of run
if (filetype == "statepoint") {
openmc_path_statepoint = argv[i];
openmc_restart_run = true;
} else if (filetype == "particle restart") {
openmc_path_particle_restart = argv[i];
openmc_particle_restart_run = true;
} else {
std::stringstream msg;
msg << "Unrecognized file after restart flag: " << filetype << ".";
strcpy(openmc_err_msg, msg.str().c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
// If its a restart run check for additional source file
if (openmc_restart_run && i + 1 < argc) {
// Check if it has extension we can read
if (ends_with(argv[i+1], ".h5")) {
// Check file type is a source file
file_id = file_open(argv[i+1], 'r', true);
len = attribute_typesize(file_id, "filetype");
read_attr_string(file_id, "filetype", len, buffer);
file_close(file_id);
if (filetype != "source") {
std::string msg {"Second file after restart flag must be a source file"};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
// It is a source file
openmc_path_sourcepoint = argv[i+1];
i += 1;
} else {
// Source is in statepoint file
openmc_path_sourcepoint = openmc_path_statepoint;
}
} else {
// Source is assumed to be in statepoint file
openmc_path_sourcepoint = openmc_path_statepoint;
}
} else if (arg == "-g" || arg == "--geometry-debug") {
openmc_check_overlaps = true;
} else if (arg == "-c" || arg == "--volume") {
openmc_run_mode = RUN_MODE_VOLUME;
} else if (arg == "-s" || arg == "--threads") {
// Read number of threads
i += 1;
#ifdef _OPENMP
// Read and set number of OpenMP threads
openmc_n_threads = std::stoi(argv[i]);
if (openmc_n_threads < 1) {
std::string msg {"Number of threads must be positive."};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
omp_set_num_threads(openmc_n_threads);
#else
if (openmc_master)
warning("Ignoring number of threads specified on command line.");
#endif
} else if (arg == "-?" || arg == "-h" || arg == "--help") {
print_usage();
return OPENMC_E_UNASSIGNED;
} else if (arg == "-v" || arg == "--version") {
print_version();
return OPENMC_E_UNASSIGNED;
} else if (arg == "-t" || arg == "--track") {
openmc_write_all_tracks = true;
} else {
std::cerr << "Unknown option: " << argv[i] << '\n';
print_usage();
return OPENMC_E_UNASSIGNED;
}
last_flag = i;
}
}
// Determine directory where XML input files are
if (argc > 1 && last_flag < argc) openmc_path_input = argv[last_flag + 1];
return 0;
}
} // namespace openmc

20
src/initialize.h Normal file
View file

@ -0,0 +1,20 @@
#ifndef INITIALIZE_H
#define INITIALIZE_H
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
extern "C" void print_usage();
extern "C" void print_version();
namespace openmc {
int parse_command_line(int argc, char* argv[]);
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm);
#endif
}
#endif // INITIALIZE_H

51
src/main.F90
View file

@ -1,51 +0,0 @@
program main
use, intrinsic :: ISO_C_BINDING
use constants
use message_passing
use openmc_api, only: openmc_init, openmc_finalize, openmc_run, &
openmc_plot_geometry, openmc_calculate_volumes
use particle_restart, only: run_particle_restart
use settings, only: run_mode
implicit none
integer(C_INT) :: err
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
#endif
! Initialize run -- when run with MPI, pass communicator
#ifdef OPENMC_MPI
#ifdef OPENMC_MPIF08
call openmc_init(MPI_COMM_WORLD % MPI_VAL)
#else
call openmc_init(MPI_COMM_WORLD)
#endif
#else
call openmc_init()
#endif
! start problem based on mode
select case (run_mode)
case (MODE_FIXEDSOURCE, MODE_EIGENVALUE)
err = openmc_run()
case (MODE_PLOTTING)
call openmc_plot_geometry()
case (MODE_PARTICLE)
if (master) call run_particle_restart()
case (MODE_VOLUME)
call openmc_calculate_volumes()
end select
! finalize run
call openmc_finalize()
#ifdef OPENMC_MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
end program main

51
src/main.cpp Normal file
View file

@ -0,0 +1,51 @@
#include "error.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "openmc.h"
int main(int argc, char* argv[]) {
int err;
// Initialize run -- when run with MPI, pass communicator
#ifdef OPENMC_MPI
MPI_Comm world {MPI_COMM_WORLD};
err = openmc_init(argc, argv, &world);
#else
err = openmc_init(argc, argv, nullptr);
#endif
if (err == -1) {
// This happens for the -h and -v flags
return 0;
} else if (err) {
openmc::fatal_error(openmc_err_msg);
}
// start problem based on mode
switch (openmc_run_mode) {
case RUN_MODE_FIXEDSOURCE:
case RUN_MODE_EIGENVALUE:
err = openmc_run();
break;
case RUN_MODE_PLOTTING:
err = openmc_plot_geometry();
break;
case RUN_MODE_PARTICLE:
if (openmc_master) err = openmc_particle_restart();
break;
case RUN_MODE_VOLUME:
err = openmc_calculate_volumes();
break;
}
if (err) openmc::fatal_error(openmc_err_msg);
// Finalize and free up memory
err = openmc_finalize();
if (err) openmc::fatal_error(openmc_err_msg);
// If MPI is in use and enabled, terminate it
#ifdef OPENMC_MPI
MPI_Finalize();
#endif
}

1
src/math.F90
View file

@ -216,7 +216,6 @@ contains
real(C_DOUBLE) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is
! easier to work with
real(C_DOUBLE) :: k1, k2, k3, k4 ! Variables for R_m_n calculation
real(C_DOUBLE) :: sqrt_norm ! normalization for radial moments
integer(C_INT) :: i,p,q ! Loop counters
real(C_DOUBLE), parameter :: SQRT_N_1(0:10) = [&

2
src/mesh_header.F90
View file

@ -2,8 +2,6 @@ module mesh_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use dict_header, only: DictIntInt
use error

8
src/message_passing.F90
View file

@ -1,5 +1,7 @@
module message_passing
use, intrinsic :: ISO_C_BINDING
#ifdef OPENMC_MPI
#ifdef OPENMC_MPIF08
use mpi_f08
@ -12,9 +14,9 @@ module message_passing
! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial
! run. In this case, these variables are still used at times.
integer :: n_procs = 1 ! number of processes
integer :: rank = 0 ! rank of process
logical :: master = .true. ! master process?
integer(C_INT), bind(C, name='openmc_n_procs') :: n_procs = 1 ! number of processes
integer(C_INT), bind(C, name='openmc_rank') :: rank = 0 ! rank of process
logical(C_BOOL), bind(C, name='openmc_master') :: master = .true. ! master process?
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
#ifdef OPENMC_MPIF08
type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank

15
src/message_passing.cpp Normal file
View file

@ -0,0 +1,15 @@
#include "message_passing.h"
namespace openmc {
namespace mpi {
int rank {0};
int n_procs {1};
#ifdef OPENMC_MPI
MPI_Comm intracomm;
MPI_Datatype bank;
#endif
} // namespace mpi
} // namespace openmc

22
src/message_passing.h Normal file
View file

@ -0,0 +1,22 @@
#ifndef MESSAGE_PASSING_H
#define MESSAGE_PASSING_H
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
namespace openmc {
namespace mpi {
extern int rank;
extern int n_procs;
#ifdef OPENMC_MPI
extern MPI_Datatype bank;
extern MPI_Comm intracomm;
#endif
} // namespace mpi
} // namespace openmc
#endif // MESSAGE_PASSING_H

4
src/mgxs_data.F90
View file

@ -132,6 +132,8 @@ contains
! Add name to dictionary
call already_read % add(name)
call close_group(xsdata_group)
end if
end do NUCLIDE_LOOP
end do MATERIAL_LOOP
@ -159,6 +161,8 @@ contains
end do NUCLIDE_LOOP2
end do MATERIAL_LOOP3
call file_close(file_id)
end subroutine read_mgxs
!===============================================================================

26
src/mgxs_header.F90
View file

@ -3,8 +3,6 @@ module mgxs_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: find, sort
use constants, only: MAX_WORD_LEN, ZERO, ONE, TWO, PI
use error, only: fatal_error
@ -248,7 +246,6 @@ contains
type(VectorInt), intent(out) :: temps_to_read ! Temperatures to read
integer, intent(out) :: order_dim ! Scattering data order size
integer(SIZE_T) :: name_len
integer(HID_T) :: kT_group
character(MAX_WORD_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
@ -259,8 +256,7 @@ contains
integer :: ipol, iazi
! Get name of dataset from group
name_len = len(this % name)
this % name = get_name(xs_id, name_len)
this % name = get_name(xs_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -551,6 +547,8 @@ contains
else
call fatal_error("beta must be provided as a 1D or 2D array")
end if
call close_dataset(xsdata)
else
temp_beta = ZERO
end if
@ -681,6 +679,8 @@ contains
call fatal_error("nu-fission must be provided as a 1D or 2D &
&array")
end if
call close_dataset(xsdata)
end if
! If chi-prompt provided, set chi-prompt
@ -784,6 +784,8 @@ contains
call fatal_error("chi-delayed must be provided as a 1D or 2D &
&array")
end if
call close_dataset(xsdata)
end if
! If prompt-nu-fission present, set prompt-nu-fission
@ -840,6 +842,8 @@ contains
call fatal_error("prompt-nu-fission must be provided as a 1D &
&or 2D array")
end if
call close_dataset(xsdata)
end if
! If delayed-nu-fission provided, set delayed-nu-fission. If
@ -965,6 +969,8 @@ contains
call fatal_error("delayed-nu-fission must be provided as a &
&1D, 2D, or 3D array")
end if
call close_dataset(xsdata)
end if
! Deallocate temporary beta array
@ -1350,6 +1356,8 @@ contains
else
call fatal_error("beta must be provided as a 3D or 4D array")
end if
call close_dataset(xsdata)
else
temp_beta = ZERO
end if
@ -1521,6 +1529,8 @@ contains
call fatal_error("nu-fission must be provided as a 3D or &
&4D array")
end if
call close_dataset(xsdata)
end if
! If chi-prompt provided, set chi-prompt
@ -1655,6 +1665,8 @@ contains
call fatal_error("chi-delayed must be provided as a 3D or 4D &
&array")
end if
call close_dataset(xsdata)
end if
! If prompt-nu-fission present, set prompt-nu-fission
@ -1724,6 +1736,8 @@ contains
call fatal_error("prompt-nu-fission must be provided as a 3D &
&or 4D array")
end if
call close_dataset(xsdata)
end if
! If delayed-nu-fission provided, set delayed-nu-fission. If
@ -1872,6 +1886,8 @@ contains
call fatal_error("delayed-nu-fission must be provided as a &
&3D, 4D, or 5D array")
end if
call close_dataset(xsdata)
end if
! Deallocate temporary beta array

2
src/multipole_header.F90
View file

@ -1,7 +1,5 @@
module multipole_header
use hdf5
use constants
use dict_header, only: DictIntInt
use error, only: fatal_error

8
src/nuclide_header.F90
View file

@ -3,8 +3,6 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: sort, find, binary_search
use constants
use dict_header, only: DictIntInt, DictCharInt
@ -277,7 +275,7 @@ contains
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
logical(C_BOOL), intent(in) :: master ! if this is the master proc
integer, intent(in) :: i_nuclide ! Nuclide index in nuclides
integer :: i
@ -292,7 +290,6 @@ contains
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
character(MAX_WORD_LEN) :: temp_str
@ -306,8 +303,7 @@ contains
type(VectorInt) :: index_inelastic_scatter
! Get name of nuclide from group
name_len = len(this % name)
this % name = get_name(group_id, name_len)
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))

4
src/output.F90
View file

@ -162,7 +162,7 @@ contains
! information
!===============================================================================
subroutine print_version()
subroutine print_version() bind(C)
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I2,".",I1)') &
@ -182,7 +182,7 @@ contains
! PRINT_USAGE displays information about command line usage of OpenMC
!===============================================================================
subroutine print_usage()
subroutine print_usage() bind(C)
if (master) then
write(OUTPUT_UNIT,*) 'Usage: openmc [options] [directory]'

4
src/particle_header.F90
View file

@ -1,7 +1,5 @@
module particle_header
use hdf5, only: HID_T
use bank_header, only: Bank, source_bank
use constants
use error, only: fatal_error, warning
@ -312,7 +310,7 @@ contains
!$omp critical (WriteParticleRestart)
! Create file
file_id = file_create(filename)
file_id = file_open(filename, 'w')
associate (src => source_bank(current_work))
! Write filetype and version info

15
src/particle_restart.F90
View file

@ -5,7 +5,7 @@ module particle_restart
use bank_header, only: Bank
use constants
use error, only: write_message
use hdf5_interface, only: file_open, file_close, read_dataset
use hdf5_interface, only: file_open, file_close, read_dataset, HID_T
use mgxs_header, only: energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: print_particle
@ -16,24 +16,25 @@ module particle_restart
use tally_header, only: n_tallies
use tracking, only: transport
use hdf5, only: HID_T
implicit none
private
public :: run_particle_restart
public :: openmc_particle_restart
contains
!===============================================================================
! RUN_PARTICLE_RESTART is the main routine that runs the particle restart
! OPENMC_PARTICLE_RESTART is the main routine that runs the particle restart
!===============================================================================
subroutine run_particle_restart()
function openmc_particle_restart() result(err) bind(C)
integer(C_INT) :: err
integer(8) :: particle_seed
integer :: previous_run_mode
type(Particle) :: p
err = 0
! Set verbosity high
verbosity = 10
@ -66,7 +67,7 @@ contains
deallocate(micro_xs)
end subroutine run_particle_restart
end function openmc_particle_restart
!===============================================================================
! READ_PARTICLE_RESTART reads the particle restart file

70
src/plot.F90
View file

@ -2,8 +2,6 @@ module plot
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use error, only: fatal_error, write_message
use geometry, only: find_cell, check_cell_overlap
@ -32,7 +30,8 @@ contains
! RUN_PLOT controls the logic for making one or many plots
!===============================================================================
subroutine openmc_plot_geometry() bind(C)
function openmc_plot_geometry() result(err) bind(C)
integer(C_INT) :: err
integer :: i ! loop index for plots
@ -52,7 +51,8 @@ contains
end associate
end do
end subroutine openmc_plot_geometry
err = 0
end function openmc_plot_geometry
!===============================================================================
! POSITION_RGB computes the red/green/blue values for a given plot with the
@ -340,23 +340,46 @@ contains
subroutine create_voxel(pl)
type(ObjectPlot), intent(in) :: pl
integer :: x, y, z ! voxel location indices
integer(C_INT) :: x, y, z ! voxel location indices
integer :: rgb(3) ! colors (red, green, blue) from 0-255
integer :: id ! id of cell or material
integer :: hdf5_err
integer, target :: data(pl%pixels(3),pl%pixels(2))
integer(C_INT), target :: data(pl%pixels(3),pl%pixels(2))
integer(HID_T) :: file_id
integer(HID_T) :: dspace
integer(HID_T) :: memspace
integer(HID_T) :: dset
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
integer(HSIZE_T) :: offset(3)
real(8) :: vox(3) ! x, y, and z voxel widths
real(8) :: ll(3) ! lower left starting point for each sweep direction
type(Particle) :: p
type(ProgressBar) :: progress
type(c_ptr) :: f_ptr
interface
subroutine voxel_init(file_id, dims, dspace, dset, memspace) bind(C)
import HID_T, HSIZE_T
integer(HID_T), value :: file_id
integer(HSIZE_T), intent(in) :: dims(*)
integer(HID_T), intent(out) :: dspace
integer(HID_T), intent(out) :: dset
integer(HID_T), intent(out) :: memspace
end subroutine voxel_init
subroutine voxel_write_slice(x, dspace, dset, memspace, buf) bind(C)
import C_INT, HID_T, C_PTR
integer(C_INT), value :: x
integer(HID_T), value :: dspace
integer(HID_T), value :: dset
integer(HID_T), value :: memspace
type(C_PTR), value :: buf
end subroutine voxel_write_slice
subroutine voxel_finalize(dspace, dset, memspace) bind(C)
import HID_T
integer(HID_T), value :: dspace
integer(HID_T), value :: dset
integer(HID_T), value :: memspace
end subroutine voxel_finalize
end interface
! compute voxel widths in each direction
vox = pl % width/dble(pl % pixels)
@ -371,7 +394,7 @@ contains
p % coord(1) % universe = root_universe
! Open binary plot file for writing
file_id = file_create(pl%path_plot)
file_id = file_open(pl%path_plot, 'w')
! write header info
call write_attribute(file_id, "filetype", 'voxel')
@ -390,20 +413,8 @@ contains
! Create dataset for voxel data -- note that the dimensions are reversed
! since we want the order in the file to be z, y, x
dims(:) = [pl%pixels(3), pl%pixels(2), pl%pixels(1)]
call h5screate_simple_f(3, dims, dspace, hdf5_err)
call h5dcreate_f(file_id, "data", H5T_NATIVE_INTEGER, dspace, dset, hdf5_err)
! Create another dataspace for 2D array in memory
dims_slab(1) = pl%pixels(3)
dims_slab(2) = pl%pixels(2)
dims_slab(3) = 1
call h5screate_simple_f(2, dims_slab(1:2), memspace, hdf5_err)
! Initialize offset and get pointer to data
offset(:) = 0
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims_slab, hdf5_err)
f_ptr = c_loc(data)
dims(:) = pl % pixels
call voxel_init(file_id, dims, dspace, dset, memspace)
! move to center of voxels
ll = ll + vox / TWO
@ -433,15 +444,10 @@ contains
p % coord(1) % xyz(3) = ll(3)
! Write to HDF5 dataset
offset(3) = x - 1
call h5soffset_simple_f(dspace, offset, hdf5_err)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, &
mem_space_id=memspace, file_space_id=dspace)
call voxel_write_slice(x, dspace, dset, memspace, c_loc(data))
end do
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call voxel_finalize(dspace, dset, memspace)
call file_close(file_id)
end subroutine create_voxel

42
src/plot.cpp Normal file
View file

@ -0,0 +1,42 @@
#include "plot.h"
namespace openmc {
void
voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
hid_t* memspace)
{
// Create dataspace/dataset for voxel data
*dspace = H5Screate_simple(3, dims, nullptr);
*dset = H5Dcreate(file_id, "data", H5T_NATIVE_INT, *dspace, H5P_DEFAULT,
H5P_DEFAULT, H5P_DEFAULT);
// Create dataspace for a slice of the voxel
hsize_t dims_slice[2] {dims[1], dims[2]};
*memspace = H5Screate_simple(2, dims_slice, nullptr);
// Select hyperslab in dataspace
hsize_t start[3] {0, 0, 0};
hsize_t count[3] {1, dims[1], dims[2]};
H5Sselect_hyperslab(*dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
}
void
voxel_write_slice(int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf)
{
hssize_t offset[3] {x - 1, 0, 0};
H5Soffset_simple(dspace, offset);
H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
}
void
voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
{
H5Dclose(dset);
H5Sclose(dspace);
H5Sclose(memspace);
}
} // namespace openmc

15
src/plot.h Normal file
View file

@ -0,0 +1,15 @@
#ifndef PLOT_H
#define PLOT_H
#include "hdf5.h"
namespace openmc {
extern "C" void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
hid_t* dset, hid_t* memspace);
extern "C" void voxel_write_slice(int x, hid_t dspace, hid_t dset,
hid_t memspace, void* buf);
extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
} // namespace openmc
#endif // PLOT_H

4
src/product_header.F90
View file

@ -1,13 +1,11 @@
module product_header
use hdf5, only: HID_T
use angleenergy_header, only: AngleEnergyContainer
use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, &
EMISSION_TOTAL, NEUTRON, PHOTON
use endf_header, only: Tabulated1D, Function1D, Polynomial
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, close_dataset, read_dataset
open_dataset, close_dataset, read_dataset, HID_T
use random_lcg, only: prn
use secondary_correlated, only: CorrelatedAngleEnergy
use secondary_kalbach, only: KalbachMann

6
src/reaction_header.F90
View file

@ -1,11 +1,7 @@
module reaction_header
use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5gget_info_f, h5lget_name_by_idx_f, &
H5_INDEX_NAME_F, H5_ITER_INC_F
use constants, only: MAX_WORD_LEN
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape, get_groups
use hdf5_interface
use product_header, only: ReactionProduct
use stl_vector, only: VectorInt
use string, only: to_str, starts_with

9
src/sab_header.F90
View file

@ -8,10 +8,7 @@ module sab_header
use dict_header, only: DictIntInt, DictCharInt
use distribution_univariate, only: Tabular
use error, only: warning, fatal_error
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_datasets, object_exists, &
get_name
use hdf5_interface
use random_lcg, only: prn
use secondary_correlated, only: CorrelatedAngleEnergy
use settings
@ -103,7 +100,6 @@ contains
integer :: n_energy, n_energy_out, n_mu
integer :: i_closest
integer :: n_temperature
integer(SIZE_T) :: name_len
integer(HID_T) :: T_group
integer(HID_T) :: elastic_group
integer(HID_T) :: inelastic_group
@ -123,8 +119,7 @@ contains
type(VectorInt) :: temps_to_read
! Get name of table from group
name_len = len(this % name)
this % name = get_name(group_id, name_len)
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))

5
src/secondary_correlated.F90
View file

@ -1,13 +1,10 @@
module secondary_correlated
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, ONE, HALF, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer, Tabular
use hdf5_interface, only: get_shape, read_attribute, open_dataset, &
read_dataset, close_dataset
use hdf5_interface
use random_lcg, only: prn
!===============================================================================

5
src/secondary_kalbach.F90
View file

@ -1,12 +1,9 @@
module secondary_kalbach
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, HALF, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use hdf5_interface, only: read_attribute, read_dataset, open_dataset, &
close_dataset, get_shape
use hdf5_interface
use random_lcg, only: prn
!===============================================================================

4
src/secondary_nbody.F90
View file

@ -1,10 +1,8 @@
module secondary_nbody
use hdf5, only: HID_T
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO, PI
use hdf5_interface, only: read_attribute
use hdf5_interface, only: read_attribute, HID_T
use math, only: maxwell_spectrum
use random_lcg, only: prn

4
src/secondary_uncorrelated.F90
View file

@ -1,7 +1,5 @@
module secondary_uncorrelated
use hdf5, only: HID_T
use angle_distribution, only: AngleDistribution
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO, MAX_WORD_LEN
@ -9,7 +7,7 @@ module secondary_uncorrelated
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, DiscretePhoton
use error, only: warning
use hdf5_interface, only: read_attribute, open_group, close_group, &
object_exists
object_exists, HID_T
use random_lcg, only: prn
!===============================================================================

14
src/settings.F90
View file

@ -75,16 +75,16 @@ module settings
real(8) :: weight_survive = ONE
! Mode to run in (fixed source, eigenvalue, plotting, etc)
integer(C_INT), bind(C) :: run_mode = NONE
integer(C_INT), bind(C, name='openmc_run_mode') :: run_mode = NONE
! Restart run
logical :: restart_run = .false.
logical(C_BOOL), bind(C, name='openmc_restart_run') :: restart_run = .false.
! The verbosity controls how much information will be printed to the screen
! and in logs
integer(C_INT), bind(C) :: verbosity = 7
integer(C_INT), bind(C, name='openmc_verbosity') :: verbosity = 7
logical :: check_overlaps = .false.
logical(C_BOOL), bind(C, name='openmc_check_overlaps') :: check_overlaps = .false.
! Trace for single particle
integer :: trace_batch
@ -92,11 +92,13 @@ module settings
integer(8) :: trace_particle
! Particle tracks
logical :: write_all_tracks = .false.
logical(C_BOOL), bind(C, name='openmc_write_all_tracks') :: &
write_all_tracks = .false.
integer, allocatable :: track_identifiers(:,:)
! Particle restart run
logical :: particle_restart_run = .false.
logical(C_BOOL), bind(C, name='openmc_particle_restart_run') :: &
particle_restart_run = .false.
! Write out initial source
logical :: write_initial_source = .false.

38
src/simulation.F90
View file

@ -44,7 +44,6 @@ module simulation
implicit none
private
public :: openmc_next_batch
public :: openmc_run
public :: openmc_simulation_init
public :: openmc_simulation_finalize
@ -54,25 +53,6 @@ module simulation
contains
!===============================================================================
! 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.
!===============================================================================
function openmc_run() result(err) bind(C)
integer(C_INT) :: err
integer(C_INT) :: status
call openmc_simulation_init()
do
err = openmc_next_batch(status)
if (status /= 0 .or. err < 0) exit
end do
call openmc_simulation_finalize()
end function openmc_run
!===============================================================================
! OPENMC_NEXT_BATCH
!===============================================================================
@ -332,6 +312,7 @@ contains
subroutine finalize_batch()
integer(C_INT) :: err
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
#endif
@ -367,7 +348,7 @@ contains
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call openmc_statepoint_write()
err = openmc_statepoint_write()
end if
! Write out source point if it's been specified for this batch
@ -416,9 +397,13 @@ contains
! INITIALIZE_SIMULATION
!===============================================================================
subroutine openmc_simulation_init() bind(C)
function openmc_simulation_init() result(err) bind(C)
integer(C_INT) :: err
integer :: i
err = 0
! Skip if simulation has already been initialized
if (simulation_initialized) return
@ -476,14 +461,15 @@ contains
! Set flag indicating initialization is done
simulation_initialized = .true.
end subroutine openmc_simulation_init
end function openmc_simulation_init
!===============================================================================
! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
! execution time and results
!===============================================================================
subroutine openmc_simulation_finalize() bind(C)
function openmc_simulation_finalize() result(err) bind(C)
integer(C_INT) :: err
integer :: i ! loop index
#ifdef OPENMC_MPI
@ -499,6 +485,8 @@ contains
#endif
#endif
err = 0
! Skip if simulation was never run
if (.not. simulation_initialized) return
@ -573,7 +561,7 @@ contains
need_depletion_rx = .false.
simulation_initialized = .false.
end subroutine openmc_simulation_finalize
end function openmc_simulation_finalize
!===============================================================================
! CALCULATE_WORK determines how many particles each processor should simulate

18
src/simulation.cpp Normal file
View file

@ -0,0 +1,18 @@
#include "openmc.h"
// 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.
int openmc_run() {
openmc_simulation_init();
int err = 0;
int status = 0;
while (status == 0 && err == 0) {
err = openmc_next_batch(&status);
}
openmc_simulation_finalize();
return err;
}

13
src/simulation_header.F90
View file

@ -23,7 +23,8 @@ module simulation_header
integer :: current_batch ! current batch
integer :: current_gen ! current generation within a batch
integer :: total_gen = 0 ! total number of generations simulated
logical(C_BOOL), bind(C) :: simulation_initialized = .false.
logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: &
simulation_initialized = .false.
logical :: need_depletion_rx ! need to calculate depletion reaction rx?
! ============================================================================
@ -31,8 +32,8 @@ module simulation_header
logical :: satisfy_triggers = .false. ! whether triggers are satisfied
integer(8) :: work ! number of particles per processor
integer(8), allocatable :: work_index(:) ! starting index in source bank for each process
integer(C_INT64_T), bind(C, name='openmc_work') :: work ! number of particles per processor
integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process
integer(8) :: current_work ! index in source bank of current history simulated
! ============================================================================
@ -40,8 +41,8 @@ module simulation_header
! Temporary k-effective values
type(VectorReal) :: k_generation ! single-generation estimates of k
real(C_DOUBLE), bind(C) :: keff = ONE ! average k over active batches
real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k
real(C_DOUBLE), bind(C, name='openmc_keff') :: keff = ONE ! average k over active batches
real(C_DOUBLE), bind(C, name='openmc_keff_std') :: keff_std ! standard deviation of average k
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
@ -57,7 +58,7 @@ module simulation_header
! PARALLEL PROCESSING VARIABLES
#ifdef _OPENMP
integer :: n_threads = NONE ! number of OpenMP threads
integer(C_INT), bind(C, name='openmc_n_threads') :: n_threads = NONE ! number of OpenMP threads
integer :: thread_id ! ID of a given thread
#endif

13
src/source.F90
View file

@ -1,6 +1,5 @@
module source
use hdf5, only: HID_T
#ifdef OPENMC_MPI
use message_passing
#endif
@ -12,7 +11,7 @@ module source
use distribution_multivariate, only: SpatialBox
use error, only: fatal_error
use geometry, only: find_cell
use hdf5_interface, only: file_create, file_open, file_close, read_dataset
use hdf5_interface
use math
use message_passing, only: rank
use mgxs_header, only: rev_energy_bins, num_energy_groups
@ -55,16 +54,16 @@ contains
file_id = file_open(path_source, 'r', parallel=.true.)
! Read the file type
call read_dataset(filetype, file_id, "filetype")
call read_attribute(filetype, file_id, "filetype")
! Check to make sure this is a source file
if (filetype /= 'source') then
if (filetype /= 'source' .and. filetype /= 'statepoint') then
call fatal_error("Specified starting source file not a source file &
&type.")
end if
! Read in the source bank
call read_source_bank(file_id)
call read_source_bank(file_id, work_index, source_bank)
! Close file
call file_close(file_id)
@ -88,8 +87,8 @@ contains
if (write_initial_source) then
call write_message('Writing out initial source...', 5)
filename = trim(path_output) // 'initial_source.h5'
file_id = file_create(filename, parallel=.true.)
call write_source_bank(file_id)
file_id = file_open(filename, 'w', parallel=.true.)
call write_source_bank(file_id, work_index, source_bank)
call file_close(file_id)
end if

231
src/state_point.F90
View file

@ -13,8 +13,7 @@ module state_point
use, intrinsic :: ISO_C_BINDING
use hdf5
use bank_header, only: Bank
use cmfd_header
use constants
use eigenvalue, only: openmc_get_keff
@ -37,14 +36,31 @@ module state_point
implicit none
interface
subroutine write_source_bank(group_id, work_index, bank_) bind(C)
import HID_T, C_INT64_T, Bank
integer(HID_T), value :: group_id
integer(C_INT64_T), intent(in) :: work_index(*)
type(Bank), intent(in) :: bank_(*)
end subroutine write_source_bank
subroutine read_source_bank(group_id, work_index, bank_) bind(C)
import HID_T, C_INT64_T, Bank
integer(HID_T), value :: group_id
integer(C_INT64_T), intent(in) :: work_index(*)
type(Bank), intent(out) :: bank_(*)
end subroutine read_source_bank
end interface
contains
!===============================================================================
! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk
!===============================================================================
subroutine openmc_statepoint_write(filename) bind(C)
function openmc_statepoint_write(filename) result(err) bind(C)
type(C_PTR), intent(in), optional :: filename
integer(C_INT) :: err
integer :: i, j, k
integer :: i_xs
@ -55,12 +71,12 @@ contains
integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_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(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: filename_
err = 0
if (present(filename)) then
call c_f_pointer(filename, string, [MAX_FILE_LEN])
filename_ = to_f_string(string)
@ -76,7 +92,7 @@ contains
if (master) then
! Create statepoint file
file_id = file_create(filename_)
file_id = file_open(filename_, 'w')
! Write file type
call write_attribute(file_id, "filetype", "statepoint")
@ -443,7 +459,7 @@ contains
call file_close(file_id)
end if
end subroutine openmc_statepoint_write
end function openmc_statepoint_write
!===============================================================================
! WRITE_SOURCE_POINT
@ -476,8 +492,8 @@ contains
! Create separate source file
if (master .or. parallel) then
file_id = file_create(filename, parallel=.true.)
call write_dataset(file_id, "filetype", 'source')
file_id = file_open(filename, 'w', parallel=.true.)
call write_attribute(file_id, "filetype", 'source')
end if
else
filename = trim(path_output) // 'statepoint.' // &
@ -485,11 +501,11 @@ contains
filename = trim(filename) // '.h5'
if (master .or. parallel) then
file_id = file_open(filename, 'w', parallel=.true.)
file_id = file_open(filename, 'a', parallel=.true.)
end if
end if
call write_source_bank(file_id)
call write_source_bank(file_id, work_index, source_bank)
if (master .or. parallel) call file_close(file_id)
end if
@ -498,11 +514,11 @@ contains
filename = trim(path_output) // 'source' // '.h5'
call write_message("Creating source file " // trim(filename) // "...", 5)
if (master .or. parallel) then
file_id = file_create(filename, parallel=.true.)
call write_dataset(file_id, "filetype", 'source')
file_id = file_open(filename, 'w', parallel=.true.)
call write_attribute(file_id, "filetype", 'source')
end if
call write_source_bank(file_id)
call write_source_bank(file_id, work_index, source_bank)
if (master .or. parallel) call file_close(file_id)
end if
@ -815,14 +831,10 @@ contains
! Open source file
file_id = file_open(path_source_point, 'r', parallel=.true.)
! Read file type
call read_dataset(int_array(1), file_id, "filetype")
end if
! Write out source
call read_source_bank(file_id)
call read_source_bank(file_id, work_index, source_bank)
end if
@ -831,187 +843,4 @@ contains
end subroutine load_state_point
!===============================================================================
! WRITE_SOURCE_BANK writes OpenMC source_bank data
!===============================================================================
subroutine write_source_bank(group_id)
use bank_header, only: Bank
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: memspace ! memory space handle
integer(HSIZE_T) :: offset(1) ! source data offset
integer(HSIZE_T) :: dims(1)
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#else
integer :: i
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
type(Bank), allocatable, target :: temp_source(:)
#endif
#endif
#ifdef PHDF5
! Set size of total dataspace for all procs and rank
dims(1) = n_particles
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err)
! Create another data space but for each proc individually
dims(1) = work
call h5screate_simple_f(1, dims, memspace, hdf5_err)
! Select hyperslab for this dataspace
offset(1) = work_index(rank)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up the property list for parallel writing
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
! Write data to file in parallel
call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace, &
xfer_prp=plist)
! Close all ids
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5pclose_f(plist, hdf5_err)
#else
if (master) then
! Create dataset big enough to hold all source sites
dims(1) = n_particles
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, &
dspace, dset, hdf5_err)
! Save source bank sites since the souce_bank array is overwritten below
#ifdef OPENMC_MPI
allocate(temp_source(work))
temp_source(:) = source_bank(:)
#endif
do i = 0, n_procs - 1
! Create memory space
dims(1) = work_index(i+1) - work_index(i)
call h5screate_simple_f(1, dims, memspace, hdf5_err)
#ifdef OPENMC_MPI
! Receive source sites from other processes
if (i > 0) then
call MPI_RECV(source_bank, int(dims(1)), MPI_BANK, i, i, &
mpi_intracomm, MPI_STATUS_IGNORE, mpi_err)
end if
#endif
! Select hyperslab for this dataspace
call h5dget_space_f(dset, dspace, hdf5_err)
offset(1) = work_index(i)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up pointer to data and write data to hyperslab
f_ptr = c_loc(source_bank)
call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace)
call h5sclose_f(memspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
end do
! Close all ids
call h5dclose_f(dset, hdf5_err)
! Restore state of source bank
#ifdef OPENMC_MPI
source_bank(:) = temp_source(:)
deallocate(temp_source)
#endif
else
#ifdef OPENMC_MPI
call MPI_SEND(source_bank, int(work), MPI_BANK, 0, rank, &
mpi_intracomm, mpi_err)
#endif
end if
#endif
end subroutine write_source_bank
!===============================================================================
! READ_SOURCE_BANK reads OpenMC source_bank data
!===============================================================================
subroutine read_source_bank(group_id)
use bank_header, only: Bank
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data space handle
integer(HID_T) :: memspace ! memory space handle
integer(HSIZE_T) :: dims(1) ! dimensions on one processor
integer(HSIZE_T) :: dims_all(1) ! overall dimensions
integer(HSIZE_T) :: maxdims(1) ! maximum dimensions
integer(HSIZE_T) :: offset(1) ! offset of data
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
! Open the dataset
call h5dopen_f(group_id, "source_bank", dset, hdf5_err)
! Create another data space but for each proc individually
dims(1) = work
call h5screate_simple_f(1, dims, memspace, hdf5_err)
! Make sure source bank is big enough
call h5dget_space_f(dset, dspace, hdf5_err)
call h5sget_simple_extent_dims_f(dspace, dims_all, maxdims, hdf5_err)
if (size(source_bank, KIND=HSIZE_T) > dims_all(1)) then
call fatal_error("Number of source sites in source file is less than &
&number of source particles per generation.")
end if
! Select hyperslab for each process
offset(1) = work_index(rank)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
#ifdef PHDF5
! Read data in parallel
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err)
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace, &
xfer_prp=plist)
call h5pclose_f(plist, hdf5_err)
#else
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace)
#endif
! Close all ids
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5dclose_f(dset, hdf5_err)
end subroutine read_source_bank
end module state_point

165
src/state_point.cpp Normal file
View file

@ -0,0 +1,165 @@
#include "state_point.h"
#include <algorithm>
#include <vector>
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "error.h"
#include "message_passing.h"
#include "openmc.h"
namespace openmc {
hid_t h5banktype() {
// Create type for array of 3 reals
hsize_t dims[] {3};
hid_t triplet = H5Tarray_create(H5T_NATIVE_DOUBLE, 1, dims);
// Create bank datatype
hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct Bank));
H5Tinsert(banktype, "wgt", HOFFSET(Bank, wgt), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "xyz", HOFFSET(Bank, xyz), triplet);
H5Tinsert(banktype, "uvw", HOFFSET(Bank, uvw), triplet);
H5Tinsert(banktype, "E", HOFFSET(Bank, E), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "delayed_group", HOFFSET(Bank, delayed_group), H5T_NATIVE_INT);
H5Tclose(triplet);
return banktype;
}
void
write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
{
hid_t banktype = h5banktype();
#ifdef PHDF5
// Set size of total dataspace for all procs and rank
hsize_t dims[] {static_cast<hsize_t>(n_particles)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
// Create another data space but for each proc individually
hsize_t count[] {static_cast<hsize_t>(openmc_work)};
hid_t memspace = H5Screate_simple(1, count, nullptr);
// Select hyperslab for this dataspace
hsize_t start[] {static_cast<hsize_t>(work_index[openmc::mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Set up the property list for parallel writing
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
// Write data to file in parallel
H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank);
// Free resources
H5Sclose(dspace);
H5Sclose(memspace);
H5Dclose(dset);
H5Pclose(plist);
#else
if (openmc_master) {
// Create dataset big enough to hold all source sites
hsize_t dims[] {static_cast<hsize_t>(n_particles)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
// Save source bank sites since the souce_bank array is overwritten below
#ifdef OPENMC_MPI
std::vector<Bank> temp_source {source_bank, source_bank + openmc_work};
#endif
for (int i = 0; i < openmc::mpi::n_procs; ++i) {
// Create memory space
hsize_t count[] {static_cast<hsize_t>(work_index[i+1] - work_index[i])};
hid_t memspace = H5Screate_simple(1, count, nullptr);
#ifdef OPENMC_MPI
// Receive source sites from other processes
if (i > 0)
MPI_Recv(source_bank, count[0], openmc::mpi::bank, i, i,
openmc::mpi::intracomm, MPI_STATUS_IGNORE);
#endif
// Select hyperslab for this dataspace
dspace = H5Dget_space(dset);
hsize_t start[] {static_cast<hsize_t>(work_index[i])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Write data to hyperslab
H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
H5Sclose(memspace);
H5Sclose(dspace);
}
// Close all ids
H5Dclose(dset);
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank);
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank, openmc_work, openmc::mpi::bank, 0, openmc::mpi::rank,
openmc::mpi::intracomm);
#endif
}
#endif
H5Tclose(banktype);
}
void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
{
hid_t banktype = h5banktype();
// Open the dataset
hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT);
// Create another data space but for each proc individually
hsize_t dims[] {static_cast<hsize_t>(openmc_work)};
hid_t memspace = H5Screate_simple(1, dims, nullptr);
// Make sure source bank is big enough
hid_t dspace = H5Dget_space(dset);
hsize_t dims_all[1];
H5Sget_simple_extent_dims(dspace, dims_all, nullptr);
if (work_index[openmc::mpi::n_procs] > dims_all[0]) {
fatal_error("Number of source sites in source file is less "
"than number of source particles per generation.");
}
// Select hyperslab for each process
hsize_t start[] {static_cast<hsize_t>(work_index[openmc::mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, dims, nullptr);
#ifdef PHDF5
// Read data in parallel
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
H5Dread(dset, banktype, memspace, dspace, plist, source_bank);
H5Pclose(plist);
#else
H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
#endif
// Close all ids
H5Sclose(dspace);
H5Sclose(memspace);
H5Dclose(dset);
H5Tclose(banktype);
}
} // namespace openmc

17
src/state_point.h Normal file
View file

@ -0,0 +1,17 @@
#ifndef STATE_POINT_H
#define STATE_POINT_H
#include <cstdint>
#include "hdf5.h"
#include "openmc.h"
namespace openmc {
extern "C" void write_source_bank(hid_t group_id, int64_t* work_index,
Bank* source_bank);
extern "C" void read_source_bank(hid_t group_id, int64_t* work_index,
Bank* source_bank);
} // namespace openmc
#endif // STATE_POINT_H

4
src/summary.F90
View file

@ -1,7 +1,5 @@
module summary
use hdf5
use constants
use endf, only: reaction_name
use error, only: write_message
@ -38,7 +36,7 @@ contains
call write_message("Writing summary.h5 file...", 5)
! Create a new file using default properties.
file_id = file_create("summary.h5")
file_id = file_open("summary.h5", 'w')
call write_header(file_id)
call write_nuclides(file_id)

83
src/surface.cpp
View file

@ -11,6 +11,12 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
int32_t n_surfaces;
//==============================================================================
// Helper functions for reading the "coeffs" node of an XML surface element
//==============================================================================
@ -202,21 +208,21 @@ Surface::to_hdf5(hid_t group_id) const
switch(bc) {
case BC_TRANSMIT :
write_string(surf_group, "boundary_type", "transmission");
write_string(surf_group, "boundary_type", "transmission", false);
break;
case BC_VACUUM :
write_string(surf_group, "boundary_type", "vacuum");
write_string(surf_group, "boundary_type", "vacuum", false);
break;
case BC_REFLECT :
write_string(surf_group, "boundary_type", "reflective");
write_string(surf_group, "boundary_type", "reflective", false);
break;
case BC_PERIODIC :
write_string(surf_group, "boundary_type", "periodic");
write_string(surf_group, "boundary_type", "periodic", false);
break;
}
if (!name.empty()) {
write_string(surf_group, "name", name);
write_string(surf_group, "name", name, false);
}
to_hdf5_inner(surf_group);
@ -297,7 +303,7 @@ inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const
void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-plane");
write_string(group_id, "type", "x-plane", false);
std::array<double, 1> coeffs {{x0}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -362,7 +368,7 @@ inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const
void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-plane");
write_string(group_id, "type", "y-plane", false);
std::array<double, 1> coeffs {{y0}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -428,7 +434,7 @@ inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const
void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-plane");
write_string(group_id, "type", "z-plane", false);
std::array<double, 1> coeffs {{z0}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -489,7 +495,7 @@ SurfacePlane::normal(const double xyz[3], double uvw[3]) const
void SurfacePlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "plane");
write_string(group_id, "type", "plane", false);
std::array<double, 4> coeffs {{A, B, C, D}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -621,7 +627,7 @@ inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const
void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cylinder");
write_string(group_id, "type", "x-cylinder", false);
std::array<double, 3> coeffs {{y0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -655,7 +661,7 @@ inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const
void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-cylinder");
write_string(group_id, "type", "y-cylinder", false);
std::array<double, 3> coeffs {{x0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -689,7 +695,7 @@ inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const
void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-cylinder");
write_string(group_id, "type", "z-cylinder", false);
std::array<double, 3> coeffs {{x0, y0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -760,7 +766,7 @@ inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const
void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "sphere");
write_string(group_id, "type", "sphere", false);
std::array<double, 4> coeffs {{x0, y0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -877,7 +883,7 @@ inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const
void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cone");
write_string(group_id, "type", "x-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -911,7 +917,7 @@ inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const
void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-cone");
write_string(group_id, "type", "y-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -945,7 +951,7 @@ inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const
void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-cone");
write_string(group_id, "type", "z-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -1039,7 +1045,7 @@ SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const
void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "quadric");
write_string(group_id, "type", "quadric", false);
std::array<double, 10> coeffs {{A, B, C, D, E, F, G, H, J, K}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -1232,4 +1238,47 @@ read_surfaces(pugi::xml_node *node)
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
extern "C" int surface_id(Surface *surf) {return surf->id;}
extern "C" int surface_bc(Surface *surf) {return surf->bc;}
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3])
{return surf->sense(xyz, uvw);}
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{surf->reflect(xyz, uvw);}
extern "C" double
surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident)
{return surf->distance(xyz, uvw, coincident);}
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{return surf->normal(xyz, uvw);}
extern "C" void surface_to_hdf5(Surface *surf, hid_t group)
{surf->to_hdf5(group);}
extern "C" int surface_i_periodic(PeriodicSurface *surf)
{return surf->i_periodic;}
extern "C" bool
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
double uvw[3])
{return surf->periodic_translate(other, xyz, uvw);}
extern "C" void free_memory_surfaces_c()
{
for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
delete surfaces_c;
surfaces_c = nullptr;
n_surfaces = 0;
surface_dict.clear();
}
} // namespace openmc

56
src/surface.h
View file

@ -25,15 +25,14 @@ extern "C" const int BC_PERIODIC {3};
//==============================================================================
extern "C" double FP_COINCIDENT;
constexpr double INFTY{std::numeric_limits<double>::max()};
constexpr double INFTY {std::numeric_limits<double>::max()};
constexpr int C_NONE {-1};
//==============================================================================
// Global variables
//==============================================================================
// Braces force n_surfaces to be defined here, not just declared.
extern "C" {int32_t n_surfaces {0};}
extern "C" int32_t n_surfaces;
class Surface;
Surface **surfaces_c;
@ -384,44 +383,19 @@ public:
// Fortran compatibility functions
//==============================================================================
extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
extern "C" int surface_id(Surface *surf) {return surf->id;}
extern "C" int surface_bc(Surface *surf) {return surf->bc;}
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3])
{return surf->sense(xyz, uvw);}
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{surf->reflect(xyz, uvw);}
extern "C" double
surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident)
{return surf->distance(xyz, uvw, coincident);}
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{return surf->normal(xyz, uvw);}
extern "C" void surface_to_hdf5(Surface *surf, hid_t group)
{surf->to_hdf5(group);}
extern "C" int surface_i_periodic(PeriodicSurface *surf)
{return surf->i_periodic;}
extern "C" bool
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
double uvw[3])
{return surf->periodic_translate(other, xyz, uvw);}
extern "C" void free_memory_surfaces_c()
{
for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
delete surfaces_c;
surfaces_c = nullptr;
n_surfaces = 0;
surface_dict.clear();
}
extern "C" Surface* surface_pointer(int surf_ind);
extern "C" int surface_id(Surface *surf);
extern "C" int surface_bc(Surface *surf);
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
extern "C" double surface_distance(Surface *surf, double xyz[3], double uvw[3],
bool coincident);
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
extern "C" void surface_to_hdf5(Surface *surf, hid_t group);
extern "C" int surface_i_periodic(PeriodicSurface *surf);
extern "C" bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
double xyz[3], double uvw[3]);
extern "C" void free_memory_surfaces_c();
} // namespace openmc
#endif // SURFACE_H

5
src/surface_header.F90
View file

@ -1,9 +1,9 @@
module surface_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use dict_header, only: DictIntInt
use hdf5_interface
implicit none
@ -71,8 +71,7 @@ module surface_header
subroutine surface_to_hdf5_c(surf_ptr, group) &
bind(C, name='surface_to_hdf5')
use ISO_C_BINDING
use hdf5
import C_PTR, HID_T
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(HID_T), intent(in), value :: group

2
src/tallies/tally_filter.F90
View file

@ -2,8 +2,6 @@ module tally_filter
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use error
use string, only: to_f_string
use tally_filter_header

2
src/tallies/tally_filter_azimuthal.F90
View file

@ -2,8 +2,6 @@ module tally_filter_azimuthal
use, intrinsic :: ISO_C_BINDING
use hdf5
use algorithm, only: binary_search
use constants
use error, only: fatal_error

2
src/tallies/tally_filter_cell.F90
View file

@ -2,8 +2,6 @@ module tally_filter_cell
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error

2
src/tallies/tally_filter_cellborn.F90
View file

@ -2,8 +2,6 @@ module tally_filter_cellborn
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error

2
src/tallies/tally_filter_cellfrom.F90
View file

@ -2,8 +2,6 @@ module tally_filter_cellfrom
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error

2
src/tallies/tally_filter_delayedgroup.F90
View file

@ -2,8 +2,6 @@ module tally_filter_delayedgroup
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN, MAX_DELAYED_GROUPS
use error, only: fatal_error
use hdf5_interface

2
src/tallies/tally_filter_distribcell.F90
View file

@ -2,8 +2,6 @@ module tally_filter_distribcell
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use dict_header, only: EMPTY
use error

2
src/tallies/tally_filter_energy.F90
View file

@ -2,8 +2,6 @@ module tally_filter_energy
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use algorithm, only: binary_search
use constants
use error

2
src/tallies/tally_filter_energyfunc.F90
View file

@ -2,8 +2,6 @@ module tally_filter_energyfunc
use, intrinsic :: ISO_C_BINDING
use hdf5
use algorithm, only: binary_search
use constants
use error, only: fatal_error

3
src/tallies/tally_filter_header.F90
View file

@ -5,13 +5,12 @@ module tally_filter_header
use constants, only: MAX_LINE_LEN
use dict_header, only: DictIntInt
use error
use hdf5_interface, only: HID_T
use particle_header, only: Particle
use stl_vector, only: VectorInt, VectorReal
use string, only: to_str
use xml_interface, only: XMLNode
use hdf5
implicit none
private
public :: free_memory_tally_filter

2
src/tallies/tally_filter_legendre.F90
View file

@ -2,8 +2,6 @@ module tally_filter_legendre
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use error
use hdf5_interface

2
src/tallies/tally_filter_material.F90
View file

@ -2,8 +2,6 @@ module tally_filter_material
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use dict_header, only: DictIntInt, EMPTY
use error

2
src/tallies/tally_filter_mesh.F90
View file

@ -2,8 +2,6 @@ module tally_filter_mesh
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use dict_header, only: EMPTY
use error

2
src/tallies/tally_filter_meshsurface.F90
View file

@ -2,8 +2,6 @@ module tally_filter_meshsurface
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use dict_header, only: EMPTY
use error

2
src/tallies/tally_filter_mu.F90
View file

@ -2,8 +2,6 @@ module tally_filter_mu
use, intrinsic :: ISO_C_BINDING
use hdf5
use algorithm, only: binary_search
use constants, only: ONE, TWO, MAX_LINE_LEN, NO_BIN_FOUND
use error, only: fatal_error

2
src/tallies/tally_filter_polar.F90
View file

@ -2,8 +2,6 @@ module tally_filter_polar
use, intrinsic :: ISO_C_BINDING
use hdf5
use algorithm, only: binary_search
use constants
use error, only: fatal_error

2
src/tallies/tally_filter_sph_harm.F90
View file

@ -2,8 +2,6 @@ module tally_filter_sph_harm
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use error
use hdf5_interface

2
src/tallies/tally_filter_sptl_legendre.F90
View file

@ -2,8 +2,6 @@ module tally_filter_sptl_legendre
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use error
use hdf5_interface

2
src/tallies/tally_filter_surface.F90
View file

@ -2,8 +2,6 @@ module tally_filter_surface
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error

2
src/tallies/tally_filter_universe.F90
View file

@ -2,8 +2,6 @@ module tally_filter_universe
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error

2
src/tallies/tally_filter_zernike.F90
View file

@ -2,8 +2,6 @@ module tally_filter_zernike
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use error
use hdf5_interface

84
src/tallies/tally_header.F90
View file

@ -2,11 +2,10 @@ module tally_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use error
use dict_header, only: DictIntInt
use hdf5_interface, only: HID_T, HSIZE_T
use message_passing, only: n_procs
use nuclide_header, only: nuclide_dict
use settings, only: reduce_tallies, run_mode
@ -200,67 +199,42 @@ contains
class(TallyObject), intent(in) :: this
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset, dspace
integer(HID_T) :: memspace
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
integer(HSIZE_T) :: offset(3) = [1,0,0]
integer(HSIZE_T) :: n_filter, n_score
interface
subroutine write_tally_results(group_id, n_filter, n_score, results) &
bind(C)
import HID_T, HSIZE_T, C_DOUBLE
integer(HID_T), value :: group_id
integer(HSIZE_T), value :: n_filter
integer(HSIZE_T), value :: n_score
real(C_DOUBLE), intent(in) :: results(*)
end subroutine write_tally_results
end interface
! Create file dataspace
dims_slab(:) = shape(this % results)
dims_slab(1) = 2
call h5screate_simple_f(3, dims_slab, dspace, hdf5_err)
! Create memory dataspace that contains only SUM and SUM_SQ values
dims(:) = shape(this % results)
call h5screate_simple_f(3, dims, memspace, hdf5_err)
call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, &
hdf5_err)
! Create and write to dataset
call h5dcreate_f(group_id, "results", H5T_NATIVE_DOUBLE, dspace, dset, &
hdf5_err)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, &
hdf5_err, mem_space_id=memspace)
! Close identifiers
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
n_filter = size(this % results, 3)
n_score = size(this % results, 2)
call write_tally_results(group_id, n_filter, n_score, this % results)
end subroutine tally_write_results_hdf5
subroutine tally_read_results_hdf5(this, group_id)
class(TallyObject), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset, dspace
integer(HID_T) :: memspace
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
integer(HSIZE_T) :: offset(3) = [1,0,0]
integer(HSIZE_T) :: n_filter, n_score
interface
subroutine read_tally_results(group_id, n_filter, n_score, results) &
bind(C)
import HID_T, HSIZE_T, C_DOUBLE
integer(HID_T), value :: group_id
integer(HSIZE_T), value :: n_filter
integer(HSIZE_T), value :: n_score
real(C_DOUBLE), intent(out) :: results(*)
end subroutine read_tally_results
end interface
! Create file dataspace
dims_slab(:) = shape(this % results)
dims_slab(1) = 2
call h5screate_simple_f(3, dims_slab, dspace, hdf5_err)
! Create memory dataspace that contains only SUM and SUM_SQ values
dims(:) = shape(this % results)
call h5screate_simple_f(3, dims, memspace, hdf5_err)
call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, &
hdf5_err)
! Create and write to dataset
call h5dopen_f(group_id, "results", dset, hdf5_err)
call h5dread_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, &
hdf5_err, mem_space_id=memspace)
! Close identifiers
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
n_filter = size(this % results, 3)
n_score = size(this % results, 2)
call read_tally_results(group_id, n_filter, n_score, this % results)
end subroutine tally_read_results_hdf5
!===============================================================================

4
src/track_output.F90
View file

@ -5,8 +5,6 @@
module track_output
use hdf5
use constants
use hdf5_interface
use particle_header, only: Particle
@ -115,7 +113,7 @@ contains
end do
!$omp critical (FinalizeParticleTrack)
file_id = file_create(fname)
file_id = file_open(fname, 'w')
call write_attribute(file_id, 'filetype', 'track')
call write_attribute(file_id, 'version', VERSION_TRACK)
call write_attribute(file_id, 'n_particles', n_particle_tracks)

3
src/urr_header.F90
View file

@ -1,8 +1,7 @@
module urr_header
use hdf5, only: HID_T, HSIZE_T
use hdf5_interface, only: read_attribute, open_dataset, read_dataset, &
close_dataset, get_shape
close_dataset, get_shape, HID_T, HSIZE_T
implicit none

13
src/volume_calc.F90
View file

@ -2,7 +2,6 @@ module volume_calc
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
#ifdef _OPENMP
use omp_lib
#endif
@ -11,8 +10,8 @@ module volume_calc
use error, only: write_message
use geometry, only: find_cell
use geometry_header, only: universes, cells
use hdf5_interface, only: file_create, file_close, write_attribute, &
create_group, close_group, write_dataset
use hdf5_interface, only: file_open, file_close, write_attribute, &
create_group, close_group, write_dataset, HID_T
use output, only: header, time_stamp
use material_header, only: materials
use message_passing
@ -37,9 +36,10 @@ contains
! the user has specified and writes results to HDF5 files
!===============================================================================
subroutine openmc_calculate_volumes() bind(C)
function openmc_calculate_volumes() result(err) bind(C)
integer :: i, j
integer :: n
integer(C_INT) :: err
real(8), allocatable :: volume(:,:) ! volume mean/stdev in each domain
character(10) :: domain_type
character(MAX_FILE_LEN) :: filename ! filename for HDF5 file
@ -100,7 +100,8 @@ contains
call write_message("Elapsed time: " // trim(to_str(time_volume % &
get_value())) // " s", 6)
end if
end subroutine openmc_calculate_volumes
err = 0
end function openmc_calculate_volumes
!===============================================================================
! GET_VOLUME stochastically determines the volume of a set of domains along with
@ -435,7 +436,7 @@ contains
character(MAX_WORD_LEN), allocatable :: nucnames(:) ! names of nuclides
! Create HDF5 file
file_id = file_create(filename)
file_id = file_open(filename, 'w')
! Write header info
call write_attribute(file_id, "filetype", "volume")

40
src/xml_interface.cpp Normal file
View file

@ -0,0 +1,40 @@
#include "xml_interface.h"
#include <algorithm> // for std::transform
#include <sstream>
#include <string>
#include "pugixml/pugixml.hpp"
#include "error.h"
namespace openmc {
std::string
get_node_value(pugi::xml_node node, const char *name)
{
// Search for either an attribute or child tag and get the data as a char*.
const pugi::char_t *value_char;
if (node.attribute(name)) {
value_char = node.attribute(name).value();
} else if (node.child(name)) {
value_char = node.child_value(name);
} else {
std::stringstream err_msg;
err_msg << "Node \"" << name << "\" is not a member of the \""
<< node.name() << "\" XML node";
fatal_error(err_msg);
}
// Convert to lowercase string.
std::string value(value_char);
std::transform(value.begin(), value.end(), value.begin(), ::tolower);
// Remove whitespace.
value.erase(0, value.find_first_not_of(" \t\r\n"));
value.erase(value.find_last_not_of(" \t\r\n") + 1);
return value;
}
} // namespace openmc

32
src/xml_interface.h
View file

@ -1,8 +1,6 @@
#ifndef XML_INTERFACE_H
#define XML_INTERFACE_H
#include <algorithm> // for std::transform
#include <sstream>
#include <string>
#include "pugixml/pugixml.hpp"
@ -10,39 +8,13 @@
namespace openmc {
bool
inline bool
check_for_node(pugi::xml_node node, const char *name)
{
return node.attribute(name) || node.child(name);
}
std::string
get_node_value(pugi::xml_node node, const char *name)
{
// Search for either an attribute or child tag and get the data as a char*.
const pugi::char_t *value_char;
if (node.attribute(name)) {
value_char = node.attribute(name).value();
} else if (node.child(name)) {
value_char = node.child_value(name);
} else {
std::stringstream err_msg;
err_msg << "Node \"" << name << "\" is not a member of the \""
<< node.name() << "\" XML node";
fatal_error(err_msg);
}
// Convert to lowercase string.
std::string value(value_char);
std::transform(value.begin(), value.end(), value.begin(), ::tolower);
// Remove whitespace.
value.erase(0, value.find_first_not_of(" \t\r\n"));
value.erase(value.find_last_not_of(" \t\r\n") + 1);
return value;
}
std::string get_node_value(pugi::xml_node node, const char *name);
} // namespace openmc
#endif // XML_INTERFACE_H