From 65095096f9ad72981f1be9c6ad9d89e55e77ea12 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 20 Apr 2018 13:53:53 -0500 Subject: [PATCH] Convert read_source_bank to C++ and fix a bunch of bugs --- src/hdf5_interface.F90 | 18 +++++------ src/hdf5_interface.cpp | 7 +++- src/initialize.cpp | 40 +++++++++++++++-------- src/initialize.h | 14 ++++++++ src/main.cpp | 4 +++ src/message_passing.h | 2 ++ src/mgxs_data.F90 | 4 +++ src/mgxs_header.F90 | 20 ++++++++++++ src/source.F90 | 2 +- src/state_point.F90 | 73 +++++------------------------------------- src/state_point.cpp | 71 ++++++++++++++++++++++++++++++++++------ src/state_point.h | 4 ++- 12 files changed, 159 insertions(+), 100 deletions(-) create mode 100644 src/initialize.h diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 752ce062dc..d76f7b6b71 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -138,7 +138,7 @@ module hdf5_interface integer(HID_T), value :: obj_id type(C_PTR), value :: name real(C_DOUBLE), intent(out) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine read_double_c subroutine read_attr_int_c(obj_id, name, buffer) & @@ -172,7 +172,7 @@ module hdf5_interface integer(HID_T), value :: obj_id type(C_PTR), value :: name integer(C_INT), intent(out) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine read_int_c subroutine read_llong_c(obj_id, name, buffer, indep) & @@ -181,7 +181,7 @@ module hdf5_interface integer(HID_T), value :: obj_id type(C_PTR), value :: name integer(C_LONG_LONG), intent(out) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine read_llong_c subroutine read_string_c(obj_id, name, slen, buffer, indep) & @@ -191,7 +191,7 @@ module hdf5_interface type(C_PTR), value :: name integer(C_SIZE_T), value :: slen character(kind=C_CHAR), intent(out) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine read_string_c subroutine read_complex_c(obj_id, name, buffer, indep) & @@ -200,7 +200,7 @@ module hdf5_interface integer(HID_T), value :: obj_id type(C_PTR), value :: name complex(C_DOUBLE_COMPLEX), intent(out) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine read_complex_c subroutine write_attr_double_c(obj_id, ndim, dims, name, buffer) & @@ -239,7 +239,7 @@ module hdf5_interface integer(HSIZE_T), intent(in) :: dims(*) character(kind=C_CHAR), intent(in) :: name(*) real(C_DOUBLE), intent(in) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine write_double_c subroutine write_int_c(group_id, ndim, dims, name, buffer, indep) & @@ -250,7 +250,7 @@ module hdf5_interface integer(HSIZE_T), intent(in) :: dims(*) character(kind=C_CHAR), intent(in) :: name(*) integer(C_INT), intent(in) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine write_int_c subroutine write_llong_c(group_id, ndim, dims, name, buffer, indep) & @@ -261,7 +261,7 @@ module hdf5_interface integer(HSIZE_T), intent(in) :: dims(*) character(kind=C_CHAR), intent(in) :: name(*) integer(C_LONG_LONG), intent(in) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine write_llong_c subroutine write_string_c(group_id, ndim, dims, slen, name, buffer, indep) & @@ -273,7 +273,7 @@ module hdf5_interface integer(C_SIZE_T), value :: slen character(kind=C_CHAR), intent(in) :: name(*) character(kind=C_CHAR), intent(in) :: buffer(*) - logical(C_BOOL), intent(in) :: indep + logical(C_BOOL), value :: indep end subroutine write_string_c end interface diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 1502ec2e12..d81c5b34a8 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -174,6 +174,11 @@ file_open(const char* filename, char mode, bool parallel) } 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 @@ -451,7 +456,7 @@ write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name, 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}; + auto data_xfer_mode = indep ? H5FD_MPIO_INDEPENDENT : H5FD_MPIO_COLLECTIVE; // Create dataset transfer property list hid_t plist = H5Pcreate(H5P_DATASET_XFER); diff --git a/src/initialize.cpp b/src/initialize.cpp index d32a6c97f3..23025feb23 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -1,17 +1,27 @@ -#include "openmc.h" +#include "initialize.h" + +#include #include "message_passing.h" +#include "openmc.h" int openmc_init(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++ - MPI_Comm intracomm = *static_cast(intracomm); - initialize_mpi(intracomm); + openmc::initialize_mpi(comm); // Continue with rest of initialization - MPI_Fint fcomm = MPI_Comm_c2f(openmc::mpi::intracomm); + MPI_Fint fcomm = MPI_Comm_c2f(comm); openmc_init_f(&fcomm); #else openmc_init_f(nullptr); @@ -19,6 +29,7 @@ int openmc_init(const void* intracomm) return 0; } +namespace openmc { #ifdef OPENMC_MPI void initialize_mpi(MPI_Comm intracomm) @@ -30,8 +41,8 @@ void initialize_mpi(MPI_Comm intracomm) if (!MPI_Initialized(&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); + 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; @@ -40,15 +51,18 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype Bank b; - MPI_Aint disp[5]; - disp[0] = &b.wgt - &b; - disp[1] = &b.xyz - &b; - disp[2] = &b.uvw - &b; - disp[3] = &b.E - &b; - disp[4] = &b.delayed_group - &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 + +#endif // OPENMC_MPI +} // namespace openmc diff --git a/src/initialize.h b/src/initialize.h new file mode 100644 index 0000000000..430987038d --- /dev/null +++ b/src/initialize.h @@ -0,0 +1,14 @@ +#ifndef INITIALIZE_H +#define INITIALIZE_H + +#include "mpi.h" + +namespace openmc { + +#ifdef OPENMC_MPI + void initialize_mpi(MPI_Comm intracomm); +#endif + +} + +#endif // INITIALIZE_H diff --git a/src/main.cpp b/src/main.cpp index 2814b538f7..1dd20ebe04 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,5 +1,9 @@ +#ifdef OPENMC_MPI +#include "mpi.h" +#endif #include "openmc.h" + int main(int argc, char** argv) { int err; diff --git a/src/message_passing.h b/src/message_passing.h index fea092c618..67353e23a1 100644 --- a/src/message_passing.h +++ b/src/message_passing.h @@ -1,6 +1,8 @@ #ifndef MESSAGE_PASSING_H #define MESSAGE_PASSING_H +#include "mpi.h" + namespace openmc { namespace mpi { diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 1b14e6e9f0..e0631d16c5 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 e599a89ff3..76042703af 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -549,6 +549,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 @@ -679,6 +681,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 @@ -782,6 +786,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 @@ -838,6 +844,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 @@ -963,6 +971,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 @@ -1346,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 @@ -1517,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 @@ -1651,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 @@ -1720,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 @@ -1868,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/source.F90 b/src/source.F90 index 59a117465c..b02430823b 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -63,7 +63,7 @@ contains 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) diff --git a/src/state_point.F90 b/src/state_point.F90 index 181ab968b9..5b1c4d462a 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -13,8 +13,6 @@ module state_point use, intrinsic :: ISO_C_BINDING - use hdf5 - use bank_header, only: Bank use cmfd_header use constants @@ -45,6 +43,13 @@ module state_point 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 @@ -832,7 +837,7 @@ contains end if ! Write out source - call read_source_bank(file_id) + call read_source_bank(file_id, work_index, source_bank) end if @@ -841,66 +846,4 @@ contains end subroutine load_state_point -!=============================================================================== -! READ_SOURCE_BANK reads OpenMC source_bank data -!=============================================================================== - - subroutine read_source_bank(group_id) - 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 index 20338c4be5..cc6aa2a41d 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -4,14 +4,14 @@ #include #include "mpi.h" +#include "error.h" #include "message_passing.h" #include "openmc.h" namespace openmc { -void -write_source_bank(hid_t group_id, int64_t* work_index, const Bank* source_bank) -{ + +hid_t h5banktype() { // Create type for array of 3 reals hsize_t dims[] {3}; hid_t triplet = H5Tarray_create(H5T_NATIVE_DOUBLE, 1, dims); @@ -24,31 +24,41 @@ write_source_bank(hid_t group_id, int64_t* work_index, const Bank* source_bank) 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, const Bank* source_bank) +{ + hid_t banktype = h5banktype(); + #ifdef PHDF5 // Set size of total dataspace for all procs and rank - dims[0] = n_particles; + hsize_t dims[] {static_cast(n_particles)}; hid_t dspace = H5Screate_simple(1, dims, H5P_DEFAULT); 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[] {openmc_work}; + hsize_t count[] {static_cast(openmc_work)}; hid_t memspace = H5Screate_simple(1, count, H5P_DEFAULT); // Select hyperslab for this dataspace - hsize_t start[] {work_index[openmc::mpi::rank]}; + 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_f(plist, H5FD_MPIO_COLLECTIVE); + H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); // Write data to file in parallel - H5Dwrite(dset, banktype, memspace, dspace, memspace, plist, source_bank); + H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank); // Free resources H5Sclose(dspace); - h5sclose(memspace); + H5Sclose(memspace); H5Dclose(dset); H5Pclose(plist); @@ -106,7 +116,48 @@ write_source_bank(hid_t group_id, int64_t* work_index, const Bank* source_bank) #endif H5Tclose(banktype); - H5Tclose(triplet); +} + + +void read_source_bank(hid_t group_id, int64_t* work_index, struct 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, H5P_DEFAULT); + + // 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 index f0dafdb809..34186d0d97 100644 --- a/src/state_point.h +++ b/src/state_point.h @@ -9,7 +9,9 @@ namespace openmc { extern "C" void write_source_bank(hid_t group_id, int64_t* work_index, - const struct Bank* bank); + const struct Bank* source_bank); +extern "C" void read_source_bank(hid_t group_id, int64_t* work_index, + struct Bank* source_bank); } // namespace openmc #endif // STATE_POINT_H