diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f72d5b05..c34dd3770 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index 0e6a2536c..85509f6f4 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -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) diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index a0a62afca..cff77d2fa 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -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. diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index b0618818e..f5b277a84 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -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 diff --git a/include/openmc.h b/include/openmc.h index ed3a86054..38b7f2167 100644 --- a/include/openmc.h +++ b/include/openmc.h @@ -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 } diff --git a/openmc/capi/core.py b/openmc/capi/core.py index f80b5fe54..e044f9504 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -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) diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py index 1063d6463..d706112c4 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -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 diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index a1d8ffb0b..b40ff6344 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -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() diff --git a/openmc/statepoint.py b/openmc/statepoint.py index eb011d874..3b878e9a3 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -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): diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90 index 5d16f7424..a1f20e779 100644 --- a/src/angle_distribution.F90 +++ b/src/angle_distribution.F90 @@ -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 diff --git a/src/angleenergy_header.F90 b/src/angleenergy_header.F90 index 60d5443c4..9fda5109e 100644 --- a/src/angleenergy_header.F90 +++ b/src/angleenergy_header.F90 @@ -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 diff --git a/src/api.F90 b/src/api.F90 index 388659f4a..a06e2f4cc 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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 diff --git a/src/endf_header.F90 b/src/endf_header.F90 index e9e45ab75..9443efe2b 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -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 diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index 770da617c..a9578e5d8 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -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 diff --git a/src/error.F90 b/src/error.F90 index f4fb34173..d041051d2 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -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 diff --git a/src/error.h b/src/error.h index 4c3373b3e..45d8bcded 100644 --- a/src/error.h +++ b/src/error.h @@ -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 diff --git a/src/finalize.cpp b/src/finalize.cpp new file mode 100644 index 000000000..696a1e80e --- /dev/null +++ b/src/finalize.cpp @@ -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 +} diff --git a/src/finalize.h b/src/finalize.h new file mode 100644 index 000000000..e606493ca --- /dev/null +++ b/src/finalize.h @@ -0,0 +1,6 @@ +#ifndef FINALIZE_H +#define FINALIZE_H + +extern "C" void openmc_free_bank(); + +#endif // FINALIZE_H diff --git a/src/geometry.F90 b/src/geometry.F90 index 0bacd7b5d..e747b619a 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 410149d9e..40f56d88e 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -12,27 +12,22 @@ module hdf5_interface use, intrinsic :: ISO_C_BINDING - use hdf5 - use h5lt - use error, only: fatal_error #ifdef PHDF5 use message_passing, only: mpi_intracomm, MPI_INFO_NULL #endif + use string, only: to_c_string, to_f_string implicit none private - integer(HID_T), public :: hdf5_bank_t ! Compound type for Bank - integer(HID_T), public :: hdf5_integer8_t ! type for integer(8) - interface write_dataset - module procedure write_double + module procedure write_double_0D module procedure write_double_1D module procedure write_double_2D module procedure write_double_3D module procedure write_double_4D - module procedure write_integer + module procedure write_integer_0D module procedure write_integer_1D module procedure write_integer_2D module procedure write_integer_3D @@ -43,12 +38,12 @@ module hdf5_interface end interface write_dataset interface read_dataset - module procedure read_double + module procedure read_double_0D module procedure read_double_1D module procedure read_double_2D module procedure read_double_3D module procedure read_double_4D - module procedure read_integer + module procedure read_integer_0D module procedure read_integer_1D module procedure read_integer_2D module procedure read_integer_3D @@ -84,7 +79,6 @@ module hdf5_interface public :: attribute_exists public :: write_attribute public :: read_attribute - public :: file_create public :: file_open public :: file_close public :: create_group @@ -99,114 +93,220 @@ module hdf5_interface public :: get_datasets public :: get_name + integer, public, parameter :: HID_T = C_INT64_T + integer, public, parameter :: HSIZE_T = C_LONG_LONG + + interface + function attribute_typesize(obj_id, name) result(sz) bind(C) + import HID_T, C_CHAR, C_SIZE_T + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(C_SIZE_T) :: sz + end function attribute_typesize + + function dataset_typesize(dset) result(sz) bind(C) + import HID_T, C_SIZE_T + integer(HID_T), value :: dset + integer(C_SIZE_T) :: sz + end function dataset_typesize + + subroutine file_close(file_id) bind(C) + import HID_T + integer(HID_T), value :: file_id + end subroutine file_close + + subroutine get_shape_c(obj_id, dims) bind(C, name='get_shape') + import HID_T, HSIZE_T + integer(HID_T), value :: obj_id + integer(HSIZE_T), intent(out) :: dims(*) + end subroutine get_shape_c + + subroutine get_shape_attr(obj_id, name, dims) bind(C) + import HID_T, HSIZE_T, C_CHAR + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(HSIZE_T), intent(out) :: dims(*) + end subroutine get_shape_attr + + subroutine read_double_c(obj_id, name, buffer, indep) & + bind(C, name='read_double') + import HID_T, C_DOUBLE, C_BOOL, C_PTR + integer(HID_T), value :: obj_id + type(C_PTR), value :: name + real(C_DOUBLE), intent(out) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine read_double_c + + subroutine read_attr_int_c(obj_id, name, buffer) & + bind(C, name='read_attr_int') + import HID_T, C_CHAR, C_INT + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(C_INT), intent(out) :: buffer(*) + end subroutine read_attr_int_c + + subroutine read_attr_double_c(obj_id, name, buffer) & + bind(C, name='read_attr_double') + import HID_T, C_CHAR, C_DOUBLE + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + real(C_DOUBLE), intent(out) :: buffer(*) + end subroutine read_attr_double_c + + subroutine read_attr_string_c(obj_id, name, slen, buffer) & + bind(C, name='read_attr_string') + import HID_T, C_CHAR, C_SIZE_T + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(C_SIZE_T), value :: slen + character(kind=C_CHAR), intent(out) :: buffer(*) + end subroutine read_attr_string_c + + subroutine read_int_c(obj_id, name, buffer, indep) & + bind(C, name='read_int') + import HID_T, C_INT, C_BOOL, C_PTR + integer(HID_T), value :: obj_id + type(C_PTR), value :: name + integer(C_INT), intent(out) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine read_int_c + + subroutine read_llong_c(obj_id, name, buffer, indep) & + bind(C, name='read_llong') + import HID_T, C_INT, C_BOOL, C_PTR, C_LONG_LONG + integer(HID_T), value :: obj_id + type(C_PTR), value :: name + integer(C_LONG_LONG), intent(out) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine read_llong_c + + subroutine read_string_c(obj_id, name, slen, buffer, indep) & + bind(C, name='read_string') + import HID_T, C_PTR, C_SIZE_T, C_CHAR, C_BOOL + integer(HID_T), value :: obj_id + type(C_PTR), value :: name + integer(C_SIZE_T), value :: slen + character(kind=C_CHAR), intent(out) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine read_string_c + + subroutine read_complex_c(obj_id, name, buffer, indep) & + bind(C, name='read_complex') + import HID_T, C_PTR, C_DOUBLE_COMPLEX, C_BOOL + integer(HID_T), value :: obj_id + type(C_PTR), value :: name + complex(C_DOUBLE_COMPLEX), intent(out) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine read_complex_c + + subroutine write_attr_double_c(obj_id, ndim, dims, name, buffer) & + bind(C, name='write_attr_double') + import HID_T, HSIZE_T, C_INT, C_DOUBLE, C_CHAR + integer(HID_T), value :: obj_id + integer(C_INT), value :: ndim + integer(HSIZE_T), intent(in) :: dims(*) + character(kind=C_CHAR), intent(in) :: name(*) + real(C_DOUBLE), intent(in) :: buffer(*) + end subroutine write_attr_double_c + + subroutine write_attr_int_c(obj_id, ndim, dims, name, buffer) & + bind(C, name='write_attr_int') + import HID_T, HSIZE_T, C_INT, C_CHAR + integer(HID_T), value :: obj_id + integer(C_INT), value :: ndim + integer(HSIZE_T), intent(in) :: dims(*) + character(kind=C_CHAR), intent(in) :: name(*) + integer(C_INT), intent(in) :: buffer(*) + end subroutine write_attr_int_c + + subroutine write_attr_string_c(obj_id, name, buffer) & + bind(C, name='write_attr_string') + import HID_T, C_CHAR + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + character(kind=C_CHAR), intent(in) :: buffer(*) + end subroutine write_attr_string_c + + subroutine write_double_c(group_id, ndim, dims, name, buffer, indep) & + bind(C, name='write_double') + import HID_T, HSIZE_T, C_INT, C_DOUBLE, C_CHAR, C_BOOL + integer(HID_T), value :: group_id + integer(C_INT), value :: ndim + integer(HSIZE_T), intent(in) :: dims(*) + character(kind=C_CHAR), intent(in) :: name(*) + real(C_DOUBLE), intent(in) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine write_double_c + + subroutine write_int_c(group_id, ndim, dims, name, buffer, indep) & + bind(C, name='write_int') + import HID_T, HSIZE_T, C_INT, C_CHAR, C_BOOL + integer(HID_T), value :: group_id + integer(C_INT), value :: ndim + integer(HSIZE_T), intent(in) :: dims(*) + character(kind=C_CHAR), intent(in) :: name(*) + integer(C_INT), intent(in) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine write_int_c + + subroutine write_llong_c(group_id, ndim, dims, name, buffer, indep) & + bind(C, name='write_llong') + import HID_T, HSIZE_T, C_INT, C_CHAR, C_BOOL, C_LONG_LONG + integer(HID_T), value :: group_id + integer(C_INT), value :: ndim + integer(HSIZE_T), intent(in) :: dims(*) + character(kind=C_CHAR), intent(in) :: name(*) + integer(C_LONG_LONG), intent(in) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine write_llong_c + + subroutine write_string_c(group_id, ndim, dims, slen, name, buffer, indep) & + bind(C, name='write_string') + import HID_T, HSIZE_T, C_INT, C_CHAR, C_BOOL, C_SIZE_T + integer(HID_T), value :: group_id + integer(C_INT), value :: ndim + integer(HSIZE_T), intent(in) :: dims(*) + integer(C_SIZE_T), value :: slen + character(kind=C_CHAR), intent(in) :: name(*) + character(kind=C_CHAR), intent(in) :: buffer(*) + logical(C_BOOL), value :: indep + end subroutine write_string_c + end interface + contains -!=============================================================================== -! FILE_CREATE creates HDF5 file -!=============================================================================== - - function file_create(filename, parallel) result(file_id) - character(*), intent(in) :: filename ! name of file - logical, optional, intent(in) :: parallel ! whether to write in serial - integer(HID_T) :: file_id - - integer(HID_T) :: plist ! property list handle - integer :: hdf5_err ! HDF5 error code - logical :: parallel_ - - ! Check for serial option - parallel_ = .false. -#ifdef PHDF5 - if (present(parallel)) parallel_ = parallel -#endif - - if (parallel_) then - ! Setup file access property list with parallel I/O access - call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) -#ifdef PHDF5 -#ifdef OPENMC_MPIF08 - call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & - MPI_INFO_NULL%MPI_VAL, hdf5_err) -#else - call h5pset_fapl_mpio_f(plist, mpi_intracomm, MPI_INFO_NULL, hdf5_err) -#endif -#endif - - ! Create the file collectively - call h5fcreate_f(trim(filename), H5F_ACC_TRUNC_F, file_id, hdf5_err, & - access_prp = plist) - - ! Close the property list - call h5pclose_f(plist, hdf5_err) - else - ! Create the file - call h5fcreate_f(trim(filename), H5F_ACC_TRUNC_F, file_id, hdf5_err) - end if - - end function file_create - !=============================================================================== ! FILE_OPEN opens HDF5 file !=============================================================================== function file_open(filename, mode, parallel) result(file_id) - character(*), intent(in) :: filename ! name of file - character(*), intent(in) :: mode ! access mode to file - logical, optional, intent(in) :: parallel ! whether to write in serial + character(*), intent(in) :: filename ! name of file + character, value :: mode ! access mode to file + logical, optional, intent(in) :: parallel ! whether to write in serial integer(HID_T) :: file_id - logical :: parallel_ - integer(HID_T) :: plist ! property list handle - integer :: hdf5_err ! HDF5 error code - integer :: open_mode ! HDF5 open mode + character(kind=C_CHAR) :: mode_ + logical(C_BOOL) :: parallel_ - ! Check for serial option + interface + function file_open_c(name, mode, parallel) bind(C, name='file_open') result(file_id) + import HID_T, C_CHAR, C_BOOL, C_INT + character(kind=C_CHAR) :: name(*) + character(kind=C_CHAR), value :: mode + logical(C_BOOL), value :: parallel + integer(HID_T) :: file_id + end function file_open_c + end interface + + mode_ = mode parallel_ = .false. #ifdef PHDF5 if (present(parallel)) parallel_ = parallel #endif - ! Determine access type - open_mode = H5F_ACC_RDONLY_F - if (mode == 'w') open_mode = H5F_ACC_RDWR_F - - if (parallel_) then - ! Setup file access property list with parallel I/O access - call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) -#ifdef PHDF5 -#ifdef OPENMC_MPIF08 - call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & - MPI_INFO_NULL%MPI_VAL, hdf5_err) -#else - call h5pset_fapl_mpio_f(plist, mpi_intracomm, MPI_INFO_NULL, hdf5_err) -#endif -#endif - - ! Open the file collectively - call h5fopen_f(trim(filename), open_mode, file_id, hdf5_err, & - access_prp = plist) - - ! Close the property list - call h5pclose_f(plist, hdf5_err) - else - ! Open file - call h5fopen_f(trim(filename), open_mode, file_id, hdf5_err) - end if - + file_id = file_open_c(to_c_string(filename), mode, parallel_) end function file_open -!=============================================================================== -! FILE_CLOSE closes HDF5 file -!=============================================================================== - - subroutine file_close(file_id) - integer(HID_T), intent(in) :: file_id - - integer :: hdf5_err - - call h5fclose_f(file_id, hdf5_err) - end subroutine file_close - !=============================================================================== ! GET_GROUPS Gets a list of all the groups in a given location. !=============================================================================== @@ -215,37 +315,43 @@ contains integer(HID_T), intent(in) :: object_id character(len=150), allocatable, intent(out) :: names(:) - integer :: n_members, i, group_count, type - integer :: hdf5_err - character(len=150) :: name + integer :: i + integer(C_INT) :: n + character(len=150,kind=C_CHAR), target, allocatable :: names_(:) + type(C_PTR), allocatable :: name_ptrs(:) - ! Get number of members in this location - call h5gn_members_f(object_id, './', n_members, hdf5_err) + interface + function get_num_groups(group_id) result(n) bind(C) + import HID_T, C_INT + integer(HID_T), value :: group_id + integer(C_INT) :: n + end function get_num_groups + subroutine get_groups_c(group_id, name) bind(C, name='get_groups') + import HID_T, C_PTR + integer(HID_T), value :: group_id + type(C_PTR) :: name(*) + end subroutine get_groups_c + end interface - ! Get the number of groups - group_count = 0 - do i = 0, n_members - 1 - call h5gget_obj_info_idx_f(object_id, "./", i, name, type, hdf5_err) - if (type == H5G_GROUP_F) then - group_count = group_count + 1 - end if + ! Determine number of groups and allocate + n = get_num_groups(object_id) + allocate(names(n), names_(n), name_ptrs(n)) + + ! Set C pointers to beginning of each string + do i = 1, size(names) + name_ptrs(i) = c_loc(names_(i)) end do - ! Now we can allocate the storage for the ids - allocate(names(group_count)) - group_count = 0 - do i = 0, n_members - 1 - call h5gget_obj_info_idx_f(object_id, "./", i, name, type, hdf5_err) - if (type == H5G_GROUP_F) then - group_count = group_count + 1 - names(group_count) = trim(name) - end if + ! Get names of groups and copy to Fortran strings + call get_groups_c(object_id, name_ptrs) + do i = 1, size(names) + names(i) = to_f_string(names_(i)) end do end subroutine get_groups !=============================================================================== -! CHECK_ATTRIBUTE Checks to see if an attribute exists in the object +! ATTRIBUTE_EXISTS checks to see if an attribute exists in the object !=============================================================================== function attribute_exists(object_id, name) result(exists) @@ -253,11 +359,17 @@ contains character(*), intent(in) :: name ! name of group logical :: exists - integer :: hdf5_err ! HDF5 error code - - ! Check if attribute exists - call h5aexists_by_name_f(object_id, '.', trim(name), exists, hdf5_err) + interface + function attribute_exists_c(obj_id, name) result(exists) & + bind(C, name='attribute_exists') + import HID_T, C_CHAR, C_BOOL + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + logical(C_BOOL) :: exists + end function attribute_exists_c + end interface + exists = attribute_exists_c(object_id, to_c_string(name)) end function attribute_exists !=============================================================================== @@ -269,11 +381,18 @@ contains character(*), intent(in) :: name ! name of group logical :: exists - integer :: hdf5_err ! HDF5 error code + interface + function object_exists_c(obj_id, name) result(exists) & + bind(C, name='object_exists') + import HID_T, C_CHAR, C_BOOL + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(in) :: name(*) + logical(C_BOOL) :: exists + end function object_exists_c + end interface ! Check if group exists - call h5ltpath_valid_f(object_id, trim(name), .true., exists, hdf5_err) - + exists = object_exists_c(object_id, to_c_string(name)) end function object_exists !=============================================================================== @@ -284,32 +403,37 @@ contains integer(HID_T), intent(in) :: object_id character(len=150), allocatable, intent(out) :: names(:) - integer :: n_members, i, dset_count, type - integer :: hdf5_err - character(len=150) :: name + integer :: i + integer(C_INT) :: n + character(len=150,kind=C_CHAR), target, allocatable :: names_(:) + type(C_PTR), allocatable :: name_ptrs(:) + interface + function get_num_datasets(group_id) result(n) bind(C) + import HID_T, C_INT + integer(HID_T), value :: group_id + integer(C_INT) :: n + end function get_num_datasets + subroutine get_datasets_c(group_id, name) bind(C, name='get_datasets') + import HID_T, C_PTR + integer(HID_T), value :: group_id + type(C_PTR) :: name(*) + end subroutine get_datasets_c + end interface - ! Get number of members in this location - call h5gn_members_f(object_id, './', n_members, hdf5_err) + ! Determine number of datasets and allocate + n = get_num_datasets(object_id) + allocate(names(n), names_(n), name_ptrs(n)) - ! Get the number of datasets - dset_count = 0 - do i = 0, n_members - 1 - call h5gget_obj_info_idx_f(object_id, "./", i, name, type, hdf5_err) - if (type == H5G_DATASET_F ) then - dset_count = dset_count + 1 - end if + ! Set C pointers to beginning of each string + do i = 1, size(names) + name_ptrs(i) = c_loc(names_(i)) end do - ! Now we can allocate the storage for the ids - allocate(names(dset_count)) - dset_count = 0 - do i = 0, n_members - 1 - call h5gget_obj_info_idx_f(object_id, "./", i, name, type, hdf5_err) - if (type == H5G_DATASET_F ) then - dset_count = dset_count + 1 - names(dset_count) = trim(name) - end if + ! Get names of datasets and copy to Fortran strings + call get_datasets_c(object_id, name_ptrs) + do i = 1, size(names) + names(i) = to_f_string(names_(i)) end do end subroutine get_datasets @@ -318,21 +442,21 @@ contains ! GET_NAME Obtains the name of the current group in group_id !=============================================================================== - function get_name(group_id, name_len_) result(name) - integer(HID_T), intent(in) :: group_id - integer(SIZE_T), optional, intent(in) :: name_len_ + function get_name(object_id) result(name) + integer(HID_T), intent(in) :: object_id - character(len=150) :: name ! name of group - integer(SIZE_T) :: name_len, name_file_len - integer :: hdf5_err ! HDF5 error code + character(150) :: name ! name of object + character(kind=C_CHAR) :: name_(150) + interface + subroutine get_name_c(obj_id, name) bind(C, name='get_name') + import HID_T, C_CHAR + integer(HID_T), value :: obj_id + character(kind=C_CHAR), intent(out) :: name(*) + end subroutine get_name_c + end interface - if (present(name_len_)) then - name_len = name_len_ - else - name_len = 150 - end if - - call h5iget_name_f(group_id, name, name_len, name_file_len, hdf5_err) + call get_name_c(object_id, name_) + name = to_f_string(name_) end function get_name !=============================================================================== @@ -344,18 +468,17 @@ contains character(*), intent(in) :: name ! name of group integer(HID_T) :: newgroup_id - logical :: exists ! does the group exist - integer :: hdf5_err ! HDF5 error code + interface + function open_group_c(group_id, name) result(newgroup_id) & + bind(C, name='open_group') + import HID_T, C_CHAR + integer(HID_T), value :: group_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(HID_T) :: newgroup_id + end function open_group_c + end interface - ! Check if group exists - exists = object_exists(group_id, name) - - ! open group if it exists - if (exists) then - call h5gopen_f(group_id, trim(name), newgroup_id, hdf5_err) - else - call fatal_error("The group '" // trim(name) // "' does not exist.") - end if + newgroup_id = open_group_c(group_id, to_c_string(name)) end function open_group !=============================================================================== @@ -367,18 +490,17 @@ contains character(*), intent(in) :: name ! name of group integer(HID_T) :: newgroup_id - integer :: hdf5_err ! HDF5 error code - logical :: exists ! does the group exist + interface + function create_group_c(group_id, name) result(newgroup_id) & + bind(C, name='create_group') + import HID_T, C_CHAR + integer(HID_T), value :: group_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(HID_T) :: newgroup_id + end function create_group_c + end interface - ! Check if group exists - exists = object_exists(group_id, name) - - ! create group - if (exists) then - call fatal_error("The group '" // trim(name) // "' already exists.") - else - call h5gcreate_f(group_id, trim(name), newgroup_id, hdf5_err) - end if + newgroup_id = create_group_c(group_id, to_c_string(name)) end function create_group !=============================================================================== @@ -387,13 +509,14 @@ contains subroutine close_group(group_id) integer(HID_T), intent(inout) :: group_id + interface + subroutine close_group_c(group_id) bind(C, name='close_group') + import HID_T + integer(HID_T), value :: group_id + end subroutine close_group_c + end interface - integer :: hdf5_err ! HDF5 error code - - call h5gclose_f(group_id, hdf5_err) - if (hdf5_err < 0) then - call fatal_error("Unable to close HDF5 group.") - end if + call close_group_c(group_id) end subroutine close_group !=============================================================================== @@ -405,132 +528,56 @@ contains character(*), intent(in) :: name ! name of dataset integer(HID_T) :: dataset_id - logical :: exists ! does the dataset exist - integer :: hdf5_err ! HDF5 error code + interface + function open_dataset_c(group_id, name) result(dset_id) & + bind(C, name='open_dataset') + import HID_T, C_CHAR + integer(HID_T), value :: group_id + character(kind=C_CHAR), intent(in) :: name(*) + integer(HID_T) :: dset_id + end function open_dataset_c + end interface - ! Check if group exists - exists = object_exists(group_id, name) - - ! open group if it exists - if (exists) then - call h5dopen_f(group_id, trim(name), dataset_id, hdf5_err) - else - call fatal_error("The dataset '" // trim(name) // "' does not exist.") - end if + dataset_id = open_dataset_c(group_id, to_c_string(name)) end function open_dataset !=============================================================================== -! CLOSE_GROUP closes HDF5 temp_group +! CLOSE_DATASET closes an HDF5 dataset !=============================================================================== subroutine close_dataset(dataset_id) integer(HID_T), intent(inout) :: dataset_id - integer :: hdf5_err ! HDF5 error code + interface + subroutine close_dataset_c(dset_id) bind(C, name='close_dataset') + import HID_T + integer(HID_T), value :: dset_id + end subroutine close_dataset_c + end interface - call h5dclose_f(dataset_id, hdf5_err) - if (hdf5_err < 0) then - call fatal_error("Unable to close HDF5 dataset.") - end if + call close_dataset_c(dataset_id) end subroutine close_dataset !=============================================================================== -! WRITE_DOUBLE writes double precision scalar data +! WRITE_DOUBLE_ND writes double precision N-D array data !=============================================================================== - subroutine write_double(group_id, name, buffer, indep) + subroutine write_double_0D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name for data real(8), intent(in), target :: buffer ! data to write logical, intent(in), optional :: indep ! independent I/O - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr + integer(HSIZE_T) :: dims(0) + logical(C_BOOL) :: indep_ + real(C_DOUBLE) :: buffer_(1) - ! Set up independentive vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + indep_ = .false. + if (present(indep)) indep_ = indep + buffer_(1) = buffer - ! Create dataspace and dataset - call h5screate_f(H5S_SCALAR_F, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_DOUBLE, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_double - -!=============================================================================== -! READ_DOUBLE reads double precision scalar data -!=============================================================================== - - subroutine read_double(buffer, obj_id, name, indep) - real(8), target, intent(inout) :: buffer - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset_id - type(c_ptr) :: f_ptr - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_double - -!=============================================================================== -! WRITE_DOUBLE_1D writes double precision 1-D array data -!=============================================================================== + call write_double_c(group_id, 0, dims, to_c_string(name), buffer_, indep_) + end subroutine write_double_0D subroutine write_double_1D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id @@ -539,127 +586,15 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(1) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_double_1D_explicit(group_id, dims, name, buffer, indep) - else - call write_double_1D_explicit(group_id, dims, name, buffer) - end if + dims(1) = size(buffer) + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_double_c(group_id, 1, dims, to_c_string(name), buffer, indep_) end subroutine write_double_1D - subroutine write_double_1D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name ! name of data - real(8), intent(in), target :: buffer(dims(1)) ! data to write - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_DOUBLE, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_double_1D_explicit - -!=============================================================================== -! READ_DOUBLE_1D reads double precision 1-D array data -!=============================================================================== - - subroutine read_double_1D(buffer, obj_id, name, indep) - real(8), target, intent(inout) :: buffer(:) - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1) - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_double_1D_explicit(dset_id, dims, buffer, indep) - else - call read_double_1D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_double_1D - - subroutine read_double_1D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(1) - real(8), target, intent(inout) :: buffer(dims(1)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - end subroutine read_double_1D_explicit - -!=============================================================================== -! WRITE_DOUBLE_2D writes double precision 2-D array data -!=============================================================================== - subroutine write_double_2D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name of data @@ -667,127 +602,18 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(2) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_double_2D_explicit(group_id, dims, name, buffer, indep) - else - call write_double_2D_explicit(group_id, dims, name, buffer) - end if + ! Reverse shape of array since it will be written from C + dims(1) = size(buffer, 2) + dims(2) = size(buffer, 1) + + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_double_c(group_id, 2, dims, to_c_string(name), buffer, indep_) end subroutine write_double_2D - subroutine write_double_2D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(2) - character(*), intent(in) :: name ! name of data - real(8), intent(in), target :: buffer(dims(1),dims(2)) - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - call h5screate_simple_f(2, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_DOUBLE, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_double_2D_explicit - -!=============================================================================== -! READ_DOUBLE_2D reads double precision 2-D array data -!=============================================================================== - - subroutine read_double_2D(buffer, obj_id, name, indep) - real(8), target, intent(inout) :: buffer(:,:) - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(2) - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_double_2D_explicit(dset_id, dims, buffer, indep) - else - call read_double_2D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_double_2D - - subroutine read_double_2D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(2) - real(8), target, intent(inout) :: buffer(dims(1),dims(2)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - end subroutine read_double_2D_explicit - -!=============================================================================== -! WRITE_DOUBLE_3D writes double precision 3-D array data -!=============================================================================== - subroutine write_double_3D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name of data @@ -795,127 +621,19 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(3) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_double_3D_explicit(group_id, dims, name, buffer, indep) - else - call write_double_3D_explicit(group_id, dims, name, buffer) - end if + ! Reverse shape of array since it will be written from C + dims(1) = size(buffer, 3) + dims(2) = size(buffer, 2) + dims(3) = size(buffer, 1) + + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_double_c(group_id, 3, dims, to_c_string(name), buffer, indep_) end subroutine write_double_3D - subroutine write_double_3D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(3) - character(*), intent(in) :: name ! name of data - real(8), intent(in), target :: buffer(dims(1),dims(2),dims(3)) - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - call h5screate_simple_f(3, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_DOUBLE, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_double_3D_explicit - -!=============================================================================== -! READ_DOUBLE_3D reads double precision 3-D array data -!=============================================================================== - - subroutine read_double_3D(buffer, obj_id, name, indep) - real(8), target, intent(inout) :: buffer(:,:,:) - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(3) - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_double_3D_explicit(dset_id, dims, buffer, indep) - else - call read_double_3D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_double_3D - - subroutine read_double_3D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(3) - real(8), target, intent(inout) :: buffer(dims(1),dims(2),dims(3)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - end subroutine read_double_3D_explicit - -!=============================================================================== -! WRITE_DOUBLE_4D writes double precision 4-D array data -!=============================================================================== - subroutine write_double_4D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name of data @@ -923,60 +641,113 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(4) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_double_4D_explicit(group_id, dims, name, buffer, indep) - else - call write_double_4D_explicit(group_id, dims, name, buffer) - end if + ! Reverse shape of array since it will be written from C + dims(1) = size(buffer, 4) + dims(2) = size(buffer, 3) + dims(3) = size(buffer, 2) + dims(4) = size(buffer, 1) + + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_double_c(group_id, 4, dims, to_c_string(name), buffer, indep_) end subroutine write_double_4D - subroutine write_double_4D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(4) - character(*), intent(in) :: name ! name of data - real(8), intent(in), target :: buffer(dims(1),dims(2),dims(3),dims(4)) - logical, intent(in), optional :: indep ! independent I/O +!=============================================================================== +! READ_DOUBLE_ND reads double precision N-D array data +!=============================================================================== - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr + subroutine read_double_0D(buffer, obj_id, name, indep) + real(8), target, intent(inout) :: buffer + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + real(C_DOUBLE) :: buffer_(1) + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) - call h5screate_simple_f(4, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_DOUBLE, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) + indep_ = .false. + if (present(indep)) indep_ = indep - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_double_c(obj_id, c_loc(name_), buffer_, indep_) else - call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) + call read_double_c(obj_id, C_NULL_PTR, buffer_, indep_) end if + buffer = buffer_(1) + end subroutine read_double_0D - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_double_4D_explicit + subroutine read_double_1D(buffer, obj_id, name, indep) + real(8), target, intent(inout) :: buffer(:) + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O -!=============================================================================== -! READ_DOUBLE_4D reads double precision 4-D array data -!=============================================================================== + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep + + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_double_c(obj_id, c_loc(name_), buffer, indep_) + else + call read_double_c(obj_id, C_NULL_PTR, buffer, indep_) + end if + end subroutine read_double_1D + + subroutine read_double_2D(buffer, obj_id, name, indep) + real(8), target, intent(inout) :: buffer(:,:) + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O + + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep + + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_double_c(obj_id, c_loc(name_), buffer, indep_) + else + call read_double_c(obj_id, C_NULL_PTR, buffer, indep_) + end if + end subroutine read_double_2D + + subroutine read_double_3D(buffer, obj_id, name, indep) + real(8), target, intent(inout) :: buffer(:,:,:) + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O + + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep + + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_double_c(obj_id, c_loc(name_), buffer, indep_) + else + call read_double_c(obj_id, C_NULL_PTR, buffer, indep_) + end if + end subroutine read_double_3D subroutine read_double_4D(buffer, obj_id, name, indep) real(8), target, intent(inout) :: buffer(:,:,:,:) @@ -984,159 +755,42 @@ contains character(*), optional, intent(in) :: name logical, optional, intent(in) :: indep ! independent I/O - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(4) + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) + name_ = to_c_string(name) + call read_double_c(obj_id, c_loc(name_), buffer, indep_) else - dset_id = obj_id + call read_double_c(obj_id, C_NULL_PTR, buffer, indep_) end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_double_4D_explicit(dset_id, dims, buffer, indep) - else - call read_double_4D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) end subroutine read_double_4D - subroutine read_double_4D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(4) - real(8), target, intent(inout) :: buffer(dims(1),dims(2),dims(3),dims(4)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end if - end subroutine read_double_4D_explicit - !=============================================================================== -! WRITE_INTEGER writes integer precision scalar data +! WRITE_INTEGER_ND writes integer precision N-D array data !=============================================================================== - subroutine write_integer(group_id, name, buffer, indep) + subroutine write_integer_0D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name for data integer, intent(in), target :: buffer ! data to write logical, intent(in), optional :: indep ! independent I/O - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr + integer(HSIZE_T) :: dims(0) + logical(C_BOOL) :: indep_ + integer(C_INT) :: buffer_(1) - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + indep_ = .false. + if (present(indep)) indep_ = indep + buffer_(1) = buffer - ! Create dataspace and dataset - call h5screate_f(H5S_SCALAR_F, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_INTEGER, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_integer - -!=============================================================================== -! READ_INTEGER reads integer precision scalar data -!=============================================================================== - - subroutine read_integer(buffer, obj_id, name, indep) - integer, target, intent(inout) :: buffer - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset_id - type(c_ptr) :: f_ptr - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_integer - -!=============================================================================== -! WRITE_INTEGER_1D writes integer precision 1-D array data -!=============================================================================== + call write_int_c(group_id, 0, dims, to_c_string(name), buffer_, indep_) + end subroutine write_integer_0D subroutine write_integer_1D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id @@ -1145,127 +799,15 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(1) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_integer_1D_explicit(group_id, dims, name, buffer, indep) - else - call write_integer_1D_explicit(group_id, dims, name, buffer) - end if + dims(1) = size(buffer) + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_int_c(group_id, 1, dims, to_c_string(name), buffer, indep_) end subroutine write_integer_1D - subroutine write_integer_1D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name ! name of data - integer, intent(in), target :: buffer(dims(1)) ! data to write - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_INTEGER, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_integer_1D_explicit - -!=============================================================================== -! READ_INTEGER_1D reads integer precision 1-D array data -!=============================================================================== - - subroutine read_integer_1D(buffer, obj_id, name, indep) - integer, target, intent(inout) :: buffer(:) - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1) - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_integer_1D_explicit(dset_id, dims, buffer, indep) - else - call read_integer_1D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_integer_1D - - subroutine read_integer_1D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(1) - integer, target, intent(inout) :: buffer(dims(1)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - end subroutine read_integer_1D_explicit - -!=============================================================================== -! WRITE_INTEGER_2D writes integer precision 2-D array data -!=============================================================================== - subroutine write_integer_2D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name of data @@ -1273,127 +815,18 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(2) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_integer_2D_explicit(group_id, dims, name, buffer, indep) - else - call write_integer_2D_explicit(group_id, dims, name, buffer) - end if + ! Reverse shape of array since it will be written from C + dims(1) = size(buffer, 2) + dims(2) = size(buffer, 1) + + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_int_c(group_id, 2, dims, to_c_string(name), buffer, indep_) end subroutine write_integer_2D - subroutine write_integer_2D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(2) - character(*), intent(in) :: name ! name of data - integer, intent(in), target :: buffer(dims(1),dims(2)) - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - call h5screate_simple_f(2, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_INTEGER, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_integer_2D_explicit - -!=============================================================================== -! READ_INTEGER_2D reads integer precision 2-D array data -!=============================================================================== - - subroutine read_integer_2D(buffer, obj_id, name, indep) - integer, target, intent(inout) :: buffer(:,:) - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(2) - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_integer_2D_explicit(dset_id, dims, buffer, indep) - else - call read_integer_2D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_integer_2D - - subroutine read_integer_2D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(2) - integer, target, intent(inout) :: buffer(dims(1),dims(2)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - end subroutine read_integer_2D_explicit - -!=============================================================================== -! WRITE_INTEGER_3D writes integer precision 3-D array data -!=============================================================================== - subroutine write_integer_3D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name of data @@ -1401,127 +834,19 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(3) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_integer_3D_explicit(group_id, dims, name, buffer, indep) - else - call write_integer_3D_explicit(group_id, dims, name, buffer) - end if + ! Reverse shape of array since it will be written from C + dims(1) = size(buffer, 3) + dims(2) = size(buffer, 2) + dims(3) = size(buffer, 1) + + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_int_c(group_id, 3, dims, to_c_string(name), buffer, indep_) end subroutine write_integer_3D - subroutine write_integer_3D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(3) - character(*), intent(in) :: name ! name of data - integer, intent(in), target :: buffer(dims(1),dims(2),dims(3)) - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - call h5screate_simple_f(3, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_INTEGER, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_integer_3D_explicit - -!=============================================================================== -! READ_INTEGER_3D reads integer precision 3-D array data -!=============================================================================== - - subroutine read_integer_3D(buffer, obj_id, name, indep) - integer, target, intent(inout) :: buffer(:,:,:) - integer(HID_T), intent(in) :: obj_id - character(*), optional, intent(in) :: name - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(3) - - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from - if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) - else - dset_id = obj_id - end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_integer_3D_explicit(dset_id, dims, buffer, indep) - else - call read_integer_3D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - end subroutine read_integer_3D - - subroutine read_integer_3D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(3) - integer, target, intent(inout) :: buffer(dims(1),dims(2),dims(3)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - end subroutine read_integer_3D_explicit - -!=============================================================================== -! WRITE_INTEGER_4D writes integer precision 4-D array data -!=============================================================================== - subroutine write_integer_4D(group_id, name, buffer, indep) integer(HID_T), intent(in) :: group_id character(*), intent(in) :: name ! name of data @@ -1529,60 +854,113 @@ contains logical, intent(in), optional :: indep ! independent I/O integer(HSIZE_T) :: dims(4) + logical(C_BOOL) :: indep_ - dims(:) = shape(buffer) - if (present(indep)) then - call write_integer_4D_explicit(group_id, dims, name, buffer, indep) - else - call write_integer_4D_explicit(group_id, dims, name, buffer) - end if + ! Reverse shape of array since it will be written from C + dims(1) = size(buffer, 4) + dims(2) = size(buffer, 3) + dims(3) = size(buffer, 2) + dims(4) = size(buffer, 1) + + indep_ = .false. + if (present(indep)) indep_ = indep + + call write_int_c(group_id, 3, dims, to_c_string(name), buffer, indep_) end subroutine write_integer_4D - subroutine write_integer_4D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(4) - character(*), intent(in) :: name ! name of data - integer, intent(in), target :: buffer(dims(1),dims(2),dims(3),dims(4)) - logical, intent(in), optional :: indep ! independent I/O +!=============================================================================== +! READ_INTEGER_ND reads integer precision N-D array data +!=============================================================================== - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr + subroutine read_integer_0D(buffer, obj_id, name, indep) + integer, target, intent(inout) :: buffer + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + integer(C_INT) :: buffer_(1) + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) - call h5screate_simple_f(4, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), H5T_NATIVE_INTEGER, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) + indep_ = .false. + if (present(indep)) indep_ = indep - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_int_c(obj_id, c_loc(name_), buffer_, indep_) else - call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) + call read_int_c(obj_id, C_NULL_PTR, buffer_, indep_) end if + buffer = buffer_(1) + end subroutine read_integer_0D - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_integer_4D_explicit + subroutine read_integer_1D(buffer, obj_id, name, indep) + integer, target, intent(inout) :: buffer(:) + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O -!=============================================================================== -! READ_INTEGER_4D reads integer precision 4-D array data -!=============================================================================== + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep + + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_int_c(obj_id, c_loc(name_), buffer, indep_) + else + call read_int_c(obj_id, C_NULL_PTR, buffer, indep_) + end if + end subroutine read_integer_1D + + subroutine read_integer_2D(buffer, obj_id, name, indep) + integer, target, intent(inout) :: buffer(:,:) + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O + + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep + + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_int_c(obj_id, c_loc(name_), buffer, indep_) + else + call read_int_c(obj_id, C_NULL_PTR, buffer, indep_) + end if + end subroutine read_integer_2D + + subroutine read_integer_3D(buffer, obj_id, name, indep) + integer, target, intent(inout) :: buffer(:,:,:) + integer(HID_T), intent(in) :: obj_id + character(*), optional, intent(in) :: name + logical, optional, intent(in) :: indep ! independent I/O + + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep + + ! If 'name' argument is passed, obj_id is interpreted to be a group and + ! 'name' is the name of the dataset we should read from + if (present(name)) then + name_ = to_c_string(name) + call read_int_c(obj_id, c_loc(name_), buffer, indep_) + else + call read_int_c(obj_id, C_NULL_PTR, buffer, indep_) + end if + end subroutine read_integer_3D subroutine read_integer_4D(buffer, obj_id, name, indep) integer, target, intent(inout) :: buffer(:,:,:,:) @@ -1590,62 +968,22 @@ contains character(*), optional, intent(in) :: name logical, optional, intent(in) :: indep ! independent I/O - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(4) + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) + name_ = to_c_string(name) + call read_int_c(obj_id, c_loc(name_), buffer, indep_) else - dset_id = obj_id + call read_int_c(obj_id, C_NULL_PTR, buffer, indep_) end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_integer_4D_explicit(dset_id, dims, buffer, indep) - else - call read_integer_4D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) end subroutine read_integer_4D - subroutine read_integer_4D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(4) - integer, target, intent(inout) :: buffer(dims(1),dims(2),dims(3),dims(4)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end if - end subroutine read_integer_4D_explicit - !=============================================================================== ! WRITE_LONG writes long integer scalar data !=============================================================================== @@ -1656,40 +994,15 @@ contains integer(8), intent(in), target :: buffer ! data to write logical, intent(in), optional :: indep ! independent I/O - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr + integer(HSIZE_T) :: dims(0) + logical(C_BOOL) :: indep_ + integer(C_LONG_LONG) :: buffer_(1) - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + indep_ = .false. + if (present(indep)) indep_ = indep + buffer_(1) = buffer - ! Create dataspace and dataset - call h5screate_f(H5S_SCALAR_F, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), hdf5_integer8_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, hdf5_integer8_t, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) + call write_llong_c(group_id, 0, dims, to_c_string(name), buffer_, indep_) end subroutine write_long !=============================================================================== @@ -1702,42 +1015,22 @@ contains character(*), optional, intent(in) :: name logical, optional, intent(in) :: indep ! independent I/O - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset_id - type(c_ptr) :: f_ptr + integer(C_LONG_LONG) :: buffer_(1) + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) + name_ = to_c_string(name) + call read_llong_c(obj_id, c_loc(name_), buffer_, indep_) else - dset_id = obj_id + call read_llong_c(obj_id, C_NULL_PTR, buffer_, indep_) end if - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr = c_loc(buffer) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, hdf5_integer8_t, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, hdf5_integer8_t, f_ptr, hdf5_err) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) + buffer = buffer_(1) end subroutine read_long !=============================================================================== @@ -1750,58 +1043,16 @@ contains character(*), intent(in), target :: buffer ! read data to here logical, intent(in), optional :: indep ! independent I/O - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - integer(HID_T) :: filetype - integer(SIZE_T) :: i, n - character(kind=C_CHAR), allocatable, target :: temp_buffer(:) - type(c_ptr) :: f_ptr + integer(HSIZE_T) :: dims(0) + integer(C_SIZE_T) :: slen + logical(C_BOOL) :: indep_ - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + indep_ = .false. + if (present(indep)) indep_ = indep + slen = len_trim(buffer) - ! Create datatype for HDF5 file based on C char - n = len_trim(buffer) - if (n > 0) then - call h5tcopy_f(H5T_C_S1, filetype, hdf5_err) - call h5tset_size_f(filetype, n, hdf5_err) - - ! Create dataspace/dataset - call h5screate_f(H5S_SCALAR_F, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), filetype, dspace, dset, hdf5_err) - - ! Copy string to temporary buffer - allocate(temp_buffer(n)) - do i = 1, n - temp_buffer(i) = buffer(i:i) - end do - - ! Get pointer to start of string - f_ptr = c_loc(temp_buffer(1)) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dwrite_f(dset, filetype, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dwrite_f(dset, filetype, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - call h5tclose_f(filetype, hdf5_err) - end if + call write_string_c(group_id, 0, dims, slen, to_c_string(name), & + to_c_string(buffer), indep_) end subroutine write_string !=============================================================================== @@ -1814,67 +1065,32 @@ contains character(*), optional, intent(in) :: name logical, optional, intent(in) :: indep ! independent I/O - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif integer(HID_T) :: dset_id - integer(HID_T) :: filetype - integer(HID_T) :: memtype - integer(SIZE_T) :: i, n - character(kind=C_CHAR), allocatable, target :: temp_buffer(:) - type(c_ptr) :: f_ptr + integer(C_SIZE_T) :: i, n + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), allocatable, target :: buffer_(:) if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) + dset_id = open_dataset(obj_id, name) else dset_id = obj_id end if - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if + ! Allocate a C char array to get string + n = dataset_typesize(dset_id) + allocate(buffer_(n)) - ! Make sure buffer is large enough - call h5dget_type_f(dset_id, filetype, hdf5_err) - call h5tget_size_f(filetype, n, hdf5_err) - if (n > len(buffer)) then - call fatal_error("Character buffer is not long enough to & - &read HDF5 string.") - end if - - ! Get datatype in memory based on Fortran character - call h5tcopy_f(H5T_C_S1, memtype, hdf5_err) - call h5tset_size_f(memtype, n + 1, hdf5_err) - - ! Get pointer to start of string - allocate(temp_buffer(n)) - f_ptr = c_loc(temp_buffer(1)) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, memtype, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, memtype, f_ptr, hdf5_err) - end if + indep_ = .false. + if (present(indep)) indep_ = indep + call read_string_c(dset_id, C_NULL_PTR, n, buffer_, indep_) buffer = '' do i = 1, n - if (temp_buffer(i) == C_NULL_CHAR) cycle - buffer(i:i) = temp_buffer(i) + if (buffer_(i) == C_NULL_CHAR) cycle + buffer(i:i) = buffer_(i) end do - if (present(name)) call h5dclose_f(dset_id, hdf5_err) - - call h5tclose_f(filetype, hdf5_err) - call h5tclose_f(memtype, hdf5_err) + call close_dataset(dset_id) end subroutine read_string !=============================================================================== @@ -1887,74 +1103,29 @@ contains character(*), intent(in), target :: buffer(:) ! read data to here logical, intent(in), optional :: indep ! independent I/O + integer(HSIZE_T) :: i integer(HSIZE_T) :: dims(1) + integer(C_SIZE_T) :: m, n + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), allocatable :: buffer_(:) - dims(:) = shape(buffer) - if (present(indep)) then - call write_string_1D_explicit(group_id, dims, name, buffer, indep) - else - call write_string_1D_explicit(group_id, dims, name, buffer) - end if + indep_ = .false. + if (present(indep)) indep_ = indep + + ! Copy array of characters into an array of C chars with the right memory + ! layout + dims(1) = size(buffer) + m = maxval(len_trim(buffer)) + 1 + allocate(buffer_(dims(1)*m)) + do i = 0, dims(1) - 1 + n = len_trim(buffer(i+1)) + 1 + buffer_(i*m+1 : i*m+n) = to_c_string(buffer(i+1)) + end do + + call write_string_c(group_id, 1, dims, m, to_c_string(name), & + buffer_, indep_) end subroutine write_string_1D - subroutine write_string_1D_explicit(group_id, dims, name, buffer, indep) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name - character(*), intent(in), target :: buffer(dims(1)) - logical, intent(in), optional :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - integer(HID_T) :: filetype - integer(HID_T) :: memtype - integer(SIZE_T) :: n - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - ! Create datatype for HDF5 file based on C char - n = maxval(len_trim(buffer)) - call h5tcopy_f(H5T_C_S1, filetype, hdf5_err) - call h5tset_size_f(filetype, n + 1, hdf5_err) - - ! Create datatype in memory based on Fortran character - call h5tcopy_f(H5T_FORTRAN_S1, memtype, hdf5_err) - call h5tset_size_f(memtype, int(len(buffer(1)), SIZE_T), hdf5_err) - - ! Create dataspace/dataset - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), filetype, dspace, dset, hdf5_err) - - ! Get pointer to start of string - f_ptr = c_loc(buffer(1)(1:1)) - - if (using_mpio_device(group_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - if (n > 0) call h5dwrite_f(dset, memtype, f_ptr, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - if (n > 0) call h5dwrite_f(dset, memtype, f_ptr, hdf5_err) - end if - - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - call h5tclose_f(memtype, hdf5_err) - call h5tclose_f(filetype, hdf5_err) - end subroutine write_string_1D_explicit - !=============================================================================== ! READ_STRING_1D reads string 1-D array data !=============================================================================== @@ -1965,87 +1136,39 @@ contains character(*), optional, intent(in) :: name logical, optional, intent(in) :: indep ! independent I/O - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1) + integer(HID_T) :: dset_id + integer(C_SIZE_T) :: i, j, k, n, m + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), allocatable, target :: buffer_(:) - ! If 'name' argument is passed, obj_id is interpreted to be a group and - ! 'name' is the name of the dataset we should read from if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) + dset_id = open_dataset(obj_id, name) else dset_id = obj_id end if - dims(:) = shape(buffer) + ! Allocate a C char array to get strings + n = dataset_typesize(dset_id) + m = size(buffer) + allocate(buffer_(n*m)) - if (present(indep)) then - call read_string_1D_explicit(dset_id, dims, buffer, indep) - else - call read_string_1D_explicit(dset_id, dims, buffer) - end if + indep_ = .false. + if (present(indep)) indep_ = indep + call read_string_c(dset_id, C_NULL_PTR, n, buffer_, indep_) + + ! Convert null-terminated C strings into Fortran strings + do i = 1, m + buffer(i) = '' + do j = 1, n + k = (i-1)*n + j + if (buffer_(k) == C_NULL_CHAR) exit + buffer(i)(j:j) = buffer_(k) + end do + end do + + call close_dataset(dset_id) end subroutine read_string_1D - subroutine read_string_1D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), target, intent(inout) :: buffer(dims(1)) - logical, optional, intent(in) :: indep ! independent I/O - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - integer(HID_T) :: space_id - integer(HID_T) :: filetype - integer(HID_T) :: memtype - integer(SIZE_T) :: size - integer(SIZE_T) :: n - type(c_ptr) :: f_ptr - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - ! Get dataset and dataspace - call h5dget_space_f(dset_id, space_id, hdf5_err) - - ! Make sure buffer is large enough - call h5dget_type_f(dset_id, filetype, hdf5_err) - call h5tget_size_f(filetype, size, hdf5_err) - if (size > len(buffer(1)) + 1) then - call fatal_error("Character buffer is not long enough to & - &read HDF5 string array.") - end if - - ! Get datatype in memory based on Fortran character - n = len(buffer(1)) - call h5tcopy_f(H5T_FORTRAN_S1, memtype, hdf5_err) - call h5tset_size_f(memtype, n, hdf5_err) - - ! Get pointer to start of string - f_ptr = c_loc(buffer(1)(1:1)) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, memtype, f_ptr, hdf5_err, mem_space_id=space_id, & - xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, memtype, f_ptr, hdf5_err, mem_space_id=space_id) - end if - - call h5sclose_f(space_id, hdf5_err) - call h5tclose_f(filetype, hdf5_err) - call h5tclose_f(memtype, hdf5_err) - end subroutine read_string_1D_explicit - !=============================================================================== ! WRITE_ATTRIBUTE_STRING !=============================================================================== @@ -2055,41 +1178,8 @@ contains character(*), intent(in) :: name ! name of attribute character(*), intent(in), target :: buffer ! string to write - integer :: hdf5_err - integer(HID_T) :: dspace_id - integer(HID_T) :: attr_id - integer(HID_T) :: filetype - integer(SIZE_T) :: i - integer(SIZE_T) :: n - character(kind=C_CHAR), allocatable, target :: temp_buffer(:) - type(c_ptr) :: f_ptr - ! Create datatype for HDF5 file based on C char - n = len_trim(buffer) - if (n > 0) then - call h5tcopy_f(H5T_C_S1, filetype, hdf5_err) - call h5tset_size_f(filetype, n, hdf5_err) - - ! Create memory space and attribute - call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err) - call h5acreate_f(obj_id, trim(name), filetype, dspace_id, & - attr_id, hdf5_err) - - ! Copy string to temporary buffer - allocate(temp_buffer(n)) - do i = 1, n - temp_buffer(i) = buffer(i:i) - end do - - ! Write attribute - f_ptr = c_loc(buffer(1:1)) - call h5awrite_f(attr_id, filetype, f_ptr, hdf5_err) - - ! Close attribute - call h5aclose_f(attr_id, hdf5_err) - call h5sclose_f(dspace_id, hdf5_err) - call h5tclose_f(filetype, hdf5_err) - end if + call write_attr_string_c(obj_id, to_c_string(name), to_c_string(buffer)) end subroutine write_attribute_string subroutine read_attribute_double(buffer, obj_id, name) @@ -2097,14 +1187,10 @@ contains integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: attr_id - type(c_ptr) :: f_ptr + real(C_DOUBLE) :: buffer_(1) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - f_ptr = c_loc(buffer) - call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) + call read_attr_double_c(obj_id, to_c_string(name), buffer_) + buffer = buffer_(1) end subroutine read_attribute_double subroutine write_attribute_double(obj_id, name, buffer) @@ -2112,18 +1198,11 @@ contains character(*), intent(in) :: name real(8), intent(in), target :: buffer - integer :: hdf5_err - integer(HID_T) :: dspace_id - integer(HID_T) :: attr_id - type(C_PTR) :: f_ptr + integer(HSIZE_T) :: dims(0) + real(C_DOUBLE) :: buffer_(1) - call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err) - call h5acreate_f(obj_id, trim(name), H5T_NATIVE_DOUBLE, dspace_id, & - attr_id, hdf5_err) - f_ptr = c_loc(buffer) - call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) - call h5sclose_f(dspace_id, hdf5_err) + buffer_(1) = buffer + call write_attr_double_c(obj_id, 0, dims, to_c_string(name), buffer_) end subroutine write_attribute_double subroutine read_attribute_double_1D(buffer, obj_id, name) @@ -2131,39 +1210,16 @@ contains integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: space_id - integer(HID_T) :: attr_id integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: maxdims(1) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - - if (allocated(buffer)) then - dims(:) = shape(buffer) - else - call h5aget_space_f(attr_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) + if (.not. allocated(buffer)) then + call get_shape_attr(obj_id, to_c_string(name), dims) allocate(buffer(dims(1))) - call h5sclose_f(space_id, hdf5_err) end if - call read_attribute_double_1D_explicit(attr_id, dims, buffer) - call h5aclose_f(attr_id, hdf5_err) + call read_attr_double_c(obj_id, to_c_string(name), buffer) end subroutine read_attribute_double_1D - subroutine read_attribute_double_1D_explicit(attr_id, dims, buffer) - integer(HID_T), intent(in) :: attr_id - integer(HSIZE_T), intent(in) :: dims(1) - real(8), target, intent(inout) :: buffer(dims(1)) - - integer :: hdf5_err - type(c_ptr) :: f_ptr - - f_ptr = c_loc(buffer) - call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end subroutine read_attribute_double_1D_explicit - subroutine write_attribute_double_1D(obj_id, name, buffer) integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name @@ -2171,81 +1227,34 @@ contains integer(HSIZE_T) :: dims(1) - dims(:) = shape(buffer) - call write_attribute_double_1D_explicit(obj_id, dims, name, buffer) + dims(1) = size(buffer) + call write_attr_double_c(obj_id, 1, dims, to_c_string(name), buffer) end subroutine write_attribute_double_1D - subroutine write_attribute_double_1D_explicit(obj_id, dims, name, buffer) - integer(HID_T), intent(in) :: obj_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name - real(8), target, intent(in) :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: dspace_id - integer(HID_T) :: attr_id - type(C_PTR) :: f_ptr - - call h5screate_simple_f(1, dims, dspace_id, hdf5_err) - call h5acreate_f(obj_id, trim(name), H5T_NATIVE_DOUBLE, dspace_id, & - attr_id, hdf5_err) - f_ptr = c_loc(buffer) - call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) - call h5sclose_f(dspace_id, hdf5_err) - end subroutine write_attribute_double_1D_explicit - subroutine read_attribute_double_2D(buffer, obj_id, name) real(8), target, allocatable, intent(inout) :: buffer(:,:) integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: space_id - integer(HID_T) :: attr_id integer(HSIZE_T) :: dims(2) - integer(HSIZE_T) :: maxdims(2) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - - if (allocated(buffer)) then - dims(:) = shape(buffer) - else - call h5aget_space_f(attr_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) - allocate(buffer(dims(1), dims(2))) - call h5sclose_f(space_id, hdf5_err) + if (.not. allocated(buffer)) then + call get_shape_attr(obj_id, to_c_string(name), dims) + allocate(buffer(dims(2), dims(1))) end if - call read_attribute_double_2D_explicit(attr_id, dims, buffer) - call h5aclose_f(attr_id, hdf5_err) + call read_attr_double_c(obj_id, to_c_string(name), buffer) end subroutine read_attribute_double_2D - subroutine read_attribute_double_2D_explicit(attr_id, dims, buffer) - integer(HID_T), intent(in) :: attr_id - integer(HSIZE_T), intent(in) :: dims(2) - real(8), target, intent(inout) :: buffer(dims(1),dims(2)) - - integer :: hdf5_err - type(c_ptr) :: f_ptr - - f_ptr = c_loc(buffer) - call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) - end subroutine read_attribute_double_2D_explicit - subroutine read_attribute_integer(buffer, obj_id, name) integer, intent(inout), target :: buffer integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: attr_id - type(c_ptr) :: f_ptr + integer(C_INT) :: buffer_(1) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - f_ptr = c_loc(buffer) - call h5aread_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) + call read_attr_int_c(obj_id, to_c_string(name), buffer_) + buffer = buffer_(1) end subroutine read_attribute_integer subroutine write_attribute_integer(obj_id, name, buffer) @@ -2253,18 +1262,11 @@ contains character(*), intent(in) :: name integer, intent(in), target :: buffer - integer :: hdf5_err - integer(HID_T) :: dspace_id - integer(HID_T) :: attr_id - type(C_PTR) :: f_ptr + integer(HSIZE_T) :: dims(0) + integer(C_INT) :: buffer_(1) - call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err) - call h5acreate_f(obj_id, trim(name), H5T_NATIVE_INTEGER, dspace_id, & - attr_id, hdf5_err) - f_ptr = c_loc(buffer) - call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) - call h5sclose_f(dspace_id, hdf5_err) + buffer_(1) = buffer + call write_attr_int_c(obj_id, 0, dims, to_c_string(name), buffer_) end subroutine write_attribute_integer subroutine read_attribute_integer_1D(buffer, obj_id, name) @@ -2272,39 +1274,16 @@ contains integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: space_id - integer(HID_T) :: attr_id integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: maxdims(1) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - - if (allocated(buffer)) then - dims(:) = shape(buffer) - else - call h5aget_space_f(attr_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) + if (.not. allocated(buffer)) then + call get_shape_attr(obj_id, to_c_string(name), dims) allocate(buffer(dims(1))) - call h5sclose_f(space_id, hdf5_err) end if - call read_attribute_integer_1D_explicit(attr_id, dims, buffer) - call h5aclose_f(attr_id, hdf5_err) + call read_attr_int_c(obj_id, to_c_string(name), buffer) end subroutine read_attribute_integer_1D - subroutine read_attribute_integer_1D_explicit(attr_id, dims, buffer) - integer(HID_T), intent(in) :: attr_id - integer(HSIZE_T), intent(in) :: dims(1) - integer, target, intent(inout) :: buffer(dims(1)) - - integer :: hdf5_err - type(c_ptr) :: f_ptr - - f_ptr = c_loc(buffer) - call h5aread_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end subroutine read_attribute_integer_1D_explicit - subroutine write_attribute_integer_1D(obj_id, name, buffer) integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name @@ -2312,111 +1291,46 @@ contains integer(HSIZE_T) :: dims(1) - dims(:) = shape(buffer) - call write_attribute_integer_1D_explicit(obj_id, dims, name, buffer) + dims(1) = size(buffer) + call write_attr_int_c(obj_id, 1, dims, to_c_string(name), buffer) end subroutine write_attribute_integer_1D - subroutine write_attribute_integer_1D_explicit(obj_id, dims, name, buffer) - integer(HID_T), intent(in) :: obj_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name - integer, target, intent(in) :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: dspace_id - integer(HID_T) :: attr_id - type(C_PTR) :: f_ptr - - call h5screate_simple_f(1, dims, dspace_id, hdf5_err) - call h5acreate_f(obj_id, trim(name), H5T_NATIVE_INTEGER, dspace_id, & - attr_id, hdf5_err) - f_ptr = c_loc(buffer) - call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) - call h5sclose_f(dspace_id, hdf5_err) - end subroutine write_attribute_integer_1D_explicit - subroutine read_attribute_integer_2D(buffer, obj_id, name) integer, target, allocatable, intent(inout) :: buffer(:,:) integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: space_id - integer(HID_T) :: attr_id integer(HSIZE_T) :: dims(2) - integer(HSIZE_T) :: maxdims(2) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - - if (allocated(buffer)) then - dims(:) = shape(buffer) - else - call h5aget_space_f(attr_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) - allocate(buffer(dims(1), dims(2))) - call h5sclose_f(space_id, hdf5_err) + if (.not. allocated(buffer)) then + call get_shape_attr(obj_id, to_c_string(name), dims) + allocate(buffer(dims(2), dims(1))) end if - call read_attribute_integer_2D_explicit(attr_id, dims, buffer) - call h5aclose_f(attr_id, hdf5_err) + call read_attr_int_c(obj_id, to_c_string(name), buffer) end subroutine read_attribute_integer_2D - subroutine read_attribute_integer_2D_explicit(attr_id, dims, buffer) - integer(HID_T), intent(in) :: attr_id - integer(HSIZE_T), intent(in) :: dims(2) - integer, target, intent(inout) :: buffer(dims(1),dims(2)) - - integer :: hdf5_err - type(c_ptr) :: f_ptr - - f_ptr = c_loc(buffer) - call h5aread_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - end subroutine read_attribute_integer_2D_explicit - subroutine read_attribute_string(buffer, obj_id, name) character(*), intent(inout) :: buffer ! read data to here integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name ! name for data - integer :: hdf5_err - integer(HID_T) :: attr_id ! data set handle - integer(HID_T) :: filetype - integer(HID_T) :: memtype - integer(SIZE_T) :: i - integer(SIZE_T) :: size - character(kind=C_CHAR), allocatable, target :: temp_buffer(:) - type(c_ptr) :: f_ptr + integer(C_SIZE_T) :: i, n + character(kind=C_CHAR), allocatable, target :: buffer_(:) - ! Get dataset and dataspace - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) + ! Allocate a C char array to get string + n = attribute_typesize(obj_id, to_c_string(name)) + allocate(buffer_(n)) - ! Make sure buffer is large enough - call h5aget_type_f(attr_id, filetype, hdf5_err) - call h5tget_size_f(filetype, size, hdf5_err) - allocate(temp_buffer(size)) - if (size > len(buffer)) then - call fatal_error("Character buffer is not long enough to & - &read HDF5 string.") - end if + ! Read attribute + call read_attr_string_c(obj_id, to_c_string(name), n, buffer_) - ! Get datatype in memory based on Fortran character - call h5tcopy_f(H5T_C_S1, memtype, hdf5_err) - call h5tset_size_f(memtype, size + 1, hdf5_err) - - ! Get pointer to start of string - f_ptr = c_loc(temp_buffer(1)) - - call h5aread_f(attr_id, memtype, f_ptr, hdf5_err) + ! Copy back to Fortran string buffer = '' - do i = 1, size - if (temp_buffer(i) == C_NULL_CHAR) cycle - buffer(i:i) = temp_buffer(i) + do i = 1, n + if (buffer_(i) == C_NULL_CHAR) cycle + buffer(i:i) = buffer_(i) end do - - call h5aclose_f(attr_id, hdf5_err) - call h5tclose_f(filetype, hdf5_err) - call h5tclose_f(memtype, hdf5_err) end subroutine read_attribute_string subroutine read_attribute_string_1D(buffer, obj_id, name) @@ -2424,247 +1338,98 @@ contains integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer :: hdf5_err - integer(HID_T) :: space_id - integer(HID_T) :: attr_id + integer(C_SIZE_T) :: i, j, k, n, m integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: maxdims(1) + character(kind=C_CHAR), allocatable, target :: buffer_(:) - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - - if (allocated(buffer)) then - dims(:) = shape(buffer) - else - call h5aget_space_f(attr_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) + if (.not. allocated(buffer)) then + call get_shape_attr(obj_id, to_c_string(name), dims) allocate(buffer(dims(1))) - call h5sclose_f(space_id, hdf5_err) end if - call read_attribute_string_1D_explicit(attr_id, dims, buffer) - call h5aclose_f(attr_id, hdf5_err) + ! Allocate a C char array to get strings + n = attribute_typesize(obj_id, to_c_string(name)) + m = size(buffer) + allocate(buffer_((n+1)*m)) + + ! Read attribute + call read_attr_string_c(obj_id, to_c_string(name), n, buffer_) + + ! Convert null-terminated C strings into Fortran strings + do i = 1, m + buffer(i) = '' + do j = 1, n + k = (i-1)*(n+1) + j + if (buffer_(k) == C_NULL_CHAR) exit + buffer(i)(j:j) = buffer_(k) + end do + end do end subroutine read_attribute_string_1D - subroutine read_attribute_string_1D_explicit(attr_id, dims, buffer) - integer(HID_T), intent(in) :: attr_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), target, intent(inout) :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: filetype - integer(HID_T) :: memtype - integer(SIZE_T) :: size - integer(SIZE_T) :: n - type(c_ptr) :: f_ptr - - ! Make sure buffer is large enough - call h5aget_type_f(attr_id, filetype, hdf5_err) - call h5tget_size_f(filetype, size, hdf5_err) - if (size > len(buffer(1)) + 1) then - call fatal_error("Character buffer is not long enough to & - &read HDF5 string array.") - end if - - ! Get datatype in memory based on Fortran character - n = len(buffer(1)) - call h5tcopy_f(H5T_FORTRAN_S1, memtype, hdf5_err) - call h5tset_size_f(memtype, n, hdf5_err) - - ! Get pointer to start of string - f_ptr = c_loc(buffer(1)(1:1)) - - call h5aread_f(attr_id, memtype, f_ptr, hdf5_err) - - call h5tclose_f(filetype, hdf5_err) - call h5tclose_f(memtype, hdf5_err) - end subroutine read_attribute_string_1D_explicit - subroutine read_attribute_logical(buffer, obj_id, name) logical, intent(inout), target :: buffer integer(HID_T), intent(in) :: obj_id character(*), intent(in) :: name - integer, target :: int_buffer - integer :: hdf5_err - integer(HID_T) :: attr_id - type(c_ptr) :: f_ptr + integer :: tmp - call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) - f_ptr = c_loc(int_buffer) - call h5aread_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) - call h5aclose_f(attr_id, hdf5_err) - - ! Convert to Fortran logical - if (int_buffer == 0) then - buffer = .false. - else - buffer = .true. - end if + call read_attribute_integer(tmp, obj_id, name) + buffer = (tmp /= 0) end subroutine read_attribute_logical subroutine get_shape(obj_id, dims) integer(HID_T), intent(in) :: obj_id integer(HSIZE_T), intent(out) :: dims(:) - integer :: hdf5_err - integer :: type - integer(HID_T) :: space_id - integer(HSIZE_T) :: maxdims(size(dims)) + integer :: i + integer(HSIZE_T) :: dims_c(size(dims)) - call h5iget_type_f(obj_id, type, hdf5_err) - if (type == H5I_DATASET_F) then - call h5dget_space_f(obj_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) - call h5sclose_f(space_id, hdf5_err) - elseif (type == H5I_ATTR_F) then - call h5aget_space_f(obj_id, space_id, hdf5_err) - call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) - call h5sclose_f(space_id, hdf5_err) - end if + call get_shape_c(obj_id, dims_c) + do i = 1, size(dims) + dims(i) = dims_c(size(dims) - i + 1) + end do end subroutine get_shape subroutine get_ndims(obj_id, ndims) integer(HID_T), intent(in) :: obj_id integer, intent(out) :: ndims - integer :: hdf5_err - integer :: type - integer(HID_T) :: space_id + interface + function dataset_ndims(dset) result(ndims) bind(C) + import HID_T, C_INT + integer(HID_T), value :: dset + integer(C_INT) :: ndims + end function dataset_ndims + end interface - call h5iget_type_f(obj_id, type, hdf5_err) - if (type == H5I_DATASET_F) then - call h5dget_space_f(obj_id, space_id, hdf5_err) - call h5sget_simple_extent_ndims_f(space_id, ndims, hdf5_err) - call h5sclose_f(space_id, hdf5_err) - elseif (type == H5I_ATTR_F) then - call h5aget_space_f(obj_id, space_id, hdf5_err) - call h5sget_simple_extent_ndims_f(space_id, ndims, hdf5_err) - call h5sclose_f(space_id, hdf5_err) - end if + ndims = dataset_ndims(obj_id) end subroutine get_ndims - function using_mpio_device(obj_id) result(mpio) - integer(HID_T), intent(in) :: obj_id - logical :: mpio - - integer :: hdf5_err - integer(HID_T) :: driver - integer(HID_T) :: file_id - integer(HID_T) :: fapl_id - - ! Determine file that this object is part of - call h5iget_file_id_f(obj_id, file_id, hdf5_err) - - ! Get file access property list - call h5fget_access_plist_f(file_id, fapl_id, hdf5_err) - - ! Get low-level driver identifier - call h5pget_driver_f(fapl_id, driver, hdf5_err) - - ! Close file access property list access - call h5pclose_f(fapl_id, hdf5_err) - - ! Close file access -- note that this only decreases the reference count so - ! that the file is not actually closed - call h5fclose_f(file_id, hdf5_err) - - mpio = (driver == H5FD_MPIO_F) - end function using_mpio_device - !=============================================================================== ! READ_COMPLEX_2D reads double precision complex 2-D array data as output by ! the h5py HDF5 python module. !=============================================================================== subroutine read_complex_2D(buffer, obj_id, name, indep) - complex(8), target, intent(inout) :: buffer(:,:) + complex(C_DOUBLE_COMPLEX), target, intent(inout) :: buffer(:,:) integer(HID_T), intent(in) :: obj_id character(*), optional, intent(in) :: name logical, optional, intent(in) :: indep ! independent I/O - integer :: hdf5_err - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(2) + logical(C_BOOL) :: indep_ + character(kind=C_CHAR), target, allocatable :: name_(:) + + indep_ = .false. + if (present(indep)) indep_ = indep ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then - call h5dopen_f(obj_id, trim(name), dset_id, hdf5_err) + name_ = to_c_string(name) + call read_complex_c(obj_id, c_loc(name_), buffer, indep_) else - dset_id = obj_id + call read_complex_c(obj_id, C_NULL_PTR, buffer, indep_) end if - - dims(:) = shape(buffer) - - if (present(indep)) then - call read_complex_2D_explicit(dset_id, dims, buffer, indep) - else - call read_complex_2D_explicit(dset_id, dims, buffer) - end if - - if (present(name)) call h5dclose_f(dset_id, hdf5_err) end subroutine read_complex_2D - subroutine read_complex_2D_explicit(dset_id, dims, buffer, indep) - integer(HID_T), intent(in) :: dset_id - integer(HSIZE_T), intent(in) :: dims(2) - complex(8), target, intent(inout) :: buffer(dims(1), dims(2)) - logical, optional, intent(in) :: indep ! independent I/O - - real(8), target :: buffer_r(dims(1), dims(2)) - real(8), target :: buffer_i(dims(1), dims(2)) - - integer(HSIZE_T) :: i, j - - integer :: hdf5_err - integer :: data_xfer_mode -#ifdef PHDF5 - integer(HID_T) :: plist ! property list -#endif - type(c_ptr) :: f_ptr_r, f_ptr_i - - ! Components needed for complex type support - integer(HID_T) :: dtype_real - integer(HID_T) :: dtype_imag - integer(SIZE_T) :: size_double - - ! Create the complex type - call h5tget_size_f(H5T_NATIVE_DOUBLE, size_double, hdf5_err) - - ! Insert the 'r' and 'i' identifiers - call h5tcreate_f(H5T_COMPOUND_F, size_double, dtype_real, hdf5_err) - call h5tcreate_f(H5T_COMPOUND_F, size_double, dtype_imag, hdf5_err) - call h5tinsert_f(dtype_real, "r", 0_SIZE_T, H5T_NATIVE_DOUBLE, hdf5_err) - call h5tinsert_f(dtype_imag, "i", 0_SIZE_T, H5T_NATIVE_DOUBLE, hdf5_err) - - ! Set up collective vs. independent I/O - data_xfer_mode = H5FD_MPIO_COLLECTIVE_F - if (present(indep)) then - if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F - end if - - f_ptr_r = c_loc(buffer_r) - f_ptr_i = c_loc(buffer_i) - - if (using_mpio_device(dset_id)) then -#ifdef PHDF5 - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err) - call h5dread_f(dset_id, dtype_real, f_ptr_r, hdf5_err, xfer_prp=plist) - call h5dread_f(dset_id, dtype_imag, f_ptr_i, hdf5_err, xfer_prp=plist) - call h5pclose_f(plist, hdf5_err) -#endif - else - call h5dread_f(dset_id, dtype_real, f_ptr_r, hdf5_err) - call h5dread_f(dset_id, dtype_imag, f_ptr_i, hdf5_err) - end if - - ! Reconstitute the complex numbers - do i = 1, dims(1) - do j = 1, dims(2) - buffer(i, j) = cmplx(buffer_r(i,j), buffer_i(i,j), kind=8) - end do - end do - end subroutine read_complex_2D_explicit - end module hdf5_interface diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp new file mode 100644 index 000000000..8a8391bb5 --- /dev/null +++ b/src/hdf5_interface.cpp @@ -0,0 +1,658 @@ +#include "hdf5_interface.h" + +#include +#include +#include +#include + +#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 diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index 7ec9ada9b..e38a31e99 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -1,90 +1,96 @@ #ifndef HDF5_INTERFACE_H #define HDF5_INTERFACE_H -#include -#include - #include "hdf5.h" +#include "hdf5_hl.h" -#include "error.h" +#include +#include +#include +#include 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 void write_double_1D(hid_t group_id, char const *name, std::array &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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 7b699b9cb..d5ce32b31 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/initialize.cpp b/src/initialize.cpp new file mode 100644 index 000000000..728f43695 --- /dev/null +++ b/src/initialize.cpp @@ -0,0 +1,229 @@ +#include "initialize.h" + +#include +#include +#include +#include +#include + +#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(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 diff --git a/src/initialize.h b/src/initialize.h new file mode 100644 index 000000000..14283fb5f --- /dev/null +++ b/src/initialize.h @@ -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 diff --git a/src/main.F90 b/src/main.F90 deleted file mode 100644 index 372e993f5..000000000 --- a/src/main.F90 +++ /dev/null @@ -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 diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 000000000..8d8f39e24 --- /dev/null +++ b/src/main.cpp @@ -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 +} diff --git a/src/math.F90 b/src/math.F90 index 9e1f01de5..6534b66b0 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -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) = [& diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index fc13a868e..fb22d2ad3 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -2,8 +2,6 @@ module mesh_header use, intrinsic :: ISO_C_BINDING - use hdf5 - use constants use dict_header, only: DictIntInt use error diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 7391a632c..6ade0895a 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -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 diff --git a/src/message_passing.cpp b/src/message_passing.cpp new file mode 100644 index 000000000..2ff3952d7 --- /dev/null +++ b/src/message_passing.cpp @@ -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 diff --git a/src/message_passing.h b/src/message_passing.h new file mode 100644 index 000000000..14cf3a7cb --- /dev/null +++ b/src/message_passing.h @@ -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 diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 1b14e6e9f..e0631d16c 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -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 !=============================================================================== diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 9ffb83815..691922adc 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -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 diff --git a/src/multipole_header.F90 b/src/multipole_header.F90 index 7fb438bce..f2d0c9106 100644 --- a/src/multipole_header.F90 +++ b/src/multipole_header.F90 @@ -1,7 +1,5 @@ module multipole_header - use hdf5 - use constants use dict_header, only: DictIntInt use error, only: fatal_error diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8362f5d41..ebabe5cdb 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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:)) diff --git a/src/output.F90 b/src/output.F90 index b6809ca82..174292fb3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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]' diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 552f9b042..cdb97dec8 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -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 diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 2d15f2ddc..c3d25151b 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -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 diff --git a/src/plot.F90 b/src/plot.F90 index 5eb894ce8..f17ac3351 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -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 diff --git a/src/plot.cpp b/src/plot.cpp new file mode 100644 index 000000000..0dca6d08c --- /dev/null +++ b/src/plot.cpp @@ -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 diff --git a/src/plot.h b/src/plot.h new file mode 100644 index 000000000..75406b9ac --- /dev/null +++ b/src/plot.h @@ -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 diff --git a/src/product_header.F90 b/src/product_header.F90 index a69929473..a5e23f4b6 100644 --- a/src/product_header.F90 +++ b/src/product_header.F90 @@ -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 diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 12b3888c4..48973336b 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -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 diff --git a/src/sab_header.F90 b/src/sab_header.F90 index eca8e555a..615e7b281 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -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:)) diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 index a0e203f33..6f5e6cabf 100644 --- a/src/secondary_correlated.F90 +++ b/src/secondary_correlated.F90 @@ -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 !=============================================================================== diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 index f963cff3f..3315174e1 100644 --- a/src/secondary_kalbach.F90 +++ b/src/secondary_kalbach.F90 @@ -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 !=============================================================================== diff --git a/src/secondary_nbody.F90 b/src/secondary_nbody.F90 index 14ae949a9..aee4b1ff7 100644 --- a/src/secondary_nbody.F90 +++ b/src/secondary_nbody.F90 @@ -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 diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 index 40aa827da..e559c1e5a 100644 --- a/src/secondary_uncorrelated.F90 +++ b/src/secondary_uncorrelated.F90 @@ -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 !=============================================================================== diff --git a/src/settings.F90 b/src/settings.F90 index 83aebce86..aef1c011b 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -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. diff --git a/src/simulation.F90 b/src/simulation.F90 index 68debe52b..4f6163429 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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 diff --git a/src/simulation.cpp b/src/simulation.cpp new file mode 100644 index 000000000..e9bf81775 --- /dev/null +++ b/src/simulation.cpp @@ -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; +} diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 11be9bf9c..8fe5e2da1 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -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 diff --git a/src/source.F90 b/src/source.F90 index 5beecd887..368361131 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index 980a7333c..a6e44dfa1 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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 diff --git a/src/state_point.cpp b/src/state_point.cpp new file mode 100644 index 000000000..4e6b9a9f9 --- /dev/null +++ b/src/state_point.cpp @@ -0,0 +1,165 @@ +#include "state_point.h" + +#include +#include + +#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(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(openmc_work)}; + hid_t memspace = H5Screate_simple(1, count, nullptr); + + // Select hyperslab for this dataspace + hsize_t start[] {static_cast(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(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 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(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(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(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(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 diff --git a/src/state_point.h b/src/state_point.h new file mode 100644 index 000000000..459df1a67 --- /dev/null +++ b/src/state_point.h @@ -0,0 +1,17 @@ +#ifndef STATE_POINT_H +#define STATE_POINT_H + +#include + +#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 diff --git a/src/summary.F90 b/src/summary.F90 index 3aeb42178..6f8333b2c 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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) diff --git a/src/surface.cpp b/src/surface.cpp index 9abfd23d2..2c3279604 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -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 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 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 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 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 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 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 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 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 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 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 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 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 diff --git a/src/surface.h b/src/surface.h index 627a8fbae..15ac0d5bd 100644 --- a/src/surface.h +++ b/src/surface.h @@ -25,15 +25,14 @@ extern "C" const int BC_PERIODIC {3}; //============================================================================== extern "C" double FP_COINCIDENT; -constexpr double INFTY{std::numeric_limits::max()}; +constexpr double INFTY {std::numeric_limits::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 diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 9ed89a14e..ab6babb5c 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -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 diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index e484ac231..4b589cedb 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -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 diff --git a/src/tallies/tally_filter_azimuthal.F90 b/src/tallies/tally_filter_azimuthal.F90 index 272d111d1..fc7f15a17 100644 --- a/src/tallies/tally_filter_azimuthal.F90 +++ b/src/tallies/tally_filter_azimuthal.F90 @@ -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 diff --git a/src/tallies/tally_filter_cell.F90 b/src/tallies/tally_filter_cell.F90 index 8d2b93d28..00ab8fa15 100644 --- a/src/tallies/tally_filter_cell.F90 +++ b/src/tallies/tally_filter_cell.F90 @@ -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 diff --git a/src/tallies/tally_filter_cellborn.F90 b/src/tallies/tally_filter_cellborn.F90 index 450858b6b..0d7461ef8 100644 --- a/src/tallies/tally_filter_cellborn.F90 +++ b/src/tallies/tally_filter_cellborn.F90 @@ -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 diff --git a/src/tallies/tally_filter_cellfrom.F90 b/src/tallies/tally_filter_cellfrom.F90 index 16cb294d5..2c5402122 100644 --- a/src/tallies/tally_filter_cellfrom.F90 +++ b/src/tallies/tally_filter_cellfrom.F90 @@ -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 diff --git a/src/tallies/tally_filter_delayedgroup.F90 b/src/tallies/tally_filter_delayedgroup.F90 index 32010e815..1a662005b 100644 --- a/src/tallies/tally_filter_delayedgroup.F90 +++ b/src/tallies/tally_filter_delayedgroup.F90 @@ -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 diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index 3cb34efe4..6c42161ca 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -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 diff --git a/src/tallies/tally_filter_energy.F90 b/src/tallies/tally_filter_energy.F90 index d186fa62a..93da69edd 100644 --- a/src/tallies/tally_filter_energy.F90 +++ b/src/tallies/tally_filter_energy.F90 @@ -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 diff --git a/src/tallies/tally_filter_energyfunc.F90 b/src/tallies/tally_filter_energyfunc.F90 index efaabae2e..1402707b4 100644 --- a/src/tallies/tally_filter_energyfunc.F90 +++ b/src/tallies/tally_filter_energyfunc.F90 @@ -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 diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index 93d050b1a..e57423ddf 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -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 diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 index 9545b8692..b4ee7b06b 100644 --- a/src/tallies/tally_filter_legendre.F90 +++ b/src/tallies/tally_filter_legendre.F90 @@ -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 diff --git a/src/tallies/tally_filter_material.F90 b/src/tallies/tally_filter_material.F90 index 3b9f843f8..443dc6285 100644 --- a/src/tallies/tally_filter_material.F90 +++ b/src/tallies/tally_filter_material.F90 @@ -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 diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index d2acdc4c8..486569da6 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -2,8 +2,6 @@ module tally_filter_mesh use, intrinsic :: ISO_C_BINDING - use hdf5 - use constants use dict_header, only: EMPTY use error diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index d90b63664..801d4252c 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -2,8 +2,6 @@ module tally_filter_meshsurface use, intrinsic :: ISO_C_BINDING - use hdf5 - use constants use dict_header, only: EMPTY use error diff --git a/src/tallies/tally_filter_mu.F90 b/src/tallies/tally_filter_mu.F90 index c949289b6..7f8bde923 100644 --- a/src/tallies/tally_filter_mu.F90 +++ b/src/tallies/tally_filter_mu.F90 @@ -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 diff --git a/src/tallies/tally_filter_polar.F90 b/src/tallies/tally_filter_polar.F90 index 89815c102..7cbb55f8b 100644 --- a/src/tallies/tally_filter_polar.F90 +++ b/src/tallies/tally_filter_polar.F90 @@ -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 diff --git a/src/tallies/tally_filter_sph_harm.F90 b/src/tallies/tally_filter_sph_harm.F90 index 5faf2b6fd..90f04c488 100644 --- a/src/tallies/tally_filter_sph_harm.F90 +++ b/src/tallies/tally_filter_sph_harm.F90 @@ -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 diff --git a/src/tallies/tally_filter_sptl_legendre.F90 b/src/tallies/tally_filter_sptl_legendre.F90 index d6594d05a..83db9a864 100644 --- a/src/tallies/tally_filter_sptl_legendre.F90 +++ b/src/tallies/tally_filter_sptl_legendre.F90 @@ -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 diff --git a/src/tallies/tally_filter_surface.F90 b/src/tallies/tally_filter_surface.F90 index df3bf3603..ec642bc8d 100644 --- a/src/tallies/tally_filter_surface.F90 +++ b/src/tallies/tally_filter_surface.F90 @@ -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 diff --git a/src/tallies/tally_filter_universe.F90 b/src/tallies/tally_filter_universe.F90 index 742152232..0dc5aefe1 100644 --- a/src/tallies/tally_filter_universe.F90 +++ b/src/tallies/tally_filter_universe.F90 @@ -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 diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index e96a53a7f..fa205b180 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -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 diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 413464b7b..cbfa4f048 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -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 !=============================================================================== diff --git a/src/track_output.F90 b/src/track_output.F90 index 244bf182e..95af9e6b0 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -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) diff --git a/src/urr_header.F90 b/src/urr_header.F90 index cc41d43cc..0ec9e89de 100644 --- a/src/urr_header.F90 +++ b/src/urr_header.F90 @@ -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 diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index e374f2106..8264797a7 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -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") diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp new file mode 100644 index 000000000..aea1c2c17 --- /dev/null +++ b/src/xml_interface.cpp @@ -0,0 +1,40 @@ +#include "xml_interface.h" + +#include // for std::transform +#include +#include + +#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 diff --git a/src/xml_interface.h b/src/xml_interface.h index 778497562..52aa6003f 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -1,8 +1,6 @@ #ifndef XML_INTERFACE_H #define XML_INTERFACE_H -#include // for std::transform -#include #include #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