From 57991b271db4828e8fbc92fbdf78b5892f83fb7d Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 19 Jan 2018 18:12:57 -0500 Subject: [PATCH] Start moving surface->summary.h5 to C++ --- CMakeLists.txt | 5 +- src/api.F90 | 2 +- src/cmfd_execute.F90 | 8 +- src/eigenvalue.F90 | 18 ++-- src/error.F90 | 4 +- src/geometry.F90 | 1 - src/hdf5_interface.F90 | 4 +- src/hdf5_interface.h | 47 ++++++++++ src/initialize.F90 | 12 +-- src/input_xml.F90 | 80 +++-------------- src/main.F90 | 8 +- src/mesh.F90 | 4 +- src/message_passing.F90 | 6 +- src/output.F90 | 2 +- src/simulation.F90 | 8 +- src/source.F90 | 2 +- src/state_point.F90 | 18 ++-- src/summary.F90 | 92 +++----------------- src/surface_header.C | 184 +++++++++++++++++++++++++++++++++++++++- src/surface_header.F90 | 45 ---------- src/tallies/tally.F90 | 4 +- src/tallies/trigger.F90 | 2 +- src/volume_calc.F90 | 6 +- 23 files changed, 311 insertions(+), 251 deletions(-) create mode 100644 src/hdf5_interface.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 200d836fe..58fb55e79 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,14 +48,14 @@ add_definitions(-DMAX_COORD=${maxcoord}) set(MPI_ENABLED FALSE) if($ENV{FC} MATCHES "(mpi[^/]*|ftn)$") message("-- Detected MPI wrapper: $ENV{FC}") - add_definitions(-DMPI) + add_definitions(-DOPENMC_MPI) set(MPI_ENABLED TRUE) endif() # Check for Fortran 2008 MPI interface if(MPI_ENABLED AND mpif08) message("-- Using Fortran 2008 MPI bindings") - add_definitions(-DMPIF08) + add_definitions(-DOPENMC_MPIF08) endif() #=============================================================================== @@ -435,6 +435,7 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/trigger_header.F90 ) set(LIBOPENMC_CXX_SRC + src/hdf5_interface.h src/random_lcg.h src/random_lcg.cpp src/surface_header.C diff --git a/src/api.F90 b/src/api.F90 index eb9a99695..5d56b857c 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -171,7 +171,7 @@ contains ! Close FORTRAN interface. call h5close_f(err) -#ifdef MPI +#ifdef OPENMC_MPI ! Free all MPI types call MPI_TYPE_FREE(MPI_BANK, err) #endif diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 9cdb751a4..af8d57a80 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -105,7 +105,7 @@ contains real(8) :: hxyz(3) ! cell dimensions of current ijk cell real(8) :: vol ! volume of cell real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif @@ -197,7 +197,7 @@ contains end if -#ifdef MPI +#ifdef OPENMC_MPI ! Broadcast full source to all procs call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) #endif @@ -234,7 +234,7 @@ contains real(8) :: norm ! normalization factor logical :: outside ! any source sites outside mesh logical :: in_mesh ! source site is inside mesh -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err #endif @@ -291,7 +291,7 @@ contains if (.not. cmfd_feedback) return ! Broadcast weight factors to all procs -#ifdef MPI +#ifdef OPENMC_MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & mpi_intracomm, mpi_err) #endif diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 27ae4063e..7a11d1fb4 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -42,11 +42,11 @@ contains type(Bank), save, allocatable :: & & temp_sites(:) ! local array of extra sites on each node -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code integer(8) :: n ! number of sites to send/recv integer :: neighbor ! processor to send/recv data from -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Request) :: request(20) #else integer :: request(20) ! communication request for send/recving sites @@ -66,7 +66,7 @@ contains ! fission bank its own sites starts in order to ensure reproducibility by ! skipping ahead to the proper seed. -#ifdef MPI +#ifdef OPENMC_MPI start = 0_8 call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & mpi_intracomm, mpi_err) @@ -148,7 +148,7 @@ contains ! neighboring processors, we have to perform an ALLGATHER to determine the ! indices for all processors -#ifdef MPI +#ifdef OPENMC_MPI ! First do an exclusive scan to get the starting indices for start = 0_8 call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, & @@ -191,7 +191,7 @@ contains call time_bank_sample % stop() call time_bank_sendrecv % start() -#ifdef MPI +#ifdef OPENMC_MPI ! ========================================================================== ! SEND BANK SITES TO NEIGHBORS @@ -343,14 +343,14 @@ contains subroutine calculate_generation_keff() real(8) :: keff_reduced -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif ! Get keff for this generation by subtracting off the starting value keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation -#ifdef MPI +#ifdef OPENMC_MPI ! Combine values across all processors call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, & MPI_SUM, mpi_intracomm, mpi_err) @@ -584,7 +584,7 @@ contains real(8) :: total ! total weight in source bank logical :: sites_outside ! were there sites outside the ufs mesh? -#ifdef MPI +#ifdef OPENMC_MPI integer :: n ! total number of ufs mesh cells integer :: mpi_err ! MPI error code #endif @@ -608,7 +608,7 @@ contains call fatal_error("Source sites outside of the UFS mesh!") end if -#ifdef MPI +#ifdef OPENMC_MPI ! Send source fraction to all processors n = product(m % dimension) call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) diff --git a/src/error.F90 b/src/error.F90 index 27004d729..ae02e0897 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -128,7 +128,7 @@ contains integer :: line_wrap ! length of line integer :: length ! length of message integer :: indent ! length of indentation -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err #endif @@ -180,7 +180,7 @@ contains end if end do -#ifdef MPI +#ifdef OPENMC_MPI ! Abort MPI call MPI_ABORT(mpi_intracomm, code, mpi_err) #endif diff --git a/src/geometry.F90 b/src/geometry.F90 index c19a6c8dd..f1feb3ba3 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -537,7 +537,6 @@ contains real(8) :: surf_uvw(3) ! surface normal direction logical :: coincident ! is particle on surface? type(Cell), pointer :: c - class(Surface), pointer :: surf class(Lattice), pointer :: lat ! inialize distance to infinity (huge) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index e8c39b923..410149d9e 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -124,7 +124,7 @@ contains ! Setup file access property list with parallel I/O access call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) #ifdef PHDF5 -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & MPI_INFO_NULL%MPI_VAL, hdf5_err) #else @@ -174,7 +174,7 @@ contains ! Setup file access property list with parallel I/O access call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) #ifdef PHDF5 -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & MPI_INFO_NULL%MPI_VAL, hdf5_err) #else diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h new file mode 100644 index 000000000..38bb7d88b --- /dev/null +++ b/src/hdf5_interface.h @@ -0,0 +1,47 @@ +#ifndef HDF5_INTERFACE_H +#define HDF5_INTERFACE_H + +#include // For std::array +#include // For strlen + +#include "hdf5.h" + + +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); +} + +#endif //HDF5_INTERFACE_H diff --git a/src/initialize.F90 b/src/initialize.F90 index dd7894f96..2af93627e 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -52,8 +52,8 @@ contains ! Copy the communicator to a new variable. This is done to avoid changing ! the signature of this subroutine. If MPI is being used but no communicator ! was passed, assume MPI_COMM_WORLD. -#ifdef MPI -#ifdef MPIF08 +#ifdef OPENMC_MPI +#ifdef OPENMC_MPIF08 type(MPI_Comm), intent(in) :: comm ! MPI intracommunicator if (present(intracomm)) then comm % MPI_VAL = intracomm @@ -74,7 +74,7 @@ contains call time_total%start() call time_initialize%start() -#ifdef MPI +#ifdef OPENMC_MPI ! Setup MPI call initialize_mpi(comm) #endif @@ -108,7 +108,7 @@ contains end subroutine openmc_init -#ifdef MPI +#ifdef OPENMC_MPI !=============================================================================== ! INITIALIZE_MPI starts up the Message Passing Interface (MPI) and determines ! the number of processors the problem is being run with as well as the rank of @@ -116,7 +116,7 @@ contains !=============================================================================== subroutine initialize_mpi(intracomm) -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator #else integer, intent(in) :: intracomm ! MPI intracommunicator @@ -124,7 +124,7 @@ contains integer :: mpi_err ! MPI error code integer :: bank_blocks(5) ! Count for each datatype -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Datatype) :: bank_types(5) #else integer :: bank_types(5) ! Datatypes diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4994ff76b..462b9708f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -54,7 +54,7 @@ module input_xml implicit none type(C_PTR) :: node_ptr end subroutine read_surfaces - end interface + end interface contains @@ -1081,77 +1081,23 @@ contains select type(s) type is (SurfaceXPlane) - s%x0 = coeffs(1) - ! Determine outer surfaces - xmin = min(xmin, s % x0) - xmax = max(xmax, s % x0) - if (xmin == s % x0) i_xmin = i - if (xmax == s % x0) i_xmax = i + xmin = min(xmin, coeffs(1)) + xmax = max(xmax, coeffs(1)) + if (xmin == coeffs(1)) i_xmin = i + if (xmax == coeffs(1)) i_xmax = i type is (SurfaceYPlane) - s%y0 = coeffs(1) - ! Determine outer surfaces - ymin = min(ymin, s % y0) - ymax = max(ymax, s % y0) - if (ymin == s % y0) i_ymin = i - if (ymax == s % y0) i_ymax = i + ymin = min(ymin, coeffs(1)) + ymax = max(ymax, coeffs(1)) + if (ymin == coeffs(1)) i_ymin = i + if (ymax == coeffs(1)) i_ymax = i type is (SurfaceZPlane) - s%z0 = coeffs(1) - ! Determine outer surfaces - zmin = min(zmin, s % z0) - zmax = max(zmax, s % z0) - if (zmin == s % z0) i_zmin = i - if (zmax == s % z0) i_zmax = i - type is (SurfacePlane) - s%A = coeffs(1) - s%B = coeffs(2) - s%C = coeffs(3) - s%D = coeffs(4) - type is (SurfaceXCylinder) - s%y0 = coeffs(1) - s%z0 = coeffs(2) - s%r = coeffs(3) - type is (SurfaceYCylinder) - s%x0 = coeffs(1) - s%z0 = coeffs(2) - s%r = coeffs(3) - type is (SurfaceZCylinder) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%r = coeffs(3) - type is (SurfaceSphere) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r = coeffs(4) - type is (SurfaceXCone) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r2 = coeffs(4) - type is (SurfaceYCone) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r2 = coeffs(4) - type is (SurfaceZCone) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r2 = coeffs(4) - type is (SurfaceQuadric) - s%A = coeffs(1) - s%B = coeffs(2) - s%C = coeffs(3) - s%D = coeffs(4) - s%E = coeffs(5) - s%F = coeffs(6) - s%G = coeffs(7) - s%H = coeffs(8) - s%J = coeffs(9) - s%K = coeffs(10) + zmin = min(zmin, coeffs(1)) + zmax = max(zmax, coeffs(1)) + if (zmin == coeffs(1)) i_zmin = i + if (zmax == coeffs(1)) i_zmax = i end select ! No longer need coefficients diff --git a/src/main.F90 b/src/main.F90 index d8e52643c..120bdd1e0 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -9,13 +9,13 @@ program main implicit none -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif ! Initialize run -- when run with MPI, pass communicator -#ifdef MPI -#ifdef MPIF08 +#ifdef OPENMC_MPI +#ifdef OPENMC_MPIF08 call openmc_init(MPI_COMM_WORLD % MPI_VAL) #else call openmc_init(MPI_COMM_WORLD) @@ -39,7 +39,7 @@ program main ! finalize run call openmc_finalize() -#ifdef MPI +#ifdef OPENMC_MPI ! If MPI is in use and enabled, terminate it call MPI_FINALIZE(mpi_err) #endif diff --git a/src/mesh.F90 b/src/mesh.F90 index e997890a4..790dc7d88 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -34,7 +34,7 @@ contains integer :: n ! number of energy groups / size integer :: mesh_bin ! mesh bin integer :: e_bin ! energy bin -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif logical :: outside ! was any site outside mesh? @@ -86,7 +86,7 @@ contains cnt_(e_bin, mesh_bin) = cnt_(e_bin, mesh_bin) + bank_array(i) % wgt end do FISSION_SITES -#ifdef MPI +#ifdef OPENMC_MPI ! collect values from all processors n = size(cnt_) call MPI_REDUCE(cnt_, cnt, n, MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err) diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 0f2b94d1a..7391a632c 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -1,7 +1,7 @@ module message_passing -#ifdef MPI -#ifdef MPIF08 +#ifdef OPENMC_MPI +#ifdef OPENMC_MPIF08 use mpi_f08 #else use mpi @@ -16,7 +16,7 @@ module message_passing integer :: rank = 0 ! rank of process logical :: master = .true. ! master process? logical :: mpi_enabled = .false. ! is MPI in use and initialized? -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator #else diff --git a/src/output.F90 b/src/output.F90 index 9a144721d..1b7cae121 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -88,7 +88,7 @@ contains ! Write the date and time write(UNIT=OUTPUT_UNIT, FMT='(9X,"Date/Time | ",A)') time_stamp() -#ifdef MPI +#ifdef OPENMC_MPI ! Write number of processors write(UNIT=OUTPUT_UNIT, FMT='(5X,"MPI Processes | ",A)') & trim(to_str(n_procs)) diff --git a/src/simulation.F90 b/src/simulation.F90 index ef29e7832..24e9b29fc 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -320,7 +320,7 @@ contains subroutine finalize_batch() -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif @@ -342,7 +342,7 @@ contains ! Check_triggers if (master) call check_triggers() -#ifdef MPI +#ifdef OPENMC_MPI call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, & mpi_intracomm, mpi_err) #endif @@ -469,7 +469,7 @@ contains subroutine openmc_simulation_finalize() bind(C) integer :: i ! loop index -#ifdef MPI +#ifdef OPENMC_MPI integer :: n ! size of arrays integer :: mpi_err ! MPI error code integer(8) :: temp @@ -494,7 +494,7 @@ contains ! Increment total number of generations total_gen = total_gen + current_batch*gen_per_batch -#ifdef MPI +#ifdef OPENMC_MPI ! Broadcast tally results so that each process has access to results if (allocated(tallies)) then do i = 1, size(tallies) diff --git a/src/source.F90 b/src/source.F90 index 557120e97..5beecd887 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,7 +1,7 @@ module source use hdf5, only: HID_T -#ifdef MPI +#ifdef OPENMC_MPI use message_passing #endif diff --git a/src/state_point.F90 b/src/state_point.F90 index f7a892898..ba7bcef19 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -523,7 +523,7 @@ contains integer(HID_T) :: tallies_group, tally_group real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results real(8), target :: global_temp(3,N_GLOBAL_TALLIES) -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code real(8) :: dummy ! temporary receive buffer for non-root reduces #endif @@ -543,7 +543,7 @@ contains end if -#ifdef MPI +#ifdef OPENMC_MPI ! Reduce global tallies n_bins = size(global_tallies) call MPI_REDUCE(global_tallies, global_temp, n_bins, MPI_REAL8, MPI_SUM, & @@ -586,7 +586,7 @@ contains ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable -#ifdef MPI +#ifdef OPENMC_MPI call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & MPI_SUM, 0, mpi_intracomm, mpi_err) #endif @@ -610,7 +610,7 @@ contains deallocate(dummy_tally % results) else ! Receive buffer not significant at other processors -#ifdef MPI +#ifdef OPENMC_MPI call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & 0, mpi_intracomm, mpi_err) #endif @@ -849,7 +849,7 @@ contains integer(HID_T) :: plist ! property list #else integer :: i -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code type(Bank), allocatable, target :: temp_source(:) #endif @@ -897,7 +897,7 @@ contains dspace, dset, hdf5_err) ! Save source bank sites since the souce_bank array is overwritten below -#ifdef MPI +#ifdef OPENMC_MPI allocate(temp_source(work)) temp_source(:) = source_bank(:) #endif @@ -907,7 +907,7 @@ contains dims(1) = work_index(i+1) - work_index(i) call h5screate_simple_f(1, dims, memspace, hdf5_err) -#ifdef MPI +#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, & @@ -933,12 +933,12 @@ contains call h5dclose_f(dset, hdf5_err) ! Restore state of source bank -#ifdef MPI +#ifdef OPENMC_MPI source_bank(:) = temp_source(:) deallocate(temp_source) #endif else -#ifdef MPI +#ifdef OPENMC_MPI call MPI_SEND(source_bank, int(work), MPI_BANK, 0, rank, & mpi_intracomm, mpi_err) #endif diff --git a/src/summary.F90 b/src/summary.F90 index 0f45ef1bb..ca7276407 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -24,6 +24,17 @@ module summary public :: write_summary + interface + subroutine surface_to_hdf5_c(surf_ind, group) & + bind(C, name='surface_to_hdf5') + use ISO_C_BINDING + use hdf5 + implicit none + integer(C_INT), intent(in), value :: surf_ind + integer(HID_T), intent(in), value :: group + end subroutine surface_to_hdf5_c + end interface + contains !=============================================================================== @@ -125,7 +136,6 @@ contains integer(HID_T) :: surfaces_group, surface_group integer(HID_T) :: universes_group, univ_group integer(HID_T) :: lattices_group, lattice_group - real(8), allocatable :: coeffs(:) character(:), allocatable :: region_spec type(Cell), pointer :: c class(Surface), pointer :: s @@ -245,87 +255,11 @@ contains surface_group = create_group(surfaces_group, "surface " // & trim(to_str(s%id))) + call surface_to_hdf5_c(i-1, surface_group) + ! Write name for this surface call write_dataset(surface_group, "name", s%name) - ! Write surface type - select type (s) - type is (SurfaceXPlane) - call write_dataset(surface_group, "type", "x-plane") - allocate(coeffs(1)) - coeffs(1) = s%x0 - - type is (SurfaceYPlane) - call write_dataset(surface_group, "type", "y-plane") - allocate(coeffs(1)) - coeffs(1) = s%y0 - - type is (SurfaceZPlane) - call write_dataset(surface_group, "type", "z-plane") - allocate(coeffs(1)) - coeffs(1) = s%z0 - - type is (SurfacePlane) - call write_dataset(surface_group, "type", "plane") - allocate(coeffs(4)) - coeffs(:) = [s%A, s%B, s%C, s%D] - - type is (SurfaceXCylinder) - call write_dataset(surface_group, "type", "x-cylinder") - allocate(coeffs(3)) - coeffs(:) = [s%y0, s%z0, s%r] - - type is (SurfaceYCylinder) - call write_dataset(surface_group, "type", "y-cylinder") - allocate(coeffs(3)) - coeffs(:) = [s%x0, s%z0, s%r] - - type is (SurfaceZCylinder) - call write_dataset(surface_group, "type", "z-cylinder") - allocate(coeffs(3)) - coeffs(:) = [s%x0, s%y0, s%r] - - type is (SurfaceSphere) - call write_dataset(surface_group, "type", "sphere") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r] - - type is (SurfaceXCone) - call write_dataset(surface_group, "type", "x-cone") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r2] - - type is (SurfaceYCone) - call write_dataset(surface_group, "type", "y-cone") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r2] - - type is (SurfaceZCone) - call write_dataset(surface_group, "type", "z-cone") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r2] - - type is (SurfaceQuadric) - call write_dataset(surface_group, "type", "quadric") - allocate(coeffs(10)) - coeffs(:) = [s%A, s%B, s%C, s%D, s%E, s%F, s%G, s%H, s%J, s%K] - - end select - call write_dataset(surface_group, "coefficients", coeffs) - deallocate(coeffs) - - ! Write boundary type - select case (s%bc) - case (BC_TRANSMIT) - call write_dataset(surface_group, "boundary_type", "transmission") - case (BC_VACUUM) - call write_dataset(surface_group, "boundary_type", "vacuum") - case (BC_REFLECT) - call write_dataset(surface_group, "boundary_type", "reflective") - case (BC_PERIODIC) - call write_dataset(surface_group, "boundary_type", "periodic") - end select - call close_group(surface_group) end do SURFACE_LOOP diff --git a/src/surface_header.C b/src/surface_header.C index 45fb62399..db53c8982 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -1,8 +1,12 @@ +#include // For std::array #include // For strcmp #include // For numeric_limits #include // For fabs #include "pugixml/pugixml.hpp" +#include "hdf5.h" + +#include "hdf5_interface.h" // DEBUGGING #include @@ -15,6 +19,11 @@ const double FP_COINCIDENT = 1e-12; const double INFTY = std::numeric_limits::max(); +const int BC_TRANSMIT{0}; +const int BC_VACUUM{1}; +const int BC_REFLECT{2}; +const int BC_PERIODIC{3}; + //============================================================================== // Global array of surfaces //============================================================================== @@ -127,6 +136,13 @@ public: //! @param xyz[3] A 3D Cartesian coordinate. //! @param uvw[3] This output argument provides the normal. virtual void normal(const double xyz[3], double uvw[3]) const = 0; + + //! Write all information needed to reconstruct the surface to an HDF5 group. + //! @param group_id An HDF5 group id. + virtual void to_hdf5(hid_t group_id) const = 0; + +protected: + void write_bc_to_hdf5(hid_t group_id) const; }; bool @@ -164,6 +180,25 @@ Surface::reflect(const double xyz[3], double uvw[3]) const uvw[2] -= 2.0 * projection / magnitude * norm[2]; } +void +Surface::write_bc_to_hdf5(hid_t group_id) const +{ + switch(bc) { + case BC_TRANSMIT : + write_string(group_id, "boundary_type", "transmission"); + break; + case BC_VACUUM : + write_string(group_id, "boundary_type", "vacuum"); + break; + case BC_REFLECT : + write_string(group_id, "boundary_type", "reflective"); + break; + case BC_PERIODIC : + write_string(group_id, "boundary_type", "periodic"); + break; + } +} + //============================================================================== //! A `Surface` that supports periodic boundary conditions. //! @@ -236,6 +271,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; }; @@ -261,6 +297,14 @@ inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const axis_aligned_plane_normal<0, 1, 2>(xyz, uvw); } +void SurfaceXPlane::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "x-plane"); + std::array coeffs{{x0}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const { @@ -285,7 +329,7 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], return true; } -}; +} //============================================================================== // SurfaceYPlane @@ -303,6 +347,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; }; @@ -328,6 +373,14 @@ inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const axis_aligned_plane_normal<1, 0, 2>(xyz, uvw); } +void SurfaceYPlane::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "y-plane"); + std::array coeffs{{y0}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const { @@ -352,7 +405,7 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], return true; } -}; +} //============================================================================== // SurfaceZPlane @@ -370,6 +423,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; }; @@ -395,13 +449,21 @@ inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const axis_aligned_plane_normal<2, 0, 1>(xyz, uvw); } +void SurfaceZPlane::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "z-plane"); + std::array coeffs{{z0}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const { // Assume the other plane is aligned along z. Just change the z coord. xyz[2] = z0; return false; -}; +} //============================================================================== // SurfacePlane @@ -419,6 +481,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; }; @@ -457,6 +520,14 @@ SurfacePlane::normal(const double xyz[3], double uvw[3]) const uvw[2] = C; } +void SurfacePlane::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "plane"); + std::array coeffs{{A, B, C, D}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const { @@ -564,6 +635,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) @@ -588,6 +660,15 @@ inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0); } + +void SurfaceXCylinder::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "x-cylinder"); + std::array coeffs{{y0, z0, r}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // SurfaceYCylinder //! A cylinder aligned along the y-axis. @@ -605,6 +686,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) @@ -629,6 +711,14 @@ inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0); } +void SurfaceYCylinder::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "y-cylinder"); + std::array coeffs{{x0, z0, r}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // SurfaceZCylinder //! A cylinder aligned along the z-axis. @@ -646,6 +736,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) @@ -670,6 +761,14 @@ inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0); } +void SurfaceZCylinder::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "z-cylinder"); + std::array coeffs{{x0, y0, r}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // SurfaceSphere //! A sphere. @@ -687,6 +786,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) @@ -748,6 +848,14 @@ inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const uvw[2] = 2.0 * (xyz[2] - z0); } +void SurfaceSphere::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "sphere"); + std::array coeffs{{x0, y0, z0, r}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // Generic functions for x-, y-, and z-, cones //============================================================================== @@ -848,6 +956,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) @@ -872,6 +981,14 @@ inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const axis_aligned_cone_normal<0, 1, 2>(xyz, uvw, x0, y0, z0, r_sq); } +void SurfaceXCone::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "x-cone"); + std::array coeffs{{x0, y0, z0, r_sq}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // SurfaceYCone //! A cone aligned along the y-axis. @@ -889,6 +1006,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) @@ -913,6 +1031,14 @@ inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const axis_aligned_cone_normal<1, 0, 2>(xyz, uvw, y0, x0, z0, r_sq); } +void SurfaceYCone::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "y-cone"); + std::array coeffs{{x0, y0, z0, r_sq}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // SurfaceZCone //! A cone aligned along the z-axis. @@ -930,6 +1056,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) @@ -954,6 +1081,14 @@ inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const axis_aligned_cone_normal<2, 0, 1>(xyz, uvw, z0, x0, y0, r_sq); } +void SurfaceZCone::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "z-cone"); + std::array coeffs{{x0, y0, z0, r_sq}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== // SurfaceQuadric //! A general surface described by a quadratic equation. @@ -971,6 +1106,7 @@ public: double distance(const double xyz[3], const double uvw[3], bool coincident) const; void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5(hid_t group_id) const; }; SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) @@ -1055,6 +1191,14 @@ SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const uvw[2] = 2.0*C*z + E*y + F*x + J; } +void SurfaceQuadric::to_hdf5(hid_t group_id) const +{ + write_string(group_id, "type", "quadric"); + std::array coeffs{{A, B, C, D, E, F, G, H, J, K}}; + write_double_1D(group_id, "coefficients", coeffs); + write_bc_to_hdf5(group_id); +} + //============================================================================== extern "C" void @@ -1084,6 +1228,7 @@ read_surfaces(pugi::xml_node *node) } else { std::cout << "ERROR: Found a surface with no type attribute/node" << std::endl; + //TODO: call fatal_error } if (strcmp(surf_type, "x-plane") == 0) { @@ -1125,6 +1270,33 @@ read_surfaces(pugi::xml_node *node) } else { std::cout << "Call error or handle uppercase here!" << std::endl; std::cout << surf_type << std::endl; + //TODO: call fatal_error + } + + const pugi::char_t *surf_bc{""}; + if (surf_node.attribute("boundary")) { + surf_bc = surf_node.attribute("boundary").value(); + } else if (surf_node.child("boundary")) { + surf_bc = surf_node.attribute("boundary").value(); + } + + if (strcmp(surf_bc, "transmission") == 0 + or strcmp(surf_bc, "transmit") == 0 + or strcmp(surf_bc, "") == 0) { + surfaces_c[i_surf]->bc = BC_TRANSMIT; + + } else if (strcmp(surf_bc, "vacuum") == 0) { + surfaces_c[i_surf]->bc = BC_VACUUM; + + } else if (strcmp(surf_bc, "reflective") == 0 + or strcmp(surf_bc, "reflect") == 0 + or strcmp(surf_bc, "reflecting") == 0) { + surfaces_c[i_surf]->bc = BC_REFLECT; + } else if (strcmp(surf_bc, "periodic") == 0) { + surfaces_c[i_surf]->bc = BC_PERIODIC; + } else { + std::cout << "Unknown boundary condition" << std::endl; + //TODO: call fatal_error } } } @@ -1156,6 +1328,12 @@ surface_normal(int surf_ind, double xyz[3], double uvw[3]) return surfaces_c[surf_ind]->normal(xyz, uvw); } +extern "C" void +surface_to_hdf5(int surf_ind, hid_t group) +{ + surfaces_c[surf_ind]->to_hdf5(group); +} + extern "C" bool surface_periodic(int surf_ind1, int surf_ind2, double xyz[3], double uvw[3]) { diff --git a/src/surface_header.F90 b/src/surface_header.F90 index d63290387..98e2423b6 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -36,84 +36,39 @@ module surface_header !=============================================================================== type, extends(Surface) :: SurfaceXPlane - ! x = x0 - real(8) :: x0 end type SurfaceXPlane type, extends(Surface) :: SurfaceYPlane - ! y = y0 - real(8) :: y0 end type SurfaceYPlane type, extends(Surface) :: SurfaceZPlane - ! z = z0 - real(8) :: z0 end type SurfaceZPlane type, extends(Surface) :: SurfacePlane - ! Ax + By + Cz = D - real(8) :: A - real(8) :: B - real(8) :: C - real(8) :: D end type SurfacePlane type, extends(Surface) :: SurfaceXCylinder - ! (y - y0)^2 + (z - z0)^2 = R^2 - real(8) :: y0 - real(8) :: z0 - real(8) :: r end type SurfaceXCylinder type, extends(Surface) :: SurfaceYCylinder - ! (x - x0)^2 + (z - z0)^2 = R^2 - real(8) :: x0 - real(8) :: z0 - real(8) :: r end type SurfaceYCylinder type, extends(Surface) :: SurfaceZCylinder - ! (x - x0)^2 + (y - y0)^2 = R^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: r end type SurfaceZCylinder type, extends(Surface) :: SurfaceSphere - ! (x - x0)^2 + (y - y0)^2 + (z - z0)^2 = R^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r end type SurfaceSphere type, extends(Surface) :: SurfaceXCone - ! (y - y0)^2 + (z - z0)^2 = R^2*(x - x0)^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r2 end type SurfaceXCone type, extends(Surface) :: SurfaceYCone - ! (x - x0)^2 + (z - z0)^2 = R^2*(y - y0)^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r2 end type SurfaceYCone type, extends(Surface) :: SurfaceZCone - ! (x - x0)^2 + (y - y0)^2 = R^2*(z - z0)^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r2 end type SurfaceZCone type, extends(Surface) :: SurfaceQuadric - ! Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 - real(8) :: A, B, C, D, E, F, G, H, J, K end type SurfaceQuadric integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 7684a36aa..7c234f173 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -4241,7 +4241,7 @@ contains real(C_DOUBLE) :: k_tra ! Copy of batch tracklength estimate of keff real(C_DOUBLE) :: val -#ifdef MPI +#ifdef OPENMC_MPI ! Combine tally results onto master process if (reduce_tallies) call reduce_tally_results() #endif @@ -4289,7 +4289,7 @@ contains ! REDUCE_TALLY_RESULTS collects all the results from tallies onto one processor !=============================================================================== -#ifdef MPI +#ifdef OPENMC_MPI subroutine reduce_tally_results() integer :: i diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 9e9b06d20..ee2259ed4 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -2,7 +2,7 @@ module trigger use, intrinsic :: ISO_C_BINDING -#ifdef MPI +#ifdef OPENMC_MPI use message_passing #endif diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 5985a236b..07dc5b418 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -136,7 +136,7 @@ contains integer :: total_hits ! total hits for a single domain (summed over materials) integer :: min_samples ! minimum number of samples per process integer :: remainder ! leftover samples from uneven divide -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code integer :: m ! index over materials integer :: n ! number of materials @@ -279,7 +279,7 @@ contains total_hits = 0 if (master) then -#ifdef MPI +#ifdef OPENMC_MPI do j = 1, n_procs - 1 call MPI_RECV(n, 1, MPI_INTEGER, j, 0, mpi_intracomm, & MPI_STATUS_IGNORE, mpi_err) @@ -341,7 +341,7 @@ contains end do else -#ifdef MPI +#ifdef OPENMC_MPI n = master_indices(i_domain) % size() allocate(data(2*n)) do k = 0, n - 1