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

Convert read_source_bank to C++ and fix a bunch of bugs

Browse source
This commit is contained in:
Paul Romano 2018-04-20 13:53:53 -05:00
parent 44ff22dc14
commit 65095096f9
12 changed files with 159 additions and 100 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/hdf5_interface.F90
View file

@ -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

7
src/hdf5_interface.cpp
View file

@ -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);

40
src/initialize.cpp
View file

@ -1,17 +1,27 @@
#include "openmc.h"
#include "initialize.h"
#include <cstddef>
#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<const MPI_Comm *>(intracomm);
} else {
comm = MPI_COMM_WORLD;
}
// Initialize MPI for C++
MPI_Comm intracomm = *static_cast<const MPI_Comm *>(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

14
src/initialize.h Normal file
View file

@ -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

4
src/main.cpp
View file

@ -1,5 +1,9 @@
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "openmc.h"
int main(int argc, char** argv) {
int err;

2
src/message_passing.h
View file

@ -1,6 +1,8 @@
#ifndef MESSAGE_PASSING_H
#define MESSAGE_PASSING_H
#include "mpi.h"
namespace openmc {
namespace mpi {

4
src/mgxs_data.F90
View file

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

20
src/mgxs_header.F90
View file

@ -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

2
src/source.F90
View file

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

73
src/state_point.F90
View file

@ -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

71
src/state_point.cpp
View file

@ -4,14 +4,14 @@
#include <vector>
#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<hsize_t>(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<hsize_t>(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<hsize_t>(work_index[openmc::mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Set up the property list for parallel writing
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio_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<hsize_t>(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<hsize_t>(work_index[openmc::mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, dims, nullptr);
#ifdef PHDF5
// Read data in parallel
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
H5Dread(dset, banktype, memspace, dspace, plist, source_bank);
H5Pclose(plist);
#else
H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
#endif
// Close all ids
H5Sclose(dspace);
H5Sclose(memspace);
H5Dclose(dset);
H5Tclose(banktype);
}
} // namespace openmc

4
src/state_point.h
View file

@ -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